# Unable to construct Hamiltonian

Thank you for your kind reply @maliasadi and @CatalinaAlbornoz. I will summarize the list of calculations I have done below over the past few days.

I ran the calculation systematically with 321G,631G and 631G* and cc-PVDZ. 6-31G* and ccPVDZ gave FileNotFoundError error. Initially my guess was that , in the linux shell , it is unable to write the file PSPCz_Psi4_631GD_psi4_6-31G* because in linux shell it can read it only as PSPCz_Psi4_631GD_psi4_6-31G * .out. However, when I tried to run the same molecule usign ccPVDZ, it gave the same error. To verify the observation, I also did the same Hamiltonian construction for the molecule NH3 using basis sets, STO3G,321G,6-31G,6-31G,ccPVDZ. And of them ran successfully. At this moment, I am unable to reason out why for the molecule PSPCz (45 Atoms, 200 electrons with Active Space=2 and no. of active electrons=2) is giving errors while using higher basis sets such as 6-31G or ccPVDZ. Kindly suggest how to proceed about this.

Please let me know if I am not clear.
Also, wish you a merry Christmas and happy a New Year 2022.

Thank you

Hi @raghavv, thanks for the question and sorry for the late reply.

It seems to me that, for larger basis sets, the pyscf output hdf5 file is not created at all either because the Hartree-Fock iteerations crash or you ran out of disc space. Could you please run independednt pyscf calculations for the same molecule follwoing the examples provided here and see if you can successfully compute the Hartree-Fock energy with pyscf? If you can run pyscf successfully, then we can track the problem through the pipeline. Please let me know if you face any issues with the calculations.

Dear @sjahangiri, thank you for the kind response. Sorry for the delayed response. I have been trying to construct the Hamiltonian from outside of Pennylane as you have suggested. I have not had much of a success yet. I am working on it and will update once I figure it out. If you there are any other suggestions/directions, kindly request you to share. Will keep things updated sir. Thank you very much.

Thanks @raghavv for the update. If you post the molecular geometry here, in xyz format for example, I can also give it a try.

Dear Sir @sjahangiri, Below I have given the coordinates of water which I am trying to use PySCF or any other chemistry package to construct Hamiltonian and have it written in the form of .HDF5 file format.

H -0.0399 -0.0038 0.0
O 1.5780 0.8540 0.0
H 2.7909 -0.5159 0.0

Thanks for your kind support sir.

Hi @raghavv. The following code creates the molecular Hamiltonian for water step by step. Could you please run the code for both water and PSPCz and let me know at what stage the code breaks? Please modify the molecular info given between the dashed lines for the new molecule. Thanks!

import pennylane as qml
import openfermion
from pennylane import numpy as np
from openfermionpyscf import run_pyscf

#------------------------
symbols = ["O", "H", "H"]
x = np.array([0.0000000000,  0.0000000000,  0.1501668459,
0.0000000000,  0.7580826553, -0.4856334229,
0.0000000000, -0.7580826553, -0.4856334229])
mult = 1
charge = 0
basis = 'sto-3g'
core = [0, 1, 2, 3]
active = [4, 5]
filename = './water'
#------------------------

geometry = [[s, tuple(x[3 * i : 3 * i + 3])] for i, s in enumerate(symbols)]

molecule = openfermion.MolecularData(geometry, basis, mult, charge, filename=filename)

run_pyscf(molecule, run_scf=1, verbose=0)

molecule_ = openfermion.MolecularData(filename=filename)

ht = molecule_.get_molecular_hamiltonian(occupied_indices=core, active_indices=active)

hf = openfermion.transforms.get_fermion_operator(ht) # fermion hamiltonian

hq = openfermion.transforms.jordan_wigner(hf) # qubit hamiltonian

hp = qml.qchem.convert_observable(hq) # PennyLane hamiltonian

print(hp)


The output for water should be

  (-74.19370628334303) [I0]
+ (-0.15460393042639708) [Z2]
+ (-0.15460393042639708) [Z3]
+ (0.13038589407848444) [Z0]
+ (0.13038589407848444) [Z1]
+ (0.138284189609364) [Z0 Z2]
+ (0.138284189609364) [Z1 Z3]
+ (0.14767808197773874) [Z0 Z3]
+ (0.14767808197773874) [Z1 Z2]
+ (0.14966950988795605) [Z2 Z3]
+ (0.22003977334376112) [Z0 Z1]
+ (-0.00939389236837476) [Y0 Y1 X2 X3]
+ (-0.00939389236837476) [X0 X1 Y2 Y3]
+ (0.00939389236837476) [Y0 X1 X2 Y3]
+ (0.00939389236837476) [X0 Y1 Y2 X3]

