I am using Pennylane Quantum Chemistry to run a VQE-type algorithm. However, I noticed that only a few shorter bond-length configurations are available. For example, for LiH, the maximum bond-length configuration available is 2.1 Å. Is there a way to construct a Hamiltonian for a longer bond length?
Thank you for your help! I went through the tutorial, but the difficult part is determining the coordinates. For example, how can I find the coordinates for LiH at a bond length of 3 Å?
For bond lengths that are not in the PennyLane datasets you can look them up in experimental molecule databases. I got this same question a couple of years ago and I saved the answer so I’m sharing it below. Please note that some things may have changed either in the databases or in PennyLane, so if any of the steps below don’t work or don’t make sense anymore please let me know and I can review the changes in more detail. I hope this can at least get you started, and hopefully you can help me confirm whether the steps below still work! I look forward to hearing from you on how this goes.
Note that this experimental geometry is not necessarily the optimal geometry that you will need to use for qchem.molecular_hamiltonian. This is because the optimal geometry depends on the specific basis set that you choose and the specific method that you choose.
The basis set is the set of orbitals that will be combined together. The specific combination of orbitals for a given basis set is what we will try to optimize to find the ground state energy. The qchem.molecular_hamiltonian function in PennyLane has a ‘basis’ keyword argument which has ‘sto-3g’ as the default value. You can keep this value for now.
The ‘method’ keyword argument specifies how we want to calculate the molecular hamiltonian. The default is ‘dhf’ which means Differentiable Hartree Fock. This is a PennyLane-native method that is differentiable. The other option is using the ‘pyscf’ method, which uses an OpenFermion backend but which is not differentiable. You can keep ‘dhf’ for now.
Now, how do you find the optimal geometry that you will need to use for qchem.molecular_hamiltonian? You can look for the calculated geometries in the Computational Chemistry Comparison and Benchmark DataBase. Write H2 as your chemical formula, look for the method you’re using in the rows, and then look for the column corresponding to the basis set you’re using. The CISD method will show results similar to ‘dhf’, so we will look at this row. On this row you can look for the ‘sto-3g’ basis set (1st column), which is the one we’re using. You will notice that the calculated geometry for the CISD method and the sto-3g basis set is 0.735 Angstrom.
The last step we need to do is converting this bondlength in Angstrom to atomic units, which are the units that qchem.molecular_hamiltonian takes. The conversion is: 1 Angstrom = 1.8897259885789 a.u. Since the H2 molecule is comprised of only two nuclei, we can decide that one of these nuclei will lie in the origin (coordinates [0.0, 0.0, 0.0]), and the other one will lie on coordinates [0.0, 0.0, 0.735*1.8897259885789].
Finally we define the geometry as a numpy array of the two points in space. geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.735*1.8897259885789]])
Now that we know how to find the optimal geometry, we can use qchem.molecular_hamiltonian to calculate the Hamiltonian and the eigenvalues for such Hamiltonian. The details of how the eigenvalues are classically calculated can be found in the SciPy documentation.