<jupyter_start><jupyter_text>Exercise 5 - Variational quantum eigensolver***By a naive chemist who hated Chemistry at school!******name:*** <NAME> (ironfrown)***score:*** 3***chemical accuracy:*** 2.440695***number of ansatz parameters:*** 16***Note:*** *This notebook version is slightly cleaned up with extra comments, but no functionality or behaviour altered.* ApproachThe approach taken here embraces the following principles.1. As molecular orbitals map onto qubits, we thus need to minimise the number of orbitalsused in molecular modelling (as a Hamiltonian).This was achieved by ***freezing core orbitals*** and by ***removing vacant orbitals***. We also need to investigate ***Z2 symmetries***, which allow reducing the number of qubits when Pauli operators commute with the Hamiltonian (a tapering list identifies such symmetries).2. Once the fermionic model is mapped into an equivalent qubit model, it is then possibleto find its ground state of the minimum energy. This can be achieved using the VQE algorithm, which optimises a parametrised ansatz circuit. The optimisation algorithm identifies a combination of parameter values, which generate the final circuit minimising the electronic energy, thus defining the molecular ground state. For VQE to work effectively, ***all circuit qubits need to be entangled***. 3. To minimise the number of CNOTs small ansatz had to be deployed, so ***EfficientSU2*** and ***RealAmlitudes*** ansatz algorithms were used to create very compact circuits, with ***linear entaglement*** and ***no repetiting segments***.4. As the author is a ***naive physicist / chemist*** the solution relied on ***experimentation*** and ***analytics*** rather than deep knowledge of the relevant theory.  Part 2: Final Challenge - VQE for LiH molecule In this part, you will simulate LiH molecule using the STO-3G basis with the PySCF driver.    Goal Experiment with all the parameters and then find your best ansatz. You can be as creative as you want!<jupyter_code>from qiskit_nature.drivers import PySCFDriver

molecule = 'Li 0.0 0.0 0.0; H 0.0 0.0 1.5474'
driver = PySCFDriver(atom=molecule)
qmolecule1 = driver.run()<jupyter_output><empty_output><jupyter_text>Hamiltonian and qubit reduction  Find molecule details<jupyter_code>def molstat(qmol):
    print('No of electrons: ', qmol.num_alpha + qmol.num_beta)
    print('Molecular orbitals: ', qmol.num_molecular_orbitals)
    print('Spin orbitals: ', 2 * qmol.num_molecular_orbitals)
    print('No of qubits: ', 2 * qmol.num_molecular_orbitals)
    print('Repulsion energy: ', qmol.nuclear_repulsion_energy)
    
molstat(qmolecule1)<jupyter_output><empty_output><jupyter_text>Removal of molecular orbitalsThis step aims to eliminate unnecessary orbitals from the molecular model, which later resuts in a smaller number of qubits in a quantum simulation of the molecule. We achieve this by ***freezing core orbitals*** and ***removing (hopefully) vacant orbitals***. This was achieved experimentatlly rather than by reference to the molecular properties derived from the theoretical model.+ ***Freezing core orbitals alone reduced molecular orbitals by one (-1).***+ ***Removal of [3, 4] orbitals (initially) reduced further 2 molecular orbitals (-2).***The impact of removing non-core orbitals was tested by measuring the simulated electronic energy. Depending on the selection of removed orbitals, it was possible to obtain 8, 4, 3 and 2 qubit systems with different energy levels.- It was determined that removal of orbital [0] or [1] also removed all electrons from the model, and this had a determental effect to later analysis. - Removal of orbital [5] generated an uncaught orbital error. - The only orbitals useful for further investigation were 2, 3 and 4. - Removal of either [3] or [4] alone produced high energy levels. - Removal of orbital [2] resulted in a 4 qubit system, which approximated the required energy level, still too high.In further investigation we considered orbital pairs [2, 3], [3, 4], [2, 4] and a tripple [2, 3, 4). - Removal of all [2, 3, 4] resulted in a two qubit system with a single CNOT, its low complexity did not allow its effective training to reach the threshold. - Pairs [2, 3] and [2, 4} produced 3 qubit circuits of 2 CNOTs each, which again could not be optimised below the required energy threshold. - A pair [3, 4] eventually mapped into 4 qubits and the acceptable energy level.<jupyter_code>from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem
from qiskit_nature.transformers import FreezeCoreTransformer
tr = FreezeCoreTransformer(freeze_core=True, remove_orbitals=[3, 4])
qmolecule2 = tr.transform(qmolecule1) 
molstat(qmolecule2)<jupyter_output><empty_output><jupyter_text>Creation of a quantum "problem"<jupyter_code>problem = ElectronicStructureProblem(driver, [tr])
second_q_ops = problem.second_q_ops()
main_op = second_q_ops[0] # electronic operator
print(main_op)<jupyter_output><empty_output><jupyter_text>Test symmetries and qubit reductionBy exploring some well-known symmetries, ParityMaper converter is used to perform a reduction of 2 further qubits (when possible). With an 'auto' option, this is done automatically and if any Z2 symmetries are detected they are then provided in the form of Pauli strings (+1, -1). Alternatively, the specified Z2 symmetries can be explicitly removed by relying on the list of tapering values. Here after the 'auto' removal, the list of tapering values is null and no further qubit reduction is possible.+ ***The number of quibts required for this simulation has been decreased by 2 to the total number of 4.***<jupyter_code>from qiskit_nature.mappers.second_quantization import ParityMapper, BravyiKitaevMapper, JordanWignerMapper
from qiskit_nature.converters.second_quantization.qubit_converter import QubitConverter
from qiskit.quantum_info import Pauli