Dear Sir, @sjahangiri . Thanks for sharing the code. It ran smooth for both water and PSPCz. Earlier, I had tried within Pennylane using “qchem.molecular_hamiltonian” (as shown above) and it took me a very long time to complete the job and create a corresponding HDF5 with Hamiltonian written to it. It had ran on only one core. The code you suggested ran like a gem and completed within few minutes. Thank you so much for sharing the code. I have attached the code with PSPCz molecular coordinates for your kind reference. The qubit Hamiltonian that was output is pasted below for your quick reference.

PSPCz_Hamiltonian_Externally.py (3.2 KB)

(-322.1043666152533) [I0]

• (8.87674183162617) [Z2]
• (8.87674183162617) [Z3]
• (10.218927690523312) [Z0]
• (10.218927690523312) [Z1]
• (-0.006447940283844799) [Y1 Y3]
• (-0.006447940283844799) [X1 X3]
• (0.004959753825945462) [Y0 Y2]
• (0.004959753825945462) [X0 X2]
• (0.02545094913478638) [Z0 Z2]
• (0.02545094913478638) [Z1 Z3]
• (0.02674734205767252) [Z0 Z3]
• (0.02674734205767252) [Z1 Z2]
• (0.07156347202814556) [Z0 Z1]
• (0.08691649983952954) [Z2 Z3]
• (0.14484917836670383) [Y0 Z1 Y2]
• (0.14484917836670383) [X0 Z1 X2]
• (0.14484917836670383) [Y1 Z2 Y3]
• (0.14484917836670383) [X1 Z2 X3]
• (-0.0064479402838447985) [Y0 Z1 Y2 Z3]
• (-0.0064479402838447985) [X0 Z1 X2 Z3]
• (-0.0012963929228861414) [Y0 Y1 X2 X3]
• (-0.0012963929228861414) [X0 X1 Y2 Y3]
• (0.0012963929228861414) [Y0 X1 X2 Y3]
• (0.0012963929228861414) [X0 Y1 Y2 X3]
• (0.004959753825945462) [Z0 Y1 Z2 Y3]
• (0.004959753825945462) [Z0 X1 Z2 X3]

I am not sure why using “qchem.molecular_hamiltonian” took such a long time to create the Hamiltonian or it even crashed. Here, I am assuming calling the function ‘qchem.molecular_hamiltonian’ would still go to PySCF and ask for a Hartree-Fock level of calculation to be run. Then why does it take a longer time. Can you kindly clarify.

Please let me know if I am not clear.

Thank you.

Hi @raghavv. This is a very interesting molecule!

Yes, qml.qchem.molecular_hamiltonian just calls openfermion under the hood.

I computed the PSPCz Hamiltonian with both
qml.qchem.molecular_hamiltonian and the stepwise code using sto-3g basis and 2 active electrons and 2 active orbitals. Both methods take the same time, ~11 minutes, and the size of the hdf5 files are both 4.8G. I am not sure why qml.qchem.molecular_hamiltonian takes longer time in your examples with larger basis sets.

However, there is no fundamental difference between the two approaches. We are actually considering to refactor qml.qchem.molecular_hamiltonian and ask the users to call the external libraries on their side. So, you can safely use this code example if it works better for your system.

Please feel free to let us know if you have any other questions.

1 Like

Dear Sir @sjahangiri. Sorry for the delayed response. I was trying to understand the code and the physics behind it all. As of now working on it will update here. Thank you very much for kind responses. Will keep you posted.

Dear Sir @sjahangiri, hope you are well.

