Initial State Prep. using Active Spaces for Quantum Chem

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!

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

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

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.

Hi @Christophe_Pere I have same problem. I adapted my code to the H2O problem and still getting the error. What I did is int(np.log(len(initial_state))/np.log(2)) to get the number of wires.

from pyscf import gto, scf, ci
from pennylane.qchem import import_state
import pennylane as qml
from pennylane import numpy as np
elec     = 8*2  #contados normalmente
orbs     = 5*2 #contados normalmente entonces 10 qubits para O_2
a_elec   = 2
a_orbs   = 2
delta_sz = 0
q        = 0
m        = 1




R = 1.28
# create the 
xyz_format = [["O", (0.028, 0.054, 0.0)], ["H", (0.986, 1.610, 0.0)],["H", (1.855, 0.002, 0.0)],]

s          = [  j   for i in xyz_format for j in i[1]]

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



#####################################################################
#####################################################################
#####################################################################

symbols  = [i[0] for i in xyz_format]

print(symbols)

x        = np.array(s, requires_grad=True)

core, active = qml.qchem.active_space(electrons        = elec,
                                      orbitals         = orbs,
                                      active_electrons = a_elec, 
                                      active_orbitals  = a_orbs )

def H(x):
    molecule   = qml.qchem.Molecule(symbols, x, charge=q, mult=m)
    global qubits 
    hh, qubits = qml.qchem.molecular_hamiltonian(molecule,
                                                 method           = 'openfermion',
                                                 active_electrons = a_elec,
                                                 active_orbitals  = a_orbs)
    # print(qubits)
    return hh

H(x)
# initial_state = qml.qchem.initial_state_state(electrons=elec, orbitals=orbs)
initial_state = wf_cisd



singlet, doublet = qml.qchem.excitations(electrons = a_elec , 
                                         orbitals  = qubits, 
                                         delta_sz  = delta_sz)

print('INITIAL STATE', initial_state.shape,'to wires', int(np.log(len(initial_state))/np.log(2)) )



dev = qml.device("lightning.qubit", wires=orbs)
@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    # qml.BasisState(initial_state, wires=wires)

    wires = int(np.log(len(initial_state))/np.log(2))

    qml.StatePrep(initial_state,wires=wires)
    
    [qml.SingleExcitation(phi = params[i] , wires = s) for i, s in enumerate(singlet)];
    [qml.DoubleExcitation(phi = params[j] , wires = d) for j, d in enumerate(doublet)];

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    # print(hamiltonian)
    return circuit(params, obs=hamiltonian, wires=range(qubits))

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift     = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res       = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad   = [circuit(params, obs=obs, wires=range(qubits)) for obs in grad_h]
    return np.array(grad)
print()
opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x     = qml.GradientDescentOptimizer(stepsize=0.8)

theta     = np.array(np.zeros(len(doublet)), requires_grad=True)

# store the values of the cost function
energy      = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs = 0.529177210903

for n in range(100):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond_length.append(np.linalg.norm(x[0:3] - x[3:6]) * bohr_angs)

    if n % 4 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]:.5f} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

A short version of the error is:

Cell In[10], line 127
    125 theta.requires_grad = True
    126 x.requires_grad = False
--> 127 theta, _ = opt_theta.step(cost, theta, x)
    129 # Optimize the nuclear coordinates
    130 x.requires_grad = True
...
Cell In[10], line 88
     86 hamiltonian = H(x)
     87 # print(hamiltonian)
---> 88 return circuit(params, obs=hamiltonian, wires=range(qubits))
...
Cell In[10], line 78
     72 @qml.qnode(dev, interface="autograd")
     73 def circuit(params, obs, wires):
     74     # qml.BasisState(initial_state, wires=wires)
     76     wires = int(np.log(len(initial_state))/np.log(2))
---> 78     qml.StatePrep(initial_state,wires=wires)
     80     [qml.SingleExcitation(phi = params[i] , wires = s) for i, s in enumerate(singlet)];
     81     [qml.DoubleExcitation(phi = params[j] , wires = d) for j, d in enumerate(doublet)];

Thank you in advance for your help :wink:

Hi @Raul_Guerrero ,

Can you please post your full error traceback? It’s much easier to find the source of the error when you have the full error information.