# converter = QubitConverter(mapper=JordanWignerMapper(), two_qubit_reduction=True, z2symmetry_reduction='auto')
# converter = QubitConverter(mapper=BravyiKitaevMapper(), two_qubit_reduction=True, z2symmetry_reduction='auto')
# converter = QubitConverter(mapper=ParityMapper(), two_qubit_reduction=True, z2symmetry_reduction=[1, -1])
converter = QubitConverter(mapper=ParityMapper(), two_qubit_reduction=True, z2symmetry_reduction='auto')

num_particles = (problem.molecule_data_transformed.num_alpha, problem.molecule_data_transformed.num_beta)
qubit_op = converter.convert(main_op, num_particles=num_particles)
print('No of particles: ', num_particles)

pauli_symm = converter.z2symmetries.find_Z2_symmetries(qubit_op)
print(pauli_symm, '\n')
print('No of qubits: ', qubit_op.num_qubits)<jupyter_output><empty_output><jupyter_text>Ansatz Preparation  Init statePreparation of the ansatz initial state as defined by Hartree Fock approximation (no qubit corelations considered)<jupyter_code>from qiskit_nature.circuit.library import HartreeFock
num_particles = (problem.molecule_data_transformed.num_alpha, problem.molecule_data_transformed.num_beta)
num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals
init_state = HartreeFock(num_spin_orbitals, num_particles, converter)
init_state.draw(output='mpl')<jupyter_output><empty_output><jupyter_text>AnsatzThere are several different algorithms for ansatz generation and the possibiity of creatimg custom ansatz. To reduce the number of CNOTs, it is important that the number of qubits is as small as possible, the ansatz circuit is small, and no repeated structures are included. In addition to the provided ansatz algorithms, two promising additional ones have been identified in the qiskit library.+ Qiskit ***EfficientSU2*** and ***RealAmlitudes*** ansatz algorithms were identified as capable of creating very compact circuits.<jupyter_code>from qiskit.circuit.library import TwoLocal, EfficientSU2
from qiskit_nature.circuit.library import UCCSD, PUCCD, SUCCD

# Arguments for ansatz generation
num_particles = (problem.molecule_data_transformed.num_alpha,
             problem.molecule_data_transformed.num_beta)
num_spin_orbitals = 2 * problem.molecule_data_transformed.num_molecular_orbitals
rotation_blocks = ['ry', 'rz']
entanglement_blocks = 'cx'
entanglement = 'linear' # 'sca' # 'circular' # 'full' # Linear produces min CNOTs
repetitions = 1 # No repetitions to keep ansatz small
skip_final_rotation_layer = True 

# ansatz = TwoLocal(qubit_op.num_qubits, rotation_blocks, entanglement_blocks, reps=repetitions, 
#                   entanglement=entanglement, skip_final_rotation_layer=skip_final_rotation_layer)             # Possible
# ansatz = PUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)                          # Reject
# ansatz = SUCCD(converter,num_particles,num_spin_orbitals,initial_state = init_state)                          # Reject
# ansatz = UCCSD(converter,num_particles,num_spin_orbitals,initial_state = init_state, reps=1)                  # Possible
# ansatz = RealAmplitudes(num_qubits=qubit_op.num_qubits, entanglement=entanglement, reps=repetitions)          # Possible
ansatz = EfficientSU2(num_qubits=qubit_op.num_qubits, entanglement=entanglement, reps=repetitions)              # Best

# Add the initial state
ansatz.compose(init_state, front=True, inplace=True)

display(ansatz.draw(output='mpl', scale=1))
print('Ansatz number of parameters: ', ansatz.num_parameters)
print('Ansatz ordered parameters: ', ansatz.ordered_parameters)<jupyter_output><empty_output><jupyter_text>Backend<jupyter_code>from qiskit import Aer
backend = Aer.get_backend('statevector_simulator')<jupyter_output><empty_output><jupyter_text>OptimizerSeveral optimisers have been tested, the best results were obtained with COBYLA.The electronic energy resulting from diagonalisation needs to be checked against the VQE result<jupyter_code>from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA, SLSQP
optimizer = COBYLA(maxiter=10000) # Up to 20,000 were used with some configurations
# optimizer = SPSA(maxiter=5000)
from qiskit_nature.algorithms.ground_state_solvers.minimum_eigensolver_factories import NumPyMinimumEigensolverFactory
from qiskit_nature.algorithms.ground_state_solvers import GroundStateEigensolver
import numpy as np 