I tried for quite sometime to combine the Pennylane’s mol-geom-opt (https://pennylane.ai/qml/demos/tutorial_mol_geo_opt.html) and adaptive circuit (https://pennylane.ai/qml/demos/tutorial_adaptive_circuits.html) tutorials to try to optimized molecule using Full Configuration Interaction method. After several attempts, I managed to combine the two codes , however I am unable to run optimization on the molecule. It only runs a single point calculation. I am not sure what I am missing here. It would be really helpful if you could kindly look into the code and suggest what I am missing.

Kindly let me know if I am not clear enough.

Thank youCombined_Code_That_Run-For-Question.py (7.8 KB)

Dear @raghavv,

Thank you very much for sharing your question.

I noted that you are using/mixing functionalities from PennyLane, OpenFermion and PySCF. This is typically not a good strategy, it makes the code less readable and harder to debug. All the functionalities you need for molecular geometry optimization are available in PennyLane.

I suggest the following:

1- For now, please keep the adaptive method to select the excitation operations independent from the geometry optization code. I understand from your message that this is working fine.

2- For the geometry optimization I recommend following the tutorial Optimization of molecular geometries.

• I noted you have used multiple functions to build the Hamiltonian of LiH molecule. This should be simpler. For the LiH molecule in a minimal basis sto-3g, you have 4 electrons and a total of 6 orbitals ( => 12 qubits). That would give you a relatively large Hamiltonian with 631 Pauli strings. The latter can be reduced by defining an active space. For example, you could chose 2 electrons and 3 active orbitals (=> 6 qubits, 62 Pauli strings). This can be done as follows:
import pennylane as qml

symbols = ["Li", "H"]
x = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 2.969280527], requires_grad=True)

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


Please, note that from this point on your algorithm (including the quantum circuit) to optimize the geometry of LiH, within this approximation, is identical to the one in the above mentioned tutorial.

Could you please try this and let us know whether it works for you or not? If it works, you can proceed to refine the simulation by including more orbitals and excitation operations.

Please, do not hesitate to come back to us. It is always a pleasure to help!

Best regards.

2 Likes

Dear sir @Alain_Delgado_Gran. Thank you for the kind response. I tried the method you suggested to calculate the Hamiltonian of the molecule.
I carried out the calculation for LiH and it worked. I was unable to construct Hamiltonian for the bigger molecule such as PSPCz which I had mentioned in the first comment. That is why @sjahangiri had suggested me a work around. However, as you mentioned when the Hamiltonian is borrowed from Openfermion, I guess the nuclear gradient did not run.

I am trying to optimize molecule PSPCz using Pennylane. I will keep you updated on that.

Meanwhile, I have a couple questions regarding the XYZ format that the code prints out and the molecular orbitals. I kindly request you to help me understand these.

1. The following code has the number ‘3’ in it which properly prints out the coordinate in XYZ format. I am not sure what is the significance of the no. ‘3’ in it. Also, it goes like 3i , 3i+1 and 3*i+2. I am not sure what these help us achieve, can you please help me understand.

Also looks like the coordinates being printed are in Bohr units. I would like to convert them into angstrom units. I tried something like ({x[3 * i]:.4f} *0.529177210903), but it did not work.

for i, atom in enumerate(symbols):
print(f" {atom} {x[3 * i]:.4f} {x[3 * i + 1]:.4f} {x[3 * i + 2]:.4f}")

1. Is there a way to extract molecular orbital information and their energies in Pennylane.

Thank you.

Dear @raghavv,

I am glad to hear that the geometry optimization of the LiH molecule went well.

1. The following code has the number ‘3’ in it which properly prints out the coordinate in XYZ format. I am not sure what is the significance of the no. ‘3’ in it. Also, it goes like 3 i , 3 i+1 and 3*i+2. I am not sure what these help us achieve, can you please help me understand.

This 3 comes from the fact that the array x is a one-dimensional array storing 3N coordinates [x_0, x_1, ... , x_{3N-1}] where N is the number of atoms in your molecule. Thus, if you want to print the nuclear coordinates in the XYZ format, you have to account for this. For example, you can access the coordinates of the first atom (i=0) as: x = x[0], y = x[1], z=x[2] .

For printing the coordinates in Angstroms you could simply do:

x_a = x*0.529177210903

print("%s %4s %8s %8s" % ("symbol", "x", "y", "z"))
for i, atom in enumerate(symbols):
print(f"  {atom}    {x_a[3 * i]:.4f}   {x_a[3 * i + 1]:.4f}   {x_a[3 * i + 2]:.4f}")

1. Is there a way to extract molecular orbital information and their energies in Pennylane.

Yes. You can use first the meanfield function to generate the file containing the molecular orbitals data for LiH. Next, you can use OpenFermion MolecularData class to access the energies and the wave functions of the orbitals. Please, see the example below:

import pennylane as qml
from pennylane import numpy as np
from openfermion import MolecularData

symb = ["Li", "H"]
x = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 2.969280527], requires_grad=True)

lih = qml.qchem.meanfield(symb, x)
mol = MolecularData(filename=lih)

print("Number of orbitals: ", mol.n_orbitals)
print("Orbital energies: ", mol.orbital_energies)
print("Expansion coefficients: ", mol.canonical_orbitals)