Hi @CatalinaAlbornoz, sorry for the delay in my answer. Because I was playing a lot with the code I sent to you again and the full error:

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

elec     = 8*2  #contados normalmente
orbs     = 7*2 #contados normalmente entonces 10 qubits para O_2
a_elec   = 4
a_orbs   = 4
delta_sz = 0
q        = 0
m        = 1




R = 1.28
# create the 
xyz_format = [["O", (0.      ,  0.      ,  0.119262)], ["H", (0.      ,  0.763239, -0.477047)],["H", ( 0.      , -0.763239, -0.477047)],]

s          = [  j   for i in xyz_format for j in i[1]]

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



#####################################################################
#####################################################################
#####################################################################

symbols      = [i[0] for i in xyz_format]

print(symbols)

x            = np.array(s, requires_grad=True)

core, active = qml.qchem.active_space(electrons        = elec,
                                      orbitals         = orbs,
                                      active_electrons = a_elec, 
                                      active_orbitals  = a_orbs )

def H(x):
    molecule   = qml.qchem.Molecule(symbols, x, charge=q, mult=m, unit='angstrom')
    global qubits 
    hh, qubits = qml.qchem.molecular_hamiltonian(molecule,
                                                 method           = 'openfermion',
                                                 active_electrons = a_elec,
                                                 active_orbitals  = a_orbs)
    # print(qubits)
    return hh

H(x)
# initial_state = qml.qchem.initial_state_state(electrons=elec, orbitals=orbs)
initial_state = wf_cisd

singlet, doublet = qml.qchem.excitations(electrons = a_elec , 
                                         orbitals  = qubits, 
                                         delta_sz  = delta_sz)

print('INITIAL STATE', initial_state.shape,'to wires', int(np.log(len(initial_state))/np.log(2)) )
print('act.orbs',a_orbs)
print('qubits',qubits)
print('singlet',len(singlet))
print('doublets',len(doublet))
bond1 = np.linalg.norm(x[0:3] - x[3:6])
bond2 = np.linalg.norm(x[3:6] - x[6:9])
print('initial bonds  {}, {} Ang'.format(bond1, bond2))

num_wires = int(np.log(len(initial_state))/np.log(2))
print('num_wires', num_wires)

dev = qml.device("lightning.qubit", wires=qubits)
@qml.qnode(dev, interface="autograd")
def circuit(params, obs, wires):
    wires_list = list(range(qubits))
    # qml.BasisState(initial_state, wires=wires)
    qml.StatePrep(initial_state, wires = wires_list)
    
    # [qml.SingleExcitation(phi = params[i] , wires = s) for i, s in enumerate(singlet)];
    [qml.DoubleExcitation(phi = params[j] , wires = d) for j, d in enumerate(doublet)];

    return qml.expval(obs)

def cost(params, x):
    hamiltonian = H(x)
    # print(hamiltonian)
    return circuit(params, obs=hamiltonian, wires=qubits)

def finite_diff(f, x, delta=0.01):
    """Compute the central-difference finite difference of a function"""
    gradient = []

    for i in range(len(x)):
        shift     = np.zeros_like(x)
        shift[i] += 0.5 * delta
        res       = (f(x + shift) - f(x - shift)) * delta**-1
        gradient.append(res)

    return gradient


def grad_x(params, x):
    grad_h = finite_diff(H, x)
    grad   = [circuit(params, obs=obs, wires=qubits) for obs in grad_h]
    return np.array(grad)
###########################################################
########################################################### Execution
###########################################################
opt_theta = qml.GradientDescentOptimizer(stepsize=0.4)
opt_x     = qml.GradientDescentOptimizer(stepsize=0.4)

theta     = np.array(np.zeros(len(doublet)), requires_grad=True)

# store the values of the cost function
energy      = []

# store the values of the bond length
bond_length = []

# Factor to convert from Bohrs to Angstroms
bohr_angs   =1 #0.529177210903

