Initial State Prep. using Active Spaces for Quantum Chem

PennyLane Help
1

Hello! Following the Initial state preparation for quantum chemistry tutorial, is there anyway to create the CISD initial state, but for an active-space, like of which we can specify in Pennylane’s qml.qchem.molecular_hamiltonian ? I tried to do so, but it’s giving me an error

Here is the self-contained code:

import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np
import time

from pyscf import gto, scf, ci
from pennylane.qchem import import_state
from pennylane import numpy as np

def circuit(theta, initial_state):
    # prepares reference state
    qml.StatePrep(initial_state, wires=range(qubits))
        
    # apply 2 gates for simplicty:
    qml.SingleExcitation(params[0], wires=[1, 3])
    qml.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
                             
    # returns expectation value of the ansatz prepared from this quantum circuit:   
    return qml.expval(hamiltonian)

def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0, active_electrons = active_electrons, 
                                           active_orbitals = active_orbitals)[0]

symbols = ["O", "H", "H"]
active_electrons = 2; active_orbitals = 3
qubits = active_orbitals * 2
theta = np.array([0.0] * 2, requires_grad=True)
x = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0], requires_grad=True)
hamiltonian = H(x)

# create the H2O molecule
mol = gto.M(atom=[["O", (0.028, 0.054, 0.0)],
                  ["H", (0.986, 1.610, 0.0)],
                  ["H", (1.855, 0.002, 0.0)]], charge=0)

# Create cisd initial state for the selected active-space
# prepare for active space of 2-electrons, 3-orbitals
mol.nelectron = 2
mol.build()

# perfrom restricted Hartree-Fock and then CISD
myhf = scf.RHF(mol).run()
myci = ci.CISD(myhf).run()
wf_cisd = import_state(myci, tol=1e-1)
print(f"CISD-based state vector: \n{np.round(wf_cisd.real, 4)}")

# Finally, calculate observable
print(circuit(theta, wf_cisd))

And here is the error in question:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 49
     46 print(f"CISD-based state vector: \n{np.round(wf_cisd.real, 4)}")
     48 # Finally, calculate observable
---> 49 print(circuit(theta, wf_cisd))

Cell In[1], line 12, in circuit(theta, initial_state)
     10 def circuit(theta, initial_state):
     11     # prepares reference state
---> 12     qml.StatePrep(initial_state, wires=range(qubits))
     14     # apply 2 gates for simplicty:
     15     qml.SingleExcitation(params[0], wires=[1, 3])

File ~/anaconda3/lib/python3.11/site-packages/pennylane/ops/qubit/state_preparation.py:176, in StatePrep.__init__(self, state, wires, id)
    174     state = math.reshape(state, (1, state.shape[0]))
    175 if state.shape[1] != 2 ** len(self.wires):
--> 176     raise ValueError("State vector must have shape (2**wires,) or (batch_size, 2**wires).")
    178 param = math.cast(state, np.complex128)
    179 if not math.is_abstract(param):

ValueError: State vector must have shape (2**wires,) or (batch_size, 2**wires).

Thanks so much!

2

Hi @ImranNasrullah,

The problem with your code is that you don’t have the following:

  • Correct number of qubits
  • qml.device()
  • In the circuit function, you pass theta and use params

Checking the wf_cisd parameter, I found a shape of 16384 values corresponding to 14 qubits.

Here is a code working:

import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np
import time

from pyscf import gto, scf, ci
from pennylane.qchem import import_state
from pennylane import numpy as np

num_qubits = 14
dev = qml.device("default.qubit", wires=num_qubits, )

@qml.qnode(dev, interface='auto')#, interface="autograd")
def circuit(theta, initial_state):
    # prepares reference state
    qml.StatePrep(initial_state, wires=range(num_qubits))
        
    # apply 2 gates for simplicty:
    qml.SingleExcitation(theta[0], wires=[1, 3])
    qml.DoubleExcitation(theta[1], wires=[0, 1, 4, 5])
                             
    # returns expectation value of the ansatz prepared from this quantum circuit:   
    return qml.expval(hamiltonian)

def H(x):
    return qml.qchem.molecular_hamiltonian(symbols, x, charge=0, active_electrons = active_electrons, 
                                           active_orbitals = active_orbitals)[0]

symbols = ["O", "H", "H"]
active_electrons = 2; active_orbitals = 3
qubits = active_orbitals * 2
theta = np.array([0.0] * 2, requires_grad=True)
x = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0], requires_grad=True)
hamiltonian = H(x)

# create the H2O molecule
mol = gto.M(atom=[["O", (0.028, 0.054, 0.0)],
                  ["H", (0.986, 1.610, 0.0)],
                  ["H", (1.855, 0.002, 0.0)]], charge=0)

# Create cisd initial state for the selected active-space
# prepare for active space of 2-electrons, 3-orbitals
mol.nelectron = 2
mol.build()

# perfrom restricted Hartree-Fock and then CISD
myhf = scf.RHF(mol).run()
myci = ci.CISD(myhf).run()
wf_cisd = import_state(myci, tol=1e-1)
print(f"CISD-based state vector: \n{np.round(wf_cisd.real, 4)}")

# Finally, calculate observable
print(circuit(theta, wf_cisd))

Output:

converged SCF energy = -54.6879622228754
E(RCISD) = -54.68810612450983  E_corr = -0.0001439016343809385
CISD-based state vector: 
[0. 0. 0. ... 0. 0. 0.]
-74.87230051148896
3

Hello @Christophe_Pere , thanks for the response!

I was looking to conserve the number of qubits to be 6, because VQE already gets very expensive when moving to multiple molecular interactions. Would you know of any way to modify the wf_cisd construction so that it is represented in a smaller amount of qubits?

Thanks

4

Hi @ImranNasrullah,

I understand what you mean. The number of coefficients in your wavelength is 43 when you run len(myci.ci). np.log2(43) ~ 5.43, so 6 if we make a round, like the 6 qubits you want. But the function import_state returns the state of 2^M where M is the number of spin orbitals. In H_20 you have 14 spin orbitals. So the size of the wf_cisd is 2^{14} so 16384.