FYI: New utility functions are being implemented in PennyLane to access this information natively, and also for visualizing atomic and molecular orbitals and much more! They will be released soon. Stay tuned!

Hope this helps. Do not hesitate to get back to us if you need further clarifications.

Thank you.

1 Like

Dear Sir @Alain_Delgado_Gran. Thank you for sharing the information regarding significance of ‘3’ in coordinates and helping to convert into from Bohr to Angstrom as well. Also, thanks for sharing the information on how to get the molecular orbitals.

Sometime ago, I had also tried to get molecular orbital energies using OpenFermion. I have written some lines to get the frontier molecular orbital energies and their positions. I am just sharing the lines below.

“”"
import pennylane as qml
from pennylane import numpy as np
from openfermion import MolecularData

symb = [“Li”, “H”]
x = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 2.969280527], requires_grad=True)

lih = qml.qchem.meanfield(symb, x)
mol_LiH = MolecularData(filename=lih)

print("Number of orbitals: ", mol_LiH.n_orbitals)
print("Orbital energies: ", mol_LiH.orbital_energies)
print("Expansion coefficients: ", mol_LiH.canonical_orbitals)

molecule = mol_LiH
basis = ‘STO3G’

total_MOs = (molecule.canonical_orbitals.shape[0])
print(“No. of orbitals generated for basis set”,basis, molecule.canonical_orbitals.shape[0])

num_electrons = molecule.n_electrons
occupied_orbitals = round(num_electrons/2)
print(“No. of Occupied MOs for basis set”,basis, occupied_orbitals)
unoccupied_orbitals = total_MOs - occupied_orbitals
print(“No. of Unoccupied MOs for basis set”,basis, unoccupied_orbitals)

HOMO_position = round((num_electrons/2)-1)
LUMO_position = round((num_electrons/2)+0)

print(“The HOMO and LUMO positions are”,HOMO_position,LUMO_position)

HOMO_energy = molecule.orbital_energies[HOMO_position]
LUMO_energy = molecule.orbital_energies[LUMO_position]
H_L_Gap = LUMO_energy - HOMO_energy

print(“The HOMO and LUMO energies are”,round(HOMO_energy,5),round(LUMO_energy,5),“respectively in Hartrees”)
print(“The HOMO-LUMO gap is”, round(H_L_Gap,5),“in Hartrees”)
print(“The HOMO-LUMO gap is”, round(H_L_Gap*27.2114,3),“in eV”)

“”"

Also, regarding constructing the Hamiltonian, can you kindly help me understand the reason for ’ [0] ’ in the following line.

“def H(x):
return qml.qchem.molecular_hamiltonian(symbols, x, active_electrons=2, active_orbitals=3))[0]”

I am also working on obtaining ground state energies of a few small molecules. I will share the information with you after consolidating the results.

Thank you very much for your kind help.

Hi @raghavv,

Thank you for the update.

Also, regarding constructing the Hamiltonian, can you kindly help me understand the reason for ’ [0] ’ in the following line.

Please, note that the molecular_hamiltonian function returns a tuple containing the built Hamiltonian and the number of qubits needed to represent it. Since you only need the Hamiltonian (for the example we have been discussing) the H(x) function returns the first element of the tuple.

Thank you!

1 Like

Dear Sir @Alain_Delgado_Gran, thank you for your kind response, I understand the importance of [0] now.

