!pip install qiskit

from qiskit.providers.aer import AerSimulator, QasmSimulator
from qiskit.algorithms import VQE
from qiskit.algorithms.optimizers import COBYLA 
from qiskit.circuit.library import TwoLocal
# opflow is Qiskit's module for creating operators like yours
from qiskit import *
from qiskit.opflow import OperatorBase
from qiskit.opflow import Z, X, I  # Pauli Z, X matrices and identity
import pylab
import matplotlib.pyplot as plt
import numpy as np

def store_intermediate_result(eval_count, parameters, mean, std):
  counts.append(eval_count)
  values.append(mean)

def get_gradient(values, counts):
    rvalues = []
    rcounts = []
    for i in range(1, len(values)):
        rvalues += [ (values[i - 1] - values[i]) ]
    rcounts = [ i for i in range(len(rvalues)) ]
    return rvalues, rcounts

def get_absolute(values):
    return [ (v if v > 0 else -v) for v in values]


# J = 1.0
counts = []
values = []

H =   (504.0 * I^I^I^I^I^I^I^Z) + (1008.0 * I^I^I^I^I^I^Z^I) + (2016.0 * I^I^I^I^I^Z^I^I) + (504.0 * I^I^I^I^Z^I^I^I) + (1143.7999999999997 * I^I^I^Z^I^I^I^I) + (2287.6 * I^I^Z^I^I^I^I^I) + (4575.200000000001 * I^Z^I^I^I^I^I^I) + (1143.7999999999997 * Z^I^I^I^I^I^I^I) + (98.0 * I^I^I^I^I^I^Z^Z) + (196.0 * I^I^I^I^I^Z^I^Z) + (392.0 * I^I^I^I^I^Z^Z^I) + (49.0 * I^I^I^I^Z^I^I^Z) + (98.0 * I^I^I^I^Z^I^Z^I) + (196.0 * I^I^I^I^Z^Z^I^I) + (93.1 * I^I^Z^Z^I^I^I^I) + (186.2 * I^Z^I^Z^I^I^I^I) + (372.4 * I^Z^Z^I^I^I^I^I) + (46.55 * Z^I^I^Z^I^I^I^I) + (93.1 * Z^I^Z^I^I^I^I^I) + (186.2 * Z^Z^I^I^I^I^I^I)

# you can swap this for a real quantum device and keep the rest of the code the same!
# backend = AerSimulator(method='matrix_product_state')

backend = QasmSimulator()

# COBYLA usually works well for small problems like this one
optimizer = COBYLA(maxiter=2000)

# EfficientSU2 is a standard heuristic chemistry ansatz from Qiskit's circuit library
ansatz = TwoLocal(num_qubits=8, rotation_blocks=['ry', 'rz'], entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=True)
# set the algorithm
vqe = VQE(ansatz, optimizer, quantum_instance=backend, callback=store_intermediate_result)

#run it with the Hamiltonian we defined above
result = vqe.compute_minimum_eigenvalue(H)  

#print the result (it contains lot's of information)
print(result) 

v1, c1 = get_gradient(values, counts)
plt.figure(figsize=(15,15))
plt.plot(c1,v1)
plt.xlabel('optimization step')
plt.ylabel('E')
plt.show()

counts = []
values = []


H = (504.0 * I^I^I^I^I^I^I^Z) + (1008.0 * I^I^I^I^I^I^Z^I) +  (2016.0 * I^I^I^I^I^Z^I^I) + (504.0 * I^I^I^I^Z^I^I^I) +  (1143.7999999999997 * I^I^I^Z^I^I^I^I) +  (2287.6 * I^I^Z^I^I^I^I^I) +  (4575.200000000001 * I^Z^I^I^I^I^I^I) +  (1143.7999999999997 * Z^I^I^I^I^I^I^I) +  (98.0 * I^I^I^I^I^I^Z^Z) +  (196.0 * I^I^I^I^I^Z^I^Z) +  (392.0 * I^I^I^I^I^Z^Z^I) +  (49.0 * I^I^I^I^Z^I^I^Z) +  (98.0 * I^I^I^I^Z^I^Z^I) +  (196.0 * I^I^I^I^Z^Z^I^I) +  (93.1 * I^I^Z^Z^I^I^I^I) +  (186.2 * I^Z^I^Z^I^I^I^I) +  (372.4 * I^Z^Z^I^I^I^I^I) +  (46.55 * Z^I^I^Z^I^I^I^I) +  (93.1 * Z^I^Z^I^I^I^I^I) +  (186.2 * Z^Z^I^I^I^I^I^I)
# you can swap this for a real quantum device and keep the rest of the code the same!
# backend = AerSimulator(method='matrix_product_state')

backend = QasmSimulator()

# COBYLA usually works well for small problems like this one
optimizer = COBYLA(maxiter=2000)

# EfficientSU2 is a standard heuristic chemistry ansatz from Qiskit's circuit library
ansatz = TwoLocal(num_qubits=8, rotation_blocks=['ry', 'rz'], entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=True)
# set the algorithm
vqe = VQE(ansatz, optimizer, quantum_instance=backend, callback=store_intermediate_result)

#run it with the Hamiltonian we defined above
result = vqe.compute_minimum_eigenvalue(H)  

#print the result (it contains lot's of information)
print(result) 

v2, c2 = get_gradient(values, counts)
plt.figure(figsize=(15,15))
plt.plot(c2,v2)
plt.xlabel('optimization step')
plt.ylabel('E')
plt.show()

counts = []
values = []