def exact_diagonalizer(problem, converter):
    solver = NumPyMinimumEigensolverFactory()
    calc = GroundStateEigensolver(converter, solver)
    result = calc.solve(problem)
    return result

result_exact = exact_diagonalizer(problem, converter)
exact_energy = np.real(result_exact.eigenenergies[0])
print("Exact electronic energy", exact_energy)
result_exact<jupyter_output><empty_output><jupyter_text>VQEVQE finds the ansatz parameter values for the circuit to produce its lowest electronic energy.<jupyter_code>from qiskit.algorithms import VQE
from IPython.display import display, clear_output

# Print and save the data in lists
def callback(eval_count, parameters, mean, std):  
    # Overwrites the same line when printing
    display("Evaluation: {}, Energy: {}, Std: {}".format(eval_count, mean, std))
    clear_output(wait=True)
    counts.append(eval_count)
    values.append(mean)
    params.append(parameters)
    deviation.append(std)

counts = []
values = []
params = []
deviation = []

# Set initial parameters of the ansatz
# We choose a fixed small displacement 
# So all participants start from similar starting point

try:
    initial_point = [0.01] * len(ansatz.ordered_parameters)
except:
    initial_point = [0.01] * ansatz.num_parameters

algorithm = VQE(ansatz,
                optimizer=optimizer,
                quantum_instance=backend,
                callback=callback,
                initial_point=initial_point)

result = algorithm.compute_minimum_eigenvalue(qubit_op)
print(result)<jupyter_output><empty_output><jupyter_text>Scores<jupyter_code># Store results in a dictionary
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes import Unroller

# Unroller transpile your circuit into CNOTs and U gates
pass_ = Unroller(['u', 'cx'])
pm = PassManager(pass_)
ansatz_tp = pm.run(ansatz)
cnots = ansatz_tp.count_ops()['cx']
score = cnots

accuracy_threshold = 4.0 # in mHa
energy = result.optimal_value

# TwoLocal
result_dict = {
    'optimizer': optimizer.__class__.__name__,
    'mapping': converter.mapper.__class__.__name__,
    'ansatz': ansatz.__class__.__name__,
    'rotation blocks': rotation_blocks,
    'entanglement_blocks': entanglement_blocks,
    'entanglement': entanglement,
    'repetitions': repetitions,
    'skip_final_rotation_layer': skip_final_rotation_layer,
    'energy (Ha)': energy,
    'error (mHa)': (energy-exact_energy)*1000,
    'pass': (energy-exact_energy)*1000 <= accuracy_threshold,
    '# of parameters': len(result.optimal_point),
    'final parameters': result.optimal_point,
    '# of evaluations': result.optimizer_evals,
    'optimizer time': result.optimizer_time,
    '# of qubits': int(qubit_op.num_qubits),
    '# of CNOTs': cnots,
    'score': score}

# Plot the results
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1)
ax.set_xlabel('Iterations')
ax.set_ylabel('Energy')
ax.grid()
fig.text(0.7, 0.75, f'Energy: {result.optimal_value:.3f}\nScore: {score:.0f}')
plt.title(f"{result_dict['optimizer']}-{result_dict['mapping']}\n{result_dict['ansatz']}")
ax.plot(counts, values)
ax.axhline(exact_energy, linestyle='--')
fig_title = f"\
{result_dict['optimizer']}-\
{result_dict['mapping']}-\
{result_dict['ansatz']}-\
Energy({result_dict['energy (Ha)']:.3f})-\
Score({result_dict['score']:.0f})\
.png"
fig.savefig(fig_title, dpi=300)

# Display and save the data
import pandas as pd
import os.path
filename = 'results_h2.csv'
if os.path.isfile(filename):
    result_df = pd.read_csv(filename)
    result_df = result_df.append([result_dict])
else:
    result_df = pd.DataFrame.from_dict([result_dict])
result_df.to_csv(filename)
result_df[['optimizer','ansatz', '# of qubits', '# of parameters','rotation blocks', 'entanglement_blocks',
    'entanglement', 'repetitions', 'error (mHa)', 'pass', 'score']]
# Check your answer using following code
from qc_grader import grade_ex5
freeze_core = True # change to True if you freezed core electrons
grade_ex5(ansatz,qubit_op,result,freeze_core)
# Submit your answer. You can re-submit at any time.
from qc_grader import submit_ex5
submit_ex5(ansatz,qubit_op,result,freeze_core)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09041114885946058, 'fraction_numerical': 0.010885391388754074, 'mean_word_length': 5.868571428571428, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '564714', 'license': 'Apache-2.0', 'n_tokens_mistral': 4568, 'n_tokens_neox': 4205, 'n_words': 1348, 'path': 'solutions by participants/ex5/ex5-JacobCybulski-3cnot-2.440695mHa-16params.ipynb', 'repo_name': 'fazliberkordek/ibm-quantum-challenge-2021', 'size': '14431'}