for n in range(100):

    # Optimize the circuit parameters
    theta.requires_grad = True
    x.requires_grad = False
    theta, _ = opt_theta.step(cost, theta, x)

    # Optimize the nuclear coordinates
    x.requires_grad = True
    theta.requires_grad = False
    _, x = opt_x.step(cost, theta, x, grad_fn=grad_x)

    energy.append(cost(theta, x))
    bond1 = np.linalg.norm(x[0:3] - x[3:6])* bohr_angs
    bond2 = np.linalg.norm(x[3:6] - x[6:9])* bohr_angs
    bond_length.append([bond1,bond2])
    
    if bond1 > 4 or bond2 > 4:
        raise ValueError('Bonds are diverging {},{}'.format(bond1, bond2))

    if n % 1 == 0:
        print(f"Step = {n},  E = {energy[-1]:.8f} Ha,  bond length = {bond_length[-1]} A")

    # Check maximum component of the nuclear gradient
    if np.max(grad_x(theta, x)) <= 1e-05:
        break

The error:

converged SCF energy = -74.9644048239943
E(RCISD) = -75.01467131660718  E_corr = -0.05026649261287289
myci 111
CISD-based state vector: 
[0. 0. 0. ... 0. 0. 0.]
['O', 'H', 'H']
INITIAL STATE (16384,) to wires 14
act.orbs 4
qubits 8
singlet 8
doublets 18
initial bonds  0.9685650182625841, 1.526478 Ang
num_wires 14
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[16], line 134
    132 theta.requires_grad = True
    133 x.requires_grad = False
--> 134 theta, _ = opt_theta.step(cost, theta, x)
    136 # Optimize the nuclear coordinates
    137 x.requires_grad = True

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/optimize/gradient_descent.py:93, in GradientDescentOptimizer.step(self, objective_fn, grad_fn, *args, **kwargs)
     75 def step(self, objective_fn, *args, grad_fn=None, **kwargs):
     76     """Update trainable arguments with one step of the optimizer.
     77 
     78     Args:
   (...)
     90         If single arg is provided, list [array] is replaced by array.
     91     """
---> 93     g, _ = self.compute_grad(objective_fn, args, kwargs, grad_fn=grad_fn)
     94     new_args = self.apply_grad(g, args)
     96     # unwrap from list if one argument, cleaner return

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/optimize/gradient_descent.py:122, in GradientDescentOptimizer.compute_grad(objective_fn, args, kwargs, grad_fn)
    104 r"""Compute gradient of the objective function at the given point and return it along with
    105 the objective function forward pass (if available).
    106 
   (...)
    119     will not be evaluted and instead ``None`` will be returned.
    120 """
    121 g = get_gradient(objective_fn) if grad_fn is None else grad_fn
--> 122 grad = g(*args, **kwargs)
    123 forward = getattr(g, "forward", None)
    125 num_trainable_args = sum(getattr(arg, "requires_grad", False) for arg in args)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/_grad.py:166, in grad.__call__(self, *args, **kwargs)
    163     self._forward = self._fun(*args, **kwargs)
    164     return ()
--> 166 grad_value, ans = grad_fn(*args, **kwargs)  # pylint: disable=not-callable
    167 self._forward = ans
    169 return grad_value

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/autograd/wrap_util.py:20, in unary_to_nary.<locals>.nary_operator.<locals>.nary_f(*args, **kwargs)
     18 else:
     19     x = tuple(args[i] for i in argnum)
---> 20 return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/_grad.py:184, in grad._grad_with_forward(fun, x)
    178 @staticmethod
    179 @unary_to_nary
    180 def _grad_with_forward(fun, x):
    181     """This function is a replica of ``autograd.grad``, with the only
    182     difference being that it returns both the gradient *and* the forward pass
    183     value."""
--> 184     vjp, ans = _make_vjp(fun, x)  # pylint: disable=redefined-outer-name
    186     if vspace(ans).size != 1:
    187         raise TypeError(
    188             "Grad only applies to real scalar-output functions. "
    189             "Try jacobian, elementwise_grad or holomorphic_grad."
    190         )

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/autograd/core.py:10, in make_vjp(fun, x)
      8 def make_vjp(fun, x):
      9     start_node = VJPNode.new_root()
---> 10     end_value, end_node =  trace(start_node, fun, x)
     11     if end_node is None:
     12         def vjp(g): return vspace(x).zeros()

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/autograd/tracer.py:10, in trace(start_node, fun, x)
      8 with trace_stack.new_trace() as t:
      9     start_box = new_box(x, t, start_node)
