LiH_version_check

from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.aqua.components.optimizers import COBYLA, SLSQP
from qiskit import BasicAer, Aer
import numpy as np
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry import FermionicOperator
from qiskit.aqua import QuantumInstance


def molecule_info(atom, freeze_list, remove_list):
    driver = PySCFDriver(atom = atom, 
                         unit = UnitsType.ANGSTROM,
                         charge = 0, 
                         spin = 0, 
                         basis = 'sto3g')
    molecule = driver.run()
    repulsion_energy = molecule.nuclear_repulsion_energy
    
    # for restricted HF calculations num_alpha = num_beta
    num_particles = molecule.num_alpha + molecule.num_beta
    num_spin_orbitals = molecule.num_orbitals * 2
    
    remove_list = [x % molecule.num_orbitals for x in remove_list]
    freeze_list = [x % molecule.num_orbitals for x in freeze_list]
    remove_list = [x - len(freeze_list) for x in remove_list]
    remove_list += [x + molecule.num_orbitals - len(freeze_list)  for x in remove_list]
    freeze_list += [x + molecule.num_orbitals for x in freeze_list]
    
    ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)
    energy_shift = 0
    ferOp_red = []
    if len(freeze_list) != 0:
        ferOp_red, energy_shift = ferOp.fermion_mode_freezing(freeze_list)
        num_spin_orbitals -= len(freeze_list)
        num_particles -= len(freeze_list)
    if len(remove_list) != 0:
        ferOp_red = ferOp_red.fermion_mode_elimination(remove_list)
        num_spin_orbitals -= len(remove_list)
    return ferOp, ferOp_red, energy_shift, num_spin_orbitals, num_particles, repulsion_energy, molecule

atom = 'Li 0.0 0.0 0.0; H 0.0 0.0 1.547'
freeze_list = [0]
remove_list = [-3, -2]

ferOp, ferOp_red, energy_shift, num_spin_orbitals, num_particles, repulsion_energy, molecule = molecule_info(atom, freeze_list, remove_list)
qubitOp = ferOp.mapping(map_type = 'jordan_wigner', threshold = 10**(-8))
qubitOp_red = ferOp_red.mapping(map_type = 'jordan_wigner', threshold = 10**(-8))

def uccsd_full(mapping, num_spin_orbitals, num_particles):
    z2_symmetries = None
    skip_commute_test = True
    sq_list = None
        
    initial_state = HartreeFock(
        num_spin_orbitals,
        num_particles,
        qubit_mapping = mapping, 
        sq_list = sq_list)

    var_form = UCCSD(
        num_orbitals = num_spin_orbitals,
        num_particles = num_particles,
        initial_state = initial_state,
        two_qubit_reduction = False,
        qubit_mapping = mapping,
        excitation_type = 'sd',
        z2_symmetries = z2_symmetries,
        skip_commute_test = skip_commute_test)
    
    return initial_state, var_form

initial_state, var_form = uccsd_full('jordan_wigner', qubitOp.num_qubits, molecule.num_alpha + molecule.num_beta)
initial_state_red, var_form_red = uccsd_full('jordan_wigner', qubitOp_red.num_qubits, num_particles)

backend = BasicAer.get_backend("statevector_simulator")
optimizer = SLSQP(maxiter = 200, disp = True)
initial_point = np.zeros(var_form.num_parameters)
initial_point_red = np.zeros(var_form_red.num_parameters)

algo = VQE(qubitOp, var_form, optimizer, initial_point)
res = algo.run(QuantumInstance(backend))
vqe = res.eigenvalue.real + repulsion_energy

algo = VQE(qubitOp_red, var_form_red, optimizer, initial_point_red)
res = algo.run(QuantumInstance(backend))
vqe_red = res.eigenvalue.real + repulsion_energy + energy_shift