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