# Gradient in differentiable Hartree-Fock demo

Hello, I am following the Differentiable Hartree-Fock | PennyLane Demos. Nevertheless, when applying another atomic species (but keeping the same amount of atoms, i.e., 2 atoms as in the demo) and active space, my pc is having struggles with the autograd.grad. The amount of orbitals and electrons increase compared with the demo. There is any way to reduce the memory usage of the autograd like storing the data and read it on each iteration step?

Best and thank you,
Raúl

Do you have some code that we could use to try to replicate your results? Something that we can copy-paste and run ourselves.

I touched many times the code, so I don’t have the same problem than before. The problem of memory is solved. Now, instead, I still having problems with the autograd.grad. As you can see in the code bellow there are 2 cases, H2 and C2. When I run the code for H2 the result is equal to the demo (until the cell before the last one in the demo). Here the parameters (there is only one) for the `DoubleExcitation` (one single wire) is optimizing well and also the atomic positions. When calculate the C2 the problem comes with the last two cells (I separated the codes with big ‘##’ blocks to represent the cells). Here the parameters (now they are 3) and the `DoubleExcitation` (3 wires) are not optimizing, but the atomic positions are evolving aparently well. Why the parameters are not changing or I am might missing something?

``````from pyscf import gto, scf, ci, mcscf
import pennylane as qml
from pennylane import numpy as np
import pennylane as qml

# For C2

nelecas  = 2*2
ncas     = 2*2
delta_sz = 0
q        = 0
m        = 1

symbols   = ['C', 'C']
geometry    = np.array([[0.        , 0.        , 1.        ],
[1.22975607, 0.71      , 1.        ]],

#For H2
# nelecas  = 1*2
# ncas     = 1*2
# delta_sz = 0
# q        = 0
# m        = 1
# symbols = ["H", "H"]
# # optimized geometry at the Hartree-Fock level
# geometry = np.array([[-0.672943567415407, 0.0, 0.0],
#                      [ 0.672943567415407, 0.0, 0.0]], requires_grad=True)

mol       = qml.qchem.Molecule(symbols,
geometry,
charge = q,
mult   = m)

H, qubits = qml.qchem.molecular_hamiltonian(mol,
method           = 'openfermion',
active_electrons = nelecas,
active_orbitals  = ncas)
print("Number of qubits   : ", qubits)
print("The Hamiltonian is : ", H)
singles, doubles = qml.qchem.excitations(electrons=2, orbitals=qubits)

doubles          = doubles[:3]
print('NOTE: doubles truncated to {}'.format(len(doubles))
print('singlet :', singles)
print('doubles :', doubles)
num_theta        = len(doubles)
#####################################################################
#####################################################################
#####################################################################
#####################################################################
dev = qml.device("default.qubit", wires=qubits)

def energy(mol):
def circuit(*args):
state = qml.qchem.hf_state(electrons=nelecas,orbitals=qubits)
qml.BasisState(state, wires=range(qubits))
print('FULL ARGS in circuit function', )
[print('      ',i) for i in args]
print('\nChecking the actual parameters for the qml.DoubleExcitation ')
[print('       parameters: ', args[0][0][i], 'for wire {} ='.format(i),dd) for i, dd in enumerate(doubles)]
[qml.DoubleExcitation(phi=args[0][0][i], wires=dd) for i, dd in enumerate(doubles)]

# for i, dd in enumerate(doubles):
#qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])

H = qml.qchem.molecular_hamiltonian(mol,
args=args[1:],
active_electrons = nelecas,
active_orbitals  = ncas)[0]
return qml.expval(H)
return circuit

#####################################################################
#####################################################################
#####################################################################
#####################################################################
# print(circuit_param)
for n in range(36):
print('---\niteration {}'.format(n))
args = [circuit_param, geometry]
mol = qml.qchem.Molecule(symbols, geometry)
g_param       = grad(energy(mol), argnum = 0)(*args)
print(g_param)
circuit_param = circuit_param - 0.25 * g_param[0]

mol = qml.qchem.Molecule(symbols, geometry)
forces   = -grad(energy(mol), argnum = 1)(*args)
geometry = geometry + 0.5 * forces

if n % 2 == 0:
print(f'n: {n}, E: {energy(mol)(*args):.8f}, Force-max: {abs(forces).max():.8f}',g_param[0])
``````

Best
Raúl

Update of my previous post and solution
Hi again @CatalinaAlbornoz,

I was playing again with the code and I realize that, for the case of C_2, the parameters can be updated on each iteration if the active space is made of 2 electrons and 2 orbitals (so 4 qubits). As a test I try to add an extra `DoubleExcitation` wire (wires = [0,1,3,2]) by hand to the H_2, incase it was a problem of the wires. Here the calculation it was performed as well the original (wires = [0,1,2,3]).

What could be the problem with the active space?

########################################### Solution
Finally I get the solution. I was doing an amateur mistake, in `qchem.active_space` function I was introducing by mistake only 2 electrons instead of the active_space number of electrons (nelecas) value.

Best and thank you,
Raúl

1 Like

It’s great to see that you solved it @Raul_Guerrero !

Thanks for posting the solution here

1 Like

Hi @CatalinaAlbornoz, I am having a problem with the Autograd.grad. It is giving me nan values when is calculating the forces for certain atom. Is there a way to manipulate the grad function to check where the problem is?

Best,
Raúl

Hi @Raul_Guerrero and thanks for the question. Could you please post a minimal example that gives the error, so I can look into that? Thanks!

Hi @sjahangiri, thank you in advance for helping me. Here you can see the code to replicate the error.

Thank you again.
Best,
Raúl

``````from pyscf import gto, scf, ci, mcscf
from pennylane.qchem import import_state
import pennylane as qml
from pennylane import numpy as np
from pennylane.qchem.convert import _wfdict_to_statevector
import pennylane as qml
import warnings
warnings.filterwarnings('ignore')

all_e    = 10
all_o    = 7
nelecas  = 4
ncas     = 4
delta_sz = 0
q        = 0
m        = 1
units = 'bohr'
symbols  = ["H", "O", "H"]
geometry = np.array([[-0.0399, -0.0038, 0.0],
[1.5780, 0.8540, 0.0],
[2.7909, -0.5159, 0.0]],

mol       = qml.qchem.Molecule(symbols,
geometry,
charge = q,
mult   = m,
unit   = units)
# print(mol.coordinates,mol.charge)
H, qubits = qml.qchem.molecular_hamiltonian(mol,
method           = 'dhf',
active_electrons = nelecas,
active_orbitals  = ncas)

print("Number of qubits   : ", qubits)
print("The Hamiltonian is : ", H)
singles, doubles = qml.qchem.excitations(electrons=nelecas,
orbitals=qubits)
print('singlet {}:'.format(len(singles)), singles)
print('doubles {}:'.format(len(doubles)), doubles)
num_theta        = len(doubles)  + len(singles)

core, active = qml.qchem.active_space(electrons = all_e,
orbitals  = all_o,
mult      = m,
active_electrons = nelecas,
active_orbitals  = ncas)
print('index of core orbitals : ',core,'\nindex of active orbitals: ',active)

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

def energy(mol):
def circuit(*args):
state = qml.qchem.hf_state(electrons=nelecas,orbitals=qubits)
qml.BasisState(state, wires=range(qubits))
# [qml.SingleExcitation(phi=args[0][0][i], wires=dd) for i, dd in enumerate(singles)]
[qml.DoubleExcitation(phi=args[0][0][i], wires=dd) for i, dd in enumerate(doubles)]
H = qml.qchem.molecular_hamiltonian(mol,
args=args[1:], # args represents the initial values of the differentiable parameter
method           = 'dhf',
active_electrons = nelecas,
active_orbitals  = ncas,
wires=list(range(qubits)))[0]
expval = qml.expval(H)
return expval
return circuit

full_bs  = mol.basis_set
alpha    = [np.array(aa.alpha,requires_grad = True) for aa in full_bs]
coeff    = [np.array(cc.coeff,requires_grad = True) for cc in full_bs]

for n in range(36):
print('---\niteration {}'.format(n))
args = [circuit_param, geometry, alpha, coeff]
mol  = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff, unit = units, mult=m, charge=q)
####################################################################### PARAMS
g_param         = grad(energy(mol), argnum = 0)(*args)
circuit_param   = circuit_param  - 0.25 * g_param[0]
print('     params updated', circuit_param)
####################################################################### FORCES
geometry      =  geometry  + 0.5 * forces
print('     forces updated', forces)

it_energy     = energy(mol)(*args)

print(f'n: {n}, E: {it_energy:.8f}, Force-max: {abs(forces).max()}')

``````

The output

``````Number of qubits   :  8
singlet 8: [[0, 4], [0, 6], [1, 5], [1, 7], [2, 4], [2, 6], [3, 5], [3, 7]]
doubles 18: [[0, 1, 4, 5], [0, 1, 4, 7], [0, 1, 5, 6], [0, 1, 6, 7], [0, 2, 4, 6], [0, 3, 4, 5], [0, 3, 4, 7], [0, 3, 5, 6], [0, 3, 6, 7], [1, 2, 4, 5], [1, 2, 4, 7], [1, 2, 5, 6], [1, 2, 6, 7], [1, 3, 5, 7], [2, 3, 4, 5], [2, 3, 4, 7], [2, 3, 5, 6], [2, 3, 6, 7]]
index of core orbitals :  [0, 1, 2]
index of active orbitals:  [3, 4, 5, 6]
---
iteration 0
params updated [tensor([ 2.11707715e-02, -3.08582261e-05,  3.08592969e-05,
2.99495697e-02,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  1.92008549e-02,
-7.53191631e-07,  7.53191631e-07,  1.87801609e-02,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
forces updated [[ 1.38383484e+00  8.05306831e-01 -2.41822383e-17]
[-3.67198359e-01 -2.04905546e+00  7.33617710e-18]
[-1.01663648e+00  1.24374863e+00  1.68460612e-17]]
n: 0, E: -37.45151654, Force-max: 2.049055457778269
---
iteration 1
params updated [tensor([ 1.64942544e-02,  2.43051714e-03,  1.27473936e-04,
-7.69204734e-02,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  2.19253580e-02,
3.22623371e-05, -3.16163905e-06, -2.14856854e-02,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
forces updated [[ 8.42606699e+00 -7.15002667e+00 -5.12000000e+02]
[            nan             nan             nan]
[-1.06626881e+01 -3.89410800e+00  5.17693786e-17]]
n: 1, E: -45.84641294, Force-max: nan
---
iteration 2
``````

The error

``````---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[1], line 100
98 mol  = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff, unit = units, mult=m, charge=q)
99 ####################################################################### PARAMS
--> 100 g_param         = grad(energy(mol), argnum = 0)(*args)
101 circuit_param   = circuit_param  - 0.25 * g_param[0]

18 else:
19     x = tuple(args[i] for i in argnum)
---> 20 return unary_operator(unary_f, x, *nary_op_args, **nary_op_kwargs)

21 @unary_to_nary
23     """
24     Returns a function which computes the gradient of `fun` with respect to
25     positional argument number `argnum`. The returned function takes the same
26     arguments as `fun`, but returns the gradient instead. The function `fun`
27     should be scalar-valued. The gradient has the same type as the argument."""
---> 28     vjp, ans = _make_vjp(fun, x)
29     if not vspace(ans).size == 1:
30         raise TypeError("Grad only applies to real scalar-output functions. "

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

13 else:
14     subargs = subvals(args, zip(argnum, x))
---> 15 return fun(*subargs, **kwargs)

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)

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[1], line 78
76 [qml.SingleExcitation(phi=args[0][0][i], wires=dd) for i, dd in enumerate(singles)]
77 [qml.DoubleExcitation(phi=args[0][0][i], wires=dd) for i, dd in enumerate(doubles)]
---> 78 H = qml.qchem.molecular_hamiltonian(mol,
79                                     args=args[1:], # args represents the initial values of the differentiable parameter
80                                     method           = 'dhf',
81                                     active_electrons = nelecas,
82                                     active_orbitals  = ncas,
83                                     wires=list(range(qubits)))[0]
84 expval = qml.expval(H)
85 return expval

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:376, in molecular_hamiltonian(*args, **kwargs)
254 """molecular_hamiltonian(molecule, method="dhf", active_electrons=None, active_orbitals=None,\
255 mapping="jordan_wigner", outpath=".", wires=None, args=None, convert_tol=1e12)
256 Generate the qubit Hamiltonian of a molecule.
(...)
372
373 """
375 if len(args) != 0:
--> 376     return _molecular_hamiltonian_dispatch(*args, **kwargs)
378 method = kwargs.pop("symbols", None) or kwargs.pop("molecule", None)
379 if method is not None:

File ~/miniconda3/envs/full-PL/lib/python3.12/functools.py:909, in singledispatch.<locals>.wrapper(*args, **kw)
905 if not args:
906     raise TypeError(f'{funcname} requires at least '
907                     '1 positional argument')
--> 909 return dispatch(args[0].__class__)(*args, **kw)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:409, in _(molecule, method, active_electrons, active_orbitals, mapping, outpath, wires, args, convert_tol)
397 @_molecular_hamiltonian_dispatch.register(Molecule)
398 def _(
399     molecule,
(...)
407     convert_tol=1e12,
408 ):
--> 409     return _molecular_hamiltonian(
410         molecule.symbols,
411         molecule.coordinates,
412         molecule.name,
413         molecule.charge,
414         molecule.mult,
415         molecule.basis_name,
416         method,
417         active_electrons,
418         active_orbitals,
419         mapping,
420         outpath,
421         wires,
422         molecule.alpha,
423         molecule.coeff,
424         args,
426         convert_tol,
427     )

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:554, in _molecular_hamiltonian(symbols, coordinates, name, charge, mult, basis, method, active_electrons, active_orbitals, mapping, outpath, wires, alpha, coeff, args, load_data, convert_tol)
548 core, active = qml.qchem.active_space(
549     mol.n_electrons, mol.n_orbitals, mult, active_electrons, active_orbitals
550 )
552 requires_grad = args is not None
553 h = (
--> 554     qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)(*args)
556     else qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)()
557 )
559 if active_new_opmath():
560     h_as_ps = qml.pauli.pauli_sentence(h)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:246, in diff_hamiltonian.<locals>._molecular_hamiltonian(*args)
236 def _molecular_hamiltonian(*args):
237     r"""Compute the qubit hamiltonian.
238
239     Args:
(...)
243         Hamiltonian: the qubit Hamiltonian
244     """
--> 246     h_ferm = fermionic_hamiltonian(mol, cutoff, core, active)(*args)
248     return qubit_observable(h_ferm, mapping=mapping)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:187, in fermionic_hamiltonian.<locals>._fermionic_hamiltonian(*args)
177 def _fermionic_hamiltonian(*args):
178     r"""Compute the fermionic hamiltonian.
179
180     Args:
(...)
184         FermiSentence: fermionic Hamiltonian
185     """
--> 187     core_constant, one, two = electron_integrals(mol, core, active)(*args)
189     return fermionic_observable(core_constant, one, two, cutoff)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hamiltonian.py:120, in electron_integrals.<locals>._electron_integrals(*args)
111 def _electron_integrals(*args):
112     r"""Compute the one- and two-electron integrals in the molecular orbital basis.
113
114     Args:
(...)
118         tuple[array[float]]: 1D tuple containing core constant, one- and two-electron integrals
119     """
--> 120     _, coeffs, _, h_core, repulsion_tensor = scf(mol)(*args)
121     one = qml.math.einsum("qr,rs,st->qt", coeffs.T, h_core, coeffs)
122     two = qml.math.swapaxes(
123         qml.math.einsum(
124             "ab,cd,bdeg,ef,gh->acfh", coeffs.T, coeffs.T, repulsion_tensor, coeffs, coeffs
(...)
127         3,
128     )

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/pennylane/qchem/hartree_fock.py:143, in scf.<locals>._scf(*args)
140 rng = qml.math.random.default_rng(2030)
141 s = s + qml.math.diag(rng.random(len(s)) * 1.0e-12)
--> 143 w, v = qml.math.linalg.eigh(s)
144 x = v @ qml.math.diag(1.0 / qml.math.sqrt(w)) @ v.T
146 eigvals, w_fock = qml.math.linalg.eigh(
147     x.T @ h_core @ x
148 )  # initial guess for the scf problem

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/autoray/autoray.py:81, in do(fn, like, *args, **kwargs)
79 backend = _choose_backend(fn, args, kwargs, like=like)
80 func = get_lib_fn(backend, fn)
---> 81 return func(*args, **kwargs)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/numpy/linalg/linalg.py:1487, in eigh(a, UPLO)
1484     gufunc = _umath_linalg.eigh_up
1486 signature = 'D->dD' if isComplexType(t) else 'd->dd'
-> 1487 w, vt = gufunc(a, signature=signature, extobj=extobj)
1488 w = w.astype(_realType(result_t), copy=False)
1489 vt = vt.astype(result_t, copy=False)

File ~/miniconda3/envs/full-PL/lib/python3.12/site-packages/numpy/linalg/linalg.py:118, in _raise_linalgerror_eigenvalues_nonconvergence(err, flag)
117 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
--> 118     raise LinAlgError("Eigenvalues did not converge")

LinAlgError: Eigenvalues did not converge
``````

Thanks @Raul_Guerrero, I created a simple working example that you can use to debug the original code. Please also have a look at this discussion which is a bout a similar topic. Please let me know if you have any further questions. Thanks!

``````from autograd import grad
import pennylane as qml
from pennylane import numpy as np

symbols = ["O", "H", "H"]
geometry = np.array([[0.028,  0.054,  0.10],
[0.986,  1.610, -0.10],

mol = qml.qchem.Molecule(symbols, geometry)

dev = qml.device("default.qubit")

def energy(mol):
def circuit(*args):

# note that active_electrons=2, active_orbitals=2 in this example
qml.BasisState(np.array([1, 1, 0, 0]), wires=range(4))

# note that this simple circuit has only one gate
qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])

# note that active_electrons=2, active_orbitals=2 in this example
H = qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates,
active_electrons=2, active_orbitals=2, args=args[1:])[0]

return qml.expval(H)
return circuit

# number of zeros should match the number of circuit gates

for n in range(36):

args = [circuit_param, geometry]
mol = qml.qchem.Molecule(symbols, geometry)

g_param = grad(energy(mol), argnum = 0)(*args)
circuit_param = circuit_param - 0.25 * g_param[0]

forces = -grad(energy(mol), argnum = 1)(*args)
geometry = geometry + 0.5 * forces

print(f'n: {n}, E: {energy(mol)(*args):.8f}, Force-max: {abs(forces).max():.8f}')
``````
1 Like

Hi @sjahangiri, thank you for your answer. I read the discussion, as you suggest. I was implementing the code for water (the one that you shared me) and it worked perfectly. Nevertheless, I was extending its usage to other molecules and getting similar errors than the one in my previous posts. For instance, for a C_2 molecule starting with initial separation of 2.68 a.u. (8 electrons and 8 orbitals) and I am getting nan values for the forces, so the ‘code’ is giving me the `LinAlgError: Eigenvalues did not converge`. Is there any error in my inputs of the molecule? It will be better to move to the Accelerating VQEs with quantum natural gradient | PennyLane Demos tutorial?

Raúl

Hi @Raul_Guerrero, the C2 molecule has a more complicated electronic structure and it is likely that for that Hartree-Fock iteration has not converged properly. You might change the number of active electrons and active orbitals and see if the calculations converge properly. You might also saturate the molecule with extra hydrogen atoms to reduce the complexity of the electronic structure.

1 Like

Hi @sjahangiri thank you for the advice, I will try some less complex molecule.

A new question comes to my mind, Why `method="openfermion"` is not compatible with the `autograd.grad`? I noticed that `grad` only works with `dhf` method, is because is not numerical?

Best,
Raúl

The reason why openfermion doesn’t work is because it doesn’t have automatic differentiation.

Automatic differentiation is part of what makes PennyLane so special and powerful. It’s compatible only with the `method="dhf"` , which uses our differentiable Hartree Fock workflow.

1 Like