import numpy as np
import frayedends as fe
from pyscf import fci, gto
from pyscf.geomopt import geometric_solver, as_pyscf_methodGeometry optimization with FrayedEnds and geomeTRIC
code
This tutorial is about optimizing molecular geometries using FrayedEnds and geomeTRIC. We use FrayedEnds and the fci implementation in pyscf to calculate the ground state energy and energy gradient for a given geometry and use geomeTRIC (via the geomeTRIC interface in pyscf) to adjust the molecular coordinates until the optimal geometry is reached.
To run this you need pyscf, numpy and frayedends installed.
FrayedEnds uses MADNESS (Multiresolution ADaptive Numerical Environment for Scientific Simulation) to calculate and refine molecular orbitals. To set up the MADNESS environment we create a MadWorld3D object, “thresh” is the numerical precision MADNESS follows when approximating functions using Multi-Resolution-Analysis (MRA). This world will need to be passed to all processes which use MADNESS.
world = fe.MadWorld3D(thresh=1e-4)MADNESS runtime initialized with 7 threads in the pool and affinity OFF
The main part of this implementation is the function that calculates the ground state energy and gradient of a given molecular geometry. Our function needs the MADNESS world, a MolecularGeometry object (which is this quantum chemistry package’s molecule object) and the number of orbitals which should be considered in the calculation. The rest of the input parameters specify how many iterations should be done at maximum and when the energy is considered converged.
To perform calculations with our orbital basis, the orbitals need to be represented by MADNESS using MRA. These MRA representations are C++ objects, which are transformed into “SavedFct3D” objects when passed to python (or “SavedFct2D” when doing 2D calculations). When these SavedFct3D objects are passed back to the C++ backend, the MRA representation can be loaded again.
To determine the energy, we first determine a set of molecular orbitals using the pair-natural-orbitals (PNO) algorithm described in J.S. Kottmann, F.A. Bischoff, E.F. Valeev, J. Chem. Phys. 152, 2020. Those orbitals are passed to python as “SavedFct3D” objects. The orbital refinement code, which is called later, will need to know which of these orbitals are considered frozen core orbitals and which are active. For this, it checks the orb.type attribute. In this calculation we set all orbitals to active, since (as of March 2026) frozen core orbitals are not refined and we want to refine all of our orbitals.
After normalizing the PNOs we start with the main algorithm, the orbital refinement procedure: 1. Compute one- and two-body integrals using the current orbital basis. 2. Use a many-body-method (in this case FCI) to determine the ground state energy and the one- and two-body reduced density matrices (RDMs). 3. Using the Optimization3D class, refine the orbitals. The orbital refinement needs the nuclear potential as a “SavedFct3D”, the nuclear repulsion, the current set of orbitals, as well as the one- and two-body RDMs as input. The algorithm itself is described in the “Theory” page at https://github.com/FabianLangkabel/FrayedEnds or at arXiv:2410.19116. At the end of this step the orbital basis is overwritten with the new refined orbitals, which are already orthonormal. 4. Repeat the procedure starting from step 1, and only break out of the loop once the energy has converged.
Finally we use the refined orbitals and final rdms to calculate the energy gradient with respect to the nuclear coordinates \(R_i\). During this calculation we use that our many-body-wave function and orbital basis are fully variational, which allows us to use the Hellmann-Feynman theorem. The gradient is thus given by \[\frac{ \partial E}{\partial R_{i}}=\sum_{k,l}\eta^{k}_{l}\bra{\phi_k}\frac{\partial V}{\partial R_{i}}\ket{\phi_l}+\frac{\partial c}{\partial R_i},\] where \(\{\phi_i\}\) is the orbital basis, \(\eta^k_l\) the one-body rdm, \(c\) the nuclear repulsion and \(V\) the nuclear Coulomb potential.
def energy_and_gradient(
world: fe.MadWorld3D,
molgeom: fe.MolecularGeometry,
n_orbitals: int,
maxiter_whole_alg=10,
maxiter_orbref=1,
e_convergence=1e-4,
):
geom_str = molgeom.get_geometry_string() # extracting the geometry and units as a string
# we calculate initial guess orbitals for the orbital refinement using pair natural orbitals (PNO)
# for more details on this method see: J.S. Kottmann, F.A. Bischoff, E.F. Valeev, J. Chem. Phys. 152, 2020
madpno = fe.MadPNO(world, geom_str[0], units=geom_str[1], n_orbitals=n_orbitals)
orbitals = madpno.get_orbitals() # initial guess orbitals as MRA functions, orb.type determines whether the orbital is 'active' or 'frozen_occ' (in this case all active)
world.set_function_defaults() # pno code might change some defaults of the numerical environment (world), we reset them to original values here
nuclear_repulsion = molgeom.get_nuclear_repulsion() # nuclear repulsion energy
Vnuc = molgeom.get_vnuc(world) # nuclear potential as MRA function
integrals = fe.Integrals3D(world) # class to compute different types of integrals
orbitals = integrals.orthonormalize(orbitals=orbitals) # orthonormalize pnos
current_energy = 0.0
for iteration in range(maxiter_whole_alg):
# using the guess orbitals we compute the one- and two-body integrals
# (h=<i|-\Delta/2+Vnuc|j> and g=(ij|kl) (ordering of g depends on chosen convention))
G = integrals.compute_two_body_integrals(
orbitals, ordering="chem"
) # pyscf fci code assumes chem ordering of g tensor
T = integrals.compute_kinetic_integrals(orbitals)
V = integrals.compute_potential_integrals(orbitals, Vnuc)
h = T + V
g = G
# with these integrals we can use fci to determine GS energy, one- and two-body rdms
e, fcivec = fci.direct_spin1.kernel(
h, g.elems, n_orbitals, molgeom.n_electrons
) # Computes the energy and the FCI vector
rdm1, rdm2 = fci.direct_spin1.make_rdm12(
fcivec, n_orbitals, molgeom.n_electrons
) # Computes the 1- and 2- body reduced density matrices
rdm2 = np.swapaxes(rdm2, 1, 2) # swapping axes to match convention used in orbital refinement code
print("iteration {} energy {:+2.7f}".format(iteration, e + nuclear_repulsion))
if abs(current_energy - (e + nuclear_repulsion)) < e_convergence:
break # stops the algorithm if energy is converged
current_energy = e + nuclear_repulsion
# using the rdms we now refine the orbitals, for more details see Theory at https://github.com/FabianLangkabel/FrayedEnds
# or arXiv:2410.19116
orb_refiner = fe.Optimization3D(world, Vnuc, nuclear_repulsion)
orbitals = orb_refiner.get_orbitals(
orbitals=orbitals, # all orbitals are active in this example
rdm1=rdm1,
rdm2=rdm2,
maxiter=maxiter_orbref, # maximum number of iterations the refinement algorithm does
) # orbitals are now the refined orbitals
# with the refined orbitals obtained the algorithm loops back to recompute integrals -> fci -> orbital refinement until convergence
# as soon as our energy is converged we can compute the energy gradient w.r.t. nuclear coordinates
grad = molgeom.compute_dR_dE(world, rdm1, orbitals)
return current_energy, np.array(grad)All that is left now is to use the pyscf interface for geomeTRIC to wrap the energy_and_gradient function, define a molecular geometry that should be optimized (in this case just h2, so the calculation does not take as long, but larger molecules are also possible) and set the desired number of orbitals. Finally we call the geomeTRIC optimizer.
def f(pyscf_mol):
molgeom = fe.MolecularGeometry.from_pyscf_mol(pyscf_mol, units="bohr")
e, g = energy_and_gradient(world, molgeom, n_orbitals=2)
return e, g
pyscf_mol = gto.M(atom="H 0.0 0.0 -0.5\nH 0.0 0.0 0.5", unit="angstrom") # initial guess molecule
fake_method = as_pyscf_method(pyscf_mol, f) # wrapping energy and gradient function to be compatible with geomeTRIC
new_mol = geometric_solver.optimize(fake_method) # geometry optimization with geomeTRIC
print("new geometry (a)")
print(new_mol.tostring())Output
geometric-optimize called with the following command line:
/Users/timo/miniforge3/envs/fe_test/lib/python3.11/site-packages/ipykernel_launcher.py --f=/Users/timo/Library/Jupyter/runtime/kernel-v3d98c8c9c29a3f4798fe31ce0453d96c3397fb41c.json
())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// ****** ))))))
## .%%%%%%/ (%%%%%%, ******* ))))))
.%%%%%%/ *%%%%%%, **. /******* .))))))
*%%%%%%/ (%%%%%% ********/*..,*/********* *))))))
#%%/ (%%%%%%, *********************/ )))))))
...
> Max-Grad < 4.50e-04
> RMS-Disp < 1.20e-03
> Max-Disp < 1.80e-03
> === End Optimization Info ===
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...
Geometry optimization cycle 1
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.500000 0.000000 0.000000 0.000000
H 0.000000 0.000000 0.500000 0.000000 0.000000 0.000000
WARN: Mole.unit (angstrom) is changed to Bohr
iteration 0 energy -1.1283885
iteration 1 energy -1.1294028
iteration 2 energy -1.1295112
iteration 3 energy -1.1295254
cycle 1: E = -1.12951119445 dE = -1.12951 norm(grad) = 0.108844
Step 0 : Gradient = 7.696e-02/7.696e-02 (rms/max) Energy = -1.1295111944
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.00000e-01
Geometry optimization cycle 2
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.400000 0.000000 0.000000 0.100000
H 0.000000 0.000000 0.400000 0.000000 0.000000 -0.100000
iteration 0 energy -1.1503284
iteration 1 energy -1.1509446
iteration 2 energy -1.1509989
cycle 2: E = -1.15094463084 dE = -0.0214334 norm(grad) = 0.0375677
Step 1 : Displace = 1.000e-01/1.000e-01 (rms/max) Trust = 1.000e-01 (=) Grad = 2.656e-02/2.656e-02 (rms/max) E (change) = -1.1509446308 (-2.143e-02) Quality = 0.977
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.33352e-01
Geometry optimization cycle 3
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.347292 0.000000 0.000000 0.052708
H 0.000000 0.000000 0.347292 0.000000 0.000000 -0.052708
iteration 0 energy -1.1491849
iteration 1 energy -1.1495901
iteration 2 energy -1.1496271
cycle 3: E = -1.14959007859 dE = 0.00135455 norm(grad) = 0.0653876
Step 2 : Displace = 5.271e-02/5.271e-02 (rms/max) Trust = 1.414e-01 (+) Grad = 4.624e-02/4.624e-02 (rms/max) E (change) = -1.1495900786 (+1.355e-03) Quality = -0.512
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 3.65454e-01
Geometry optimization cycle 4
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.373646 0.000000 0.000000 -0.026354
H 0.000000 0.000000 0.373646 0.000000 0.000000 0.026354
iteration 0 energy -1.1515166
iteration 1 energy -1.1520571
iteration 2 energy -1.1521021
cycle 4: E = -1.15205713809 dE = -0.00246706 norm(grad) = 0.0051986
Step 3 : Displace = 2.635e-02/2.635e-02 (rms/max) Trust = 2.635e-02 (-) Grad = 3.676e-03/3.676e-03 (rms/max) E (change) = -1.1520571381 (-2.467e-03) Quality = 0.883
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 4.27298e-01
Geometry optimization cycle 5
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.375922 0.000000 0.000000 -0.002276
H 0.000000 0.000000 0.375922 0.000000 0.000000 0.002276
iteration 0 energy -1.1515446
iteration 1 energy -1.1520850
iteration 2 energy -1.1521307
cycle 5: E = -1.15208500278 dE = -2.78647e-05 norm(grad) = 0.000901434
Step 4 : Displace = 2.276e-03/2.276e-03 (rms/max) Trust = 3.727e-02 (+) Grad = 6.374e-04/6.374e-04 (rms/max) E (change) = -1.1520850028 (-2.786e-05) Quality = 1.762
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 3.53205e-01
Geometry optimization cycle 6
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.376400 0.000000 0.000000 -0.000477
H 0.000000 0.000000 0.376400 0.000000 0.000000 0.000477
iteration 0 energy -1.1515470
iteration 1 energy -1.1520876
iteration 2 energy -1.1521334
cycle 6: E = -1.15208756153 dE = -2.55875e-06 norm(grad) = 1.61261e-05
Step 5 : Displace = 4.775e-04/4.775e-04 (rms/max) Trust = 5.271e-02 (+) Grad = 1.140e-05/1.140e-05 (rms/max) E (change) = -1.1520875615 (-2.559e-06) Quality = 4.449
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 3.46886e-01
Geometry optimization cycle 7
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 0.000000 0.000000 -0.376409 0.000000 0.000000 -0.000009
H 0.000000 0.000000 0.376409 0.000000 0.000000 0.000009
iteration 0 energy -1.1515470
iteration 1 energy -1.1520876
iteration 2 energy -1.1521334
cycle 7: E = -1.15208759768 dE = -3.61495e-08 norm(grad) = 5.69956e-08
Step 6 : Displace = 8.698e-06/8.698e-06 (rms/max) Trust = 7.454e-02 (+) Grad = 4.030e-08/4.030e-08 (rms/max) E (change) = -1.1520875977 (-3.615e-08) Quality = 192.882
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 3.46886e-01
Converged! =D
#==========================================================================#
#| If this code has benefited your research, please support us by citing: |#
#| |#
#| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |#
#| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |#
#| http://dx.doi.org/10.1063/1.4952956 |#
#==========================================================================#
Time elapsed since start of run_optimizer: 100.826 seconds
new geometry (a)
H 0.00000000 0.00000000 -0.37640864
H 0.00000000 0.00000000 0.37640864