I am also learning the the chemical reaction tutorial of Pennylane (https://pennylane.ai/qml/demos/tutorial_chemical_reactions.html) and I am just curious to know if we can calculate the zero point vibrational energy in Pennylane while calculating dissociation energy of the molecule.

Thank you.

Hi @raghavv,

In principle you can calculate the zero point vibrational energy in PennyLane while calculating the dissociation energy of the molecule.

@Alain_Delgado_Gran is looking into what functionalities are available for this.

1 Like

Hi @raghavv,

Yes, you can compute the zero-point energy correction from the vibrational frequencies of the molecule. In the harmonic approximation, they are given by the non-zero eigenvalues of the energy Hessian matrix defined as \frac{\partial^2E({\bf R})}{\partial R_i \partial R_j}. Here,E({\bf R}) is the total electronic energy of the ground or excited states of the molecule and {\bf R} denotes the optimized nuclear coordinates for a given electronic state (e.g. the ground state).

At present, we don’t have a tutorial showing how to compute the Hessian matrix elements \frac{\partial^2E({\bf R})}{\partial R_i \partial R_j} but I try to outline here how this can be done by combining different functionalities available in PennyLane. To that aim, I encourage you to check this paper by Mitarai et al. about computing analytical energy derivatives for the VQE algorithm. This contribution reports the expression to compute the second-order derivative of the energy with respect to the nuclear coordinates (See Eqs. (8) and (10) ). Evaluating this expression requires the following:

• The equilibrium geometry of the molecule and the optimized circuit parameters to prepare its (electronic) ground-state.
• Compute expectation values of the first- and second-order derivatives
of the electronic Hamiltonian (defined for the optimized geometry) with respect to the nuclear coordinates. The Hamiltonian derivatives can be computed using finite-difference approximations or more accurate automatic differentiation tools available in PennyLane.
• Compute the second-order derivative of the circuit (expectation value
of the Hamiltonian) with respect to the gate parameters (Hessian of the circuit). See example here

Please, see the attached notebook with a working example to compute the second-order energy derivative for the trihydrogen cation.

Hope this helps!

Thank you (3.7 KB)

Hi @raghavv, here’s the file @Alain_Delgado_Gran had shared but in .py format. You should be able to paste it into a Jupyter notebook and run it.

second_derivative.py (2.2 KB)

1 Like

Thank you sir @Alain_Delgado_Gran for such a beautiful explanation with appropriate references. I will go thorugh them.

Thank you @CatalinaAlbornoz for sharing the code.

I tried to run the code, upon running the line " hessian = qml.jacobian(qml.grad(energy, argnum=0))(opt_param) ", it gave me the following error.

## “”"

AttributeError Traceback (most recent call last)
/tmp/ipykernel_1076/1091415110.py in
1 # compute the Hessian of the quantum circuit
2 energy = expval(H(coord), opt_param)
----> 3 hessian = qml.jacobian(qml.grad(energy, argnum=0))(opt_param)

179
180 if len(argnum) == 1:
–> 181 return _jacobian(func, argnum[0])(*args, **kwargs)
182
183 return np.stack([_jacobian(func, arg)(*args, **kwargs) for arg in argnum]).T

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 return nary_f
22 return nary_operator

55 (out1, out2, …) then the Jacobian has shape (out1, out2, …, in1, in2, …).
56 “”"
—> 57 vjp, ans = _make_vjp(fun, x)
58 ans_vspace = vspace(ans)
59 jacobian_shape = ans_vspace.shape + vspace(x).shape

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()

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)
16 if isinstance(argnum, int):
17 x = args[argnum]

99 “”“Evaluates the gradient function, and saves the function value
100 calculated during the forward pass in :attr:.forward.”""
102 self._forward = ans

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 return nary_f
22 return nary_operator

116 difference being that it returns both the gradient and the forward pass
117 value."""
–> 118 vjp, ans = _make_vjp(fun, x)
119
120 if not vspace(ans).size == 1:

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()

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)
16 if isinstance(argnum, int):
17 x = args[argnum]

/tmp/ipykernel_1076/3145879875.py in funct(params)
2 def expval(obs, params):
3 def funct(params):
----> 4 return circuit(params, obs, wires=range(qubits))
5 return funct

~/psi4conda/envs/pennylane/lib/python3.9/site-packages/pennylane/qnode.py in call(self, *args, **kwargs)
604 if self.mutable or self.qtape is None:
605 # construct the tape
–> 606 self.construct(args, kwargs)
607
608 # Execute the tape.

~/psi4conda/envs/pennylane/lib/python3.9/site-packages/pennylane/qnode.py in construct(self, args, kwargs)
523
524 with self.qtape:
–> 525 self.qfunc_output = self.func(*args, **kwargs)
526
527 if not isinstance(self.qfunc_output, Sequence):

/tmp/ipykernel_1076/3280061623.py in circuit(params, obs, wires)
6 qml.DoubleExcitation(params[0], wires=[0, 1, 2, 3])
7 qml.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
----> 8 return qml.expval(obs)

~/psi4conda/envs/pennylane/lib/python3.9/site-packages/pennylane/measure.py in expval(op)
233 if not isinstance(op, (Observable, qml.Hamiltonian)):
234 raise qml.QuantumFunctionError(
–> 235 “{} is not an observable or Hamiltonian: cannot be used with expval”.format(op.name)
236 )
237

AttributeError: ‘tensor’ object has no attribute ‘name’

“”".
I am not sure why this error pops up. Can you help me troubleshoot this.

Thankyou.