H = (504.0 * I^I^I^I^I^I^I^Z) + (1008.0 * I^I^I^I^I^I^Z^I) + (2016.0 * I^I^I^I^I^Z^I^I) + (504.0 * I^I^I^I^Z^I^I^I) + (1143.7999999999997 * I^I^I^Z^I^I^I^I) + (2287.6 * I^I^Z^I^I^I^I^I) + (4575.200000000001 * I^Z^I^I^I^I^I^I) + (1143.7999999999997 * Z^I^I^I^I^I^I^I) + (98.0 * I^I^I^I^I^I^Z^Z)+(196.0 * I^I^I^I^I^Z^I^Z) + (392.0 * I^I^I^I^I^Z^Z^I) +  (49.0 * I^I^I^I^Z^I^I^Z) +  (98.0 * I^I^I^I^Z^I^Z^I) +  (196.0 * I^I^I^I^Z^Z^I^I) +  (93.1 * I^I^Z^Z^I^I^I^I) +  (186.2 * I^Z^I^Z^I^I^I^I)+  (372.4 * I^Z^Z^I^I^I^I^I) +  (46.55 * Z^I^I^Z^I^I^I^I) +  (93.1 * Z^I^Z^I^I^I^I^I) +  (186.2 * Z^Z^I^I^I^I^I^I)

# you can swap this for a real quantum device and keep the rest of the code the same!
# backend = AerSimulator(method='matrix_product_state')

backend = QasmSimulator()

# COBYLA usually works well for small problems like this one
optimizer = COBYLA(maxiter=2000)

# EfficientSU2 is a standard heuristic chemistry ansatz from Qiskit's circuit library
ansatz = TwoLocal(num_qubits=8, rotation_blocks=['ry', 'rz'], entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=True)
# set the algorithm
vqe = VQE(ansatz, optimizer, quantum_instance=backend, callback=store_intermediate_result)

#run it with the Hamiltonian we defined above
result = vqe.compute_minimum_eigenvalue(H)  

#print the result (it contains lot's of information)
print(result) 

v3, c3 = get_gradient(values, counts)
plt.figure(figsize=(15,15))
plt.plot(c3,v3)
plt.xlabel('optimization step')
plt.ylabel('E')
plt.show()

counts = []
values = []

H = (504.0 * I^I^I^I^I^I^I^Z) +  (1008.0 * I^I^I^I^I^I^Z^I) +  (2016.0 * I^I^I^I^I^Z^I^I) +  (504.0 * I^I^I^I^Z^I^I^I) +  (1143.7999999999997 * I^I^I^Z^I^I^I^I) +  (2287.6 * I^I^Z^I^I^I^I^I) +  (4575.200000000001 * I^Z^I^I^I^I^I^I) +  (1143.7999999999997 * Z^I^I^I^I^I^I^I)+  (98.0 * I^I^I^I^I^I^Z^Z) +  (196.0 * I^I^I^I^I^Z^I^Z) +  (392.0 * I^I^I^I^I^Z^Z^I) +  (49.0 * I^I^I^I^Z^I^I^Z)+  (98.0 * I^I^I^I^Z^I^Z^I)+  (196.0 * I^I^I^I^Z^Z^I^I) +  (93.1 * I^I^Z^Z^I^I^I^I)+  (186.2 * I^Z^I^Z^I^I^I^I)+  (372.4 * I^Z^Z^I^I^I^I^I)+  (46.55 * Z^I^I^Z^I^I^I^I)+  (93.1 * Z^I^Z^I^I^I^I^I)+  (186.2 * Z^Z^I^I^I^I^I^I)

# you can swap this for a real quantum device and keep the rest of the code the same!
# backend = AerSimulator(method='matrix_product_state')

backend = QasmSimulator()

# COBYLA usually works well for small problems like this one
optimizer = COBYLA(maxiter=2000)

# EfficientSU2 is a standard heuristic chemistry ansatz from Qiskit's circuit library
ansatz = TwoLocal(num_qubits=8, rotation_blocks=['ry', 'rz'], entanglement_blocks=None, entanglement='full', reps=1, skip_unentangled_qubits=False, skip_final_rotation_layer=True)
# set the algorithm
vqe = VQE(ansatz, optimizer, quantum_instance=backend, callback=store_intermediate_result)

#run it with the Hamiltonian we defined above
result = vqe.compute_minimum_eigenvalue(H)  

#print the result (it contains lot's of information)
print(result) 

v4, c4 = get_gradient(values, counts)
plt.figure(figsize=(15,15))
plt.plot(c4,v4)
plt.xlabel('optimization step')
plt.ylabel('E')
plt.show()

# 2nd Order Der

v1, c1 = get_gradient(v1, c1)
v2, c2 = get_gradient(v2, c2)
v3, c3 = get_gradient(v3, c3)
v4, c4 = get_gradient(v4, c4)


v1 = get_absolute(v1)
v2 = get_absolute(v2)
v3 = get_absolute(v3)
v4 = get_absolute(v4)

plt.figure(figsize=(15,15))
plt.plot(c1,v1, label="R1")
plt.plot(c2,v2, label="R2")
plt.plot(c3,v3, label="R3")
plt.plot(c4,v4, label="R4")
plt.xlabel('comparison optimization step')
plt.ylabel('E')
plt.legend()
plt.show()

plt.figure(figsize=(15,15))
plt.plot(c1[:100],v1[:100], label="R1")
plt.plot(c2[:100],v2[:100], label="R2")
plt.plot(c3[:100],v3[:100], label="R3")
plt.plot(c4[:100],v4[:100], label="R4")
plt.xlabel('comparison optimization step')
plt.ylabel('E')
plt.legend()
plt.show()