---> 10     end_box = fun(start_box)
     11     if isbox(end_box) and end_box._trace == start_box._trace:
     12         return end_box._value, end_box._node

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/autograd/wrap_util.py:15, in unary_to_nary.<locals>.nary_operator.<locals>.nary_f.<locals>.unary_f(x)
     13 else:
     14     subargs = subvals(args, zip(argnum, x))
---> 15 return fun(*subargs, **kwargs)

Cell In[16], line 93
     91 hamiltonian = H(x)
     92 # print(hamiltonian)
---> 93 return circuit(params, obs=hamiltonian, wires=qubits)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/workflow/qnode.py:1185, in QNode.__call__(self, *args, **kwargs)
   1183 if qml.capture.enabled():
   1184     return qml.capture.qnode_call(self, *args, **kwargs)
-> 1185 return self._impl_call(*args, **kwargs)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/workflow/qnode.py:1165, in QNode._impl_call(self, *args, **kwargs)
   1162     override_shots = kwargs["shots"]
   1164 # construct the tape
-> 1165 self.construct(args, kwargs)
   1167 original_grad_fn = [self.gradient_fn, self.gradient_kwargs, self.device]
   1168 self._update_gradient_fn(shots=override_shots, tape=self._tape)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/logging/decorators.py:61, in log_string_debug_func.<locals>.wrapper_entry(*args, **kwargs)
     54     s_caller = "::L".join(
     55         [str(i) for i in inspect.getouterframes(inspect.currentframe(), 2)[1][1:3]]
     56     )
     57     lgr.debug(
     58         f"Calling {f_string} from {s_caller}",
     59         **_debug_log_kwargs,
     60     )
---> 61 return func(*args, **kwargs)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/workflow/qnode.py:980, in QNode.construct(self, args, kwargs)
    978 with pldb_device_manager(self.device):
    979     with qml.queuing.AnnotatedQueue() as q:
--> 980         self._qfunc_output = self.func(*args, **kwargs)
    982 self._tape = QuantumScript.from_queue(q, shots)
    984 params = self.tape.get_parameters(trainable_only=False)

Cell In[16], line 83
     81 wires_list = list(range(qubits))
     82 # qml.BasisState(initial_state, wires=wires)
---> 83 qml.StatePrep(initial_state, wires = wires_list)
     85 # [qml.SingleExcitation(phi = params[i] , wires = s) for i, s in enumerate(singlet)];
     86 [qml.DoubleExcitation(phi = params[j] , wires = d) for j, d in enumerate(doublet)];

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/capture/capture_meta.py:89, in CaptureMeta.__call__(cls, *args, **kwargs)
     85 if enabled():
     86     # when tracing is enabled, we want to
     87     # use bind to construct the class if we want class construction to add it to the jaxpr
     88     return cls._primitive_bind_call(*args, **kwargs)
---> 89 return type.__call__(cls, *args, **kwargs)

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

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

Thank you in advance,
Raúl

Hi @Raul_Guerrero ,

As per the error message the error is generated in this line

qml.StatePrep(initial_state, wires = wires_list)

The error tells you that the State vector must have shape (2**wires,) or (batch_size, 2**wires).

The number of wires is the number of qubits.

There are 8 qubits so 2**wires=2**8=256. This means that initial_state cannot have shape 16384.

Additional note:
You don’t need to turn them into a list, you could simply write:

qml.StatePrep(initial_state, wires=range(qubits))

I’m also not sure why you need wires as an argument for circuit.

I hope this helps!

Hi @CatalinaAlbornoz , I understand, I’ll try again to fix it and let you know.

Yes, the wires thing in the circuit argument is because I start to train myself with the geometry optimization tutorial It’s just a “historical” reason.

Thanks, as always your help is very valuable.

Raúl

1 Like

Hi @CatalinaAlbornoz.
At the end I found my real error. Since I was creating a “full electron” state it was imposible to introduce this state in the Q-simulator because I was also trying to consider the active space which much more smaller and the number of wires did not coincide.

That’s why I enter my other question in this q&a.

Thank you

Thanks for adding this message here @Raul_Guerrero !

Let’s continue the conversation with Stepan on the other thread in that case :smiley:

1 Like