Alchemical geometry relaxation Open Access

Chemical compound space is conceptually defined as the infinite set of all possible chemical compounds.^1,2 Extensive knowledge of chemical compound space based on quantum mechanical calculations only is made extremely challenging by the unfathomably large number of molecules that have to be taken into consideration. In practical terms, only a small portion of it can be screened using standard computational methods. In order to improve our understanding of chemical compound space as well as to enhance molecular design efforts, faster methods to scan chemical space are needed.^3,4 This can be achieved in two possible ways: either using an exact method at a low level of theory [cheap quantum machine (QM) methods or force fields] or using an approximate method for estimating high-level theory calculations. Both quantum machine learning (QML) and alchemical perturbation density functional theory (APDFT) take advantage of data previously acquired and generate new predictions using either interpolation (QML) or extrapolation (APDFT) techniques. QML^5,6 can be used effectively to predict several quantum properties;^7–12 recent studies^13–15 succeeded in machine learning the optimized molecular geometries, but most QML methods do not predict directly geometrical structures.

APDFT¹⁶ is an alternative to QML for the rapid screening of chemical space. QML comes at the initial cost of securing a large training set, and this is not always available for the desired level of theory, the desired property of interest, or a relevant chemical subspace. In this regard, APDFT is more versatile than QML because it requires the explicit calculation of only one reference molecule. From the reference molecule calculation and through the use of alchemical perturbations, it is possible to subsequently obtain estimates of molecular properties for a combinatorially large number of target molecules.¹⁶ APDFT accuracy depends on the perturbation order¹⁷ and on the basis set choice.¹⁸ If a high-quality calculation [e.g., coupled cluster singles doubles (CCSD)] is used as a reference, the error in the third order energy prediction can be smaller than the error made using a cheaper QM method, such as Hartree–Fock (HF) and MP2.^16,18

Alchemical changes and derivatives were at first introduced as a way to rationalize energy differences between molecules.^19–21 Successive works showed that APDFT is also capable of predicting many other molecular properties, such as electron densities, dipole moments,¹⁶ ionization potentials,^22,23 HOMO energy eigenvalues,²⁴ band structures,²⁵ and deprotonation energies.^26,27 APDFT is particularly convenient if applied to highly symmetric systems where the explicit calculation of only a few derivatives (all the others can be obtained by symmetry operations) is required, e.g., in the BN doping of polycyclic aromatic hydrocarbons.^28–30 Using pseudopotentials, quantum alchemy can also be applied to crystal structures^25,31–35 and can even suggest catalyst improvements.^36–40 

One of the major remaining challenges in APDFT is to include geometry relaxation as the equilibrium molecular structure changes from the reference to the target; this would constitute a huge improvement over fixed atoms’ “vertical” predictions. To the best of our knowledge, only two studies^24,41 have tried to predict equilibrium energies, applying APDFT to “non-vertical” energy predictions (predicting energies for a perturbation that includes changes in both nuclear charges and geometry).

In this paper, we propose an alternative approach: we alchemically predict the geometry energy derivatives (gradients and Hessians) for the target molecule and, subsequently, use them to relax the geometry (see Fig. 1). Alchemical forces are defined as the second order mixed derivatives of the energy with respect to nuclear positions and nuclear charges $(∂2E∂Z∂R)$⁠. Alchemical forces appear as non-diagonal terms in the generalized Hessian⁴² (the matrix that contains all the second derivatives of the energy with respect to nuclear charges, atoms’ positions, and electron number). Other off-diagonal terms correspond to nuclear Fukui functions^43–48  $(∂2E∂Ne∂R)$ and alchemical Fukui functions^{22,23,49–52} $(∂2E∂Ne∂Z)$⁠.

These functions have already been studied, while, to the best of our knowledge, no study has yet carried out an in-depth analysis of alchemical forces. To study the predictive power and applications of alchemical predicted gradients and Hessians, we just consider the restricted Hartree–Fock case. We think that an extension to unrestricted Hartree–Fock (UHF) or density functional theory (DFT) is possible, and some of the conclusions derived for RHF may also hold for other levels of theory.

In Sec. II, we will introduce the reader to the concept of APDFT as well as to the Coupled Perturbed Hartree–Fock (CPHF) method used to calculate alchemical derivatives. We will derive an analytical formulation of RHF alchemical forces. In Sec. III, we show that it is possible to compute higher order derivatives through numerical differentiation and those derivatives can be used to predict geometrical gradients and Hessians for the target molecules; we also discuss how basis sets affect the accuracy of alchemical predictions of energy and geometry. In the results section, we give some examples of how to use the predicted gradients and Hessians to relax the geometry structures of the target molecules.

Truncating the series in Eq. (3) to a certain order will give an approximation of E^T. We call APDFTn prediction the approximation that includes all terms up to the n-th order.

As one can see from Eqs. (4) and (5), the alchemical derivatives are a functional of perturbed electronic densities. These expressions are part of APDFT (alchemical perturbation density functional theory), as already introduced in Ref. 16. They offer an orbital-free way to approximate alchemical energy differences as perturbations of the reference electronic density.

Higher order derivatives can be obtained via numerical differentiation or analytically using either Møller–Plesset perturbation theory⁴¹ or a coupled perturbed approach couple perturbed Hartree-Fock (CPHF) or coupled perturbed Kohn-Sham DFT (CPKS).⁵⁵ The contribution of nuclear–nuclear repulsion to the alchemical derivatives can be obtained analytically as the derivative of a classical charge distribution.

Coupled perturbed methods are a powerful tool to calculate derivative properties of a monodeterminant wave function (HF, DFT). They represent the standard way to calculate polarization^56–58 and can also be used for alchemical perturbations.^28,30,50,59

We used the standard index notation for molecular orbitals, i, j, k for occupied orbitals, a, b, c for virtual orbitals, and p, q, r for arbitrary orbitals. We will use Greek letters λ, μ, ν, and σ as indices of atomic orbitals and capital letters I, J, and K as atoms’ indices.

In APDFT, any molecular property can be expressed as a function of the alchemical coordinate λ.

Solving the CPHF equation for any atom I leads to the evaluation of a response matrix U^I and of the first order derivatives (ΔV^I, C^I, F^I, ϵ^I) of the operators with respect to the nuclear charge Z_I.

This formulation is particularly convenient for highly symmetrical systems, where many of the derivatives ∂A/∂Z_I are symmetrically equivalent. Only a few explicit CPHF calculations are needed, while the number of possible targets grows combinatorially with the size of the system.

We pointed out in a previous paper¹⁸ that neglecting the derivatives of the basis set coefficients with respect to nuclear charges, i.e., calculating alchemical derivatives using the basis set of the reference molecule only, the alchemical series converges to the energy of the target molecule with the basis set of the reference (E^T[R]). The difference between this energy and the true target energy was named (alchemical) basis set error: ΔE_BS ≔ E^T[R] − E^T[T].

It is worth mentioning that there is the existence of the universal Gaussian basis set (UGBS),⁶¹ an uncontracted basis set where all elements share the same basis set exponents, and thus, there is no difference between reference and target basis sets.

Alchemical forces are the mixed derivatives of the energy with respect to both nuclear charges and nuclear positions.⁴² Alchemical forces can be computed analytically either by differentiating the alchemical potential with respect to the nuclear coordinates or by differentiating the geometrical gradient with respect to the nuclear charges. In accordance with Schwarz’s theorem, both methods are valid and will be presented below.

The derivative of the one-electron density matrix $∂P∂R$ can be obtained through a CPHF calculation, independently for each molecular coordinate R, i.e., three times the number of atoms; for this reason, we think that the second approach is more convenient.

Differentiating the geometrical gradient with respect to the nuclear charge of an atom (Z_I) requires only one density derivative $∂P∂ZI$⁠, which is also needed for APDFT3 energy predictions.

The mixed derivative of the one-electron Hamiltonian [Eq. (21)] is only non-zero if the coordinate R is referred to the position of atom I.

All through our work, we used a locally modified version of PySCF⁶² in which we implemented subroutines for the analytical alchemical derivatives [Eqs. (6), (13), and (14)] as well as the alchemical force [Eqs. (25) and (26)]. If not specified differently, the pcX-2 basis set⁶³ will be used for second row elements in conjunction with the pc-2 basis set for hydrogen atoms. Basis sets’ coefficients and exponents were obtained from Basis Set Exchange.^64–66 Geometrical optimization in redundant internal coordinates was performed using a locally modified version of PyBerny,⁶⁷ where we included the transformation of the geometrical Hessian from Cartesian to internal coordinates. Root mean square deviations of atomic positions were evaluated using the program (RMSD)⁶⁸ implementing the Kabsch algorithm.⁶⁹ The molecular sketch of benzene’s B–N mutants was plotted using RDKit.^70–72 All Python code used in this work is made available, free of charge on a Zenodo repository.⁷³ 

Geometrical gradients and geometrical Hessians can be expanded as an alchemical series resulting in alchemical predictions of gradients and Hessians of the target molecules. The first alchemical prediction of the gradient is the alchemical force [Eqs. (25) and (26)].

Higher order derivatives of the gradient can be obtained through numerical differentiation of the alchemical force, while differentiation of the geometrical Hessian was only computed numerically. Derivatives were performed via a central finite difference stencil of seven points equispaced by Δλ = 0.1.

As a test case, we choose to analyze the bond stretching of some diatomic molecules and the stretching of a C–H bond in a methane molecule, and the results are plotted in Fig. 2.

The zeroth order gradient prediction is just the gradient of the reference, which is zero in its geometrical minimum, and this leads to a 100% prediction error for all molecules. Due to symmetry, the alchemical force for the transmutation N₂ → CO is zero, and therefore, the first order term does not improve the prediction. The inverse transmutation CO → N₂ has also a high first order error, which can be justified by similar reasons.

The third order prediction leads to an error of about 5%, while predictions of orders higher than 4 lead to an error, which is in any case below 1.50%. Going above the fifth order does not necessarily increase accuracy because of numerical errors.

Furthermore, if the derivatives of the atomic orbitals are neglected, the Taylor series converges to the gradient of the target calculated with the basis set of the reference. We can see this effect in the prediction of N₂ → CO, where the green line has an offset of ∼1.3%.

For Hessians, the zeroth order prediction, i.e., approximating the target Hessian with the reference’s one, has still a certain degree of accuracy. Additional terms increase the precision and going above the second order will lead to an error below 2%. The best predictions were obtained at the third and the fourth order, while after the fifth order term, the series diverge due to numerical errors in the calculation of the Hessians.

In a previous paper,¹⁸ we stated that one of the major sources of error in APDFT vertical energy predictions is the alchemical basis set error, which is explained in Sec. II D.

In the simple cases of diatomic molecules (BF ↔ CO, CO ↔ N₂) at the RHF level of theory, we tested performances of some of the most commonly used triple and quadruple ζ basis sets. We compared the family of Karlsruhe def2-nZ,^74,75 Dunning and co-worker’s correlation consistent polarized valence cc-pVnZ⁷⁶ and polarized core and valence cc-pCVnZ,⁷⁷ Jensen’s polarization consistent pc-n,^78–80 and their uncontracted variant optimized for x-ray spectroscopy pcX-n.⁶³ 

Figure 3 shows that $ΔEeqBS$ decreases with the increase in the basis set size, and this trend is approximately linear in the logarithmic plot for def2, cc-pCVnZ, and pc-n basis sets.

Core polarization is important for the alchemical invariance of basis sets: the cc-pVnZ basis functions, which are only polarized on the valence shell, have a higher $ΔEeqBS$ compared to other basis of equal size.

As a proof of this concept, we created a new basis set combining the core basis functions (the orbitals with s symmetry) of the cc-pCVQZ basis set with the valence and polarization orbitals of the smaller cc-pCVTZ. The orthonormality of this basis set labeled cc-pCV(sQ)TZ is guaranteed by spherical symmetry.

In Fig. 3, we can compare $ΔEeqBS$ for cc-pCVTZ and cc-pCV(sQ)TZ basis sets: passing from 6 to 8 s type orbitals per atom is enough to reduce the error by about one order of magnitude.

The alchemical basis set error in equilibrium bond lengths decreases with the basis set size, and the highest $ΔReqBS$ occurred using the cc-pVTZ basis set. Equilibrium bond lengths are consistently shorter as a consequence of the basis set superposition error.^81–83 

The other basis sets of the cc-pVnZ family (cc-pVQZ and cc-pV5Z) do not exhibit such behavior. $ΔReqBS$ is similar for cc-pCV(sQ)TZ and cc-pCVTZ, and both basis sets share the same valence and polarization functions. We can conclude that a good description of core and valence atomic orbitals is needed for accurate APDFT energies, while a good description of valence orbitals and the presence of polarization functions are required for accurate geometry predictions.

Best results in terms of accuracy per basis set number were achieved for both energy and geometry by the polarization consistent pcX-2 and pcX-3 basis sets. Those are uncontracted basis sets whose coefficients are optimized for x-ray spectroscopy. To effectively describe core-ionization of an atom, the basis functions used should provide a balanced representation of both the ground and core-excited state, where the effective nuclear charge changes by roughly +1. The same requirement is also true for alchemy, where the true nuclear charges change by ±1. This is a case where the computational solution to one specific problem (x-ray spectroscopy) can also be used effectively in a completely different application (APDFT). Acknowledged these conclusions on basis sets, we will use for the rest of the article the pCX-2 basis set, and since the pcX-2 basis sets are defined only for the second and third row elements, for hydrogens, we will use the pc-2 basis set.

In Sec. III B, we showed how alchemical perturbation can be used to predict geometrical gradients and Hessians. Those gradients and Hessians can be used to relax the geometry of the target molecule in several ways. We would like to show in this section that a geometrical relaxation using the alchemical predicted gradients and Hessians is possible and can be accurate depending on the APDFT order.

Morse potentials depend on four parameters: the bond dissociation energy D_e, the equilibrium distance R_e, the parameter a, which controls the width of the potential well, and V_e = V(r_e) the value of the energy at the minimum.

Knowing the energy, the first and second derivatives at a given point of the curve, we can match the number of parameters with the number conditions using an empirical value for the dissociation energy. An empirical dissociation energy can be used because the interpolated curves do not depend sensitively on D_e. In analogy to what is done in the UFF,⁸⁵ we will approximate D_e with the value of 100 kcal/mol times the total bond order [for the universal force field (UFF), the parameter is 80 kcal/mol].

For molecules bigger than diatomics, it is still possible to introduce the Morse potential interpolation as a correction to a Newton–Raphson step.

The mixed terms between different coordinates are still calculated in the paraboloid (Newton–Raphson) approximation.

A more detailed discussion of the Morse potential interpolation is given in the supplementary material of this paper.

Taking as an example the dissociation of CO, shown in Fig. 4, we can see that the quadratic approximation leads to a poor result when the starting point is far from the minimum. Since the geometries of the reference and target are usually substantially different, this condition occurs quite often in alchemical relaxation.

At the Morse potential, the whole energy curve of bond dissociations is described better, and this can lead to a more accurate relaxation, if the gradient and the Hessian are predicted accurately.

In the simple case of diatomic molecules, we can derive some considerations regarding the accuracy of the alchemical relaxation at different perturbation orders and compare the relaxations obtained using the Newton–Raphson step and the Morse potential interpolation.

Table I shows the results obtained with a NR step. In the predictions CO → BF and N₂ → BF, the error in bond lengths does not decrease with the alchemical prediction order. The reason for that is because part of the error comes from the relaxation step used and does not depend on how precisely gradients and Hessians are predicted.

Furthermore, in this method, the dependence of Hessian on the geometry is neglected, and the different geometries of the reference and target lead to a huge error in the prediction of harmonic vibrational frequencies. Table I shows that the Morse method performs better than the method of Newton–Raphson. The biggest errors were made in the predictions N₂ ↔ CO, where an alchemical change of ΔZ = ±2 is involved. All the other predictions, which have ΔZ = ±1, have an error in the bond length prediction on an average of 0.011 bohr at the third order, which decreases to 0.007 bohr at the fourth order. The error in the energy is also very low, 4.5 and 3.2 mHa for APDFT3 and APDFT4 predictions, respectively. The error in the harmonic frequency prediction is about 70 cm⁻¹; it is not accurate, but it is indicative of the value of the frequency.

Deprotonation and protonation energies are important quantities because they describe the changes in the enthalpy within an acid–base reaction and are needed in the calculation of the equilibrium constants pK_a and pK_b.^86–88 In this section, we demonstrate how APDFT can be used to predict protonation and deprotonation energies accurately. The APDFT3 error for vertical deprotonation energies can be as small as 1.4 kcal/mol,^26,27,41,89 and we would like to go beyond the fixed geometry approximation including the geometrical relaxation energy. More specifically, we exemplify our approach for the alchemical navigation of the iso-electronic ten electron series of second row hydrides CH₄ → NH₃ → H₂O → HF (see Fig. 5 for illustration). Moving from one molecule to a neighboring one, within alchemical deprotonation (protonation), a hydrogen nucleus is alchemically annihilated (created), while the nuclear charge of the heavy atom is simultaneously increased (decreased). The initial guess of the protonation site was chosen on the electrostatic potential minimum of the reference molecule; in that position, we placed a basis set for the hydrogen created. We implemented a three-point finite difference stencil (Δλ = 0.1); combining numerical and analytical differentiation, we obtained second order Hessian predictions, third order gradient prediction, and fifth order energy prediction. Combining numerical differentiation from a three-point finite difference scheme (Δλ = 0.1) with analytical differentiation [Eqs. (13), (14), and (25)], we calculated second order predictions for the Hessian, third order prediction for the gradient, and fifth order prediction for the energy. Geometry relaxation was performed using the Morse potential interpolation [Eqs. (31) and (32)].

Figure 6 shows the errors in protonations and deprotonations. The average error on energy predictions is 7 mHa for deprotonation and 5 mHa for protonation; this error is already present in the vertical predictions, as such does not come from the geometrical relaxations.

Alchemical deprotonation can predict relaxed bond lengths accurately; the results show a mean absolute error of 4 mbohrs, and in all cases, the error is lower than 10 mbohrs; for protonations, the predicted bond length of the hydrogens that were already present in the reference molecule is different from the one of the alchemically created hydrogen.

The effect of alchemical perturbation is much higher on the created hydrogens; for this reason, the error in bond length predictions is ten times higher than the error for hydrogens already present in the molecule, 56 and 4 mbohrs, respectively.

Angles are more flexible coordinates than bonds, and therefore, the prediction error of about 1° is quite substantial, and an exception is the prediction NH₃ → CH₄, where the smaller error might be due to a stiffer energy potential.

APDFT is efficient in predicting properties for multiple targets from a few derivatives when symmetric systems are chosen as the reference. An appealing application is the B–N doping of aromatic hydrocarbons, such as coronenes, fullerenes, or graphene.^28–30,42

Here, we show that it is possible to obtain approximate equilibrium energies and relaxed structures for all the 17 B–N doped mutants shown in Fig. 7 through the explicit calculation of just one CPHF derivative for benzene. The CPHF response matrix U can be obtained for all atoms via symmetry operations; in particular, we rotated the U matrix around the symmetry axis of the molecule using the Wigner D-matrix.⁹⁰ Energy prediction up to the third order can be obtained using Eqs. (6), (13), and (14); the same CPHF derivatives can be reused to calculate the alchemical force in Eq. (25). This is sufficient to give a first order prediction of the geometrical gradient for all 17 target molecules. The targets’ Hessians were approximated at the zeroth order with the Hessian of the reference. As shown in Fig. 2, this approximation has an error of ∼20% for the isoelectronic hydride series, which can still lead to qualitative improvements. In this case, the geometrical relaxation was performed using a Newton–Raphson step and not Morse because it is less sensitive to errors in gradient and Hessian. We measure geometry prediction errors through the root mean square deviation from the true target molecules’ minima (RMSD⁶⁹). It is important to remark on the energy ordering of the mutants. Every B–N substitution corresponds to a decrease in the energy of about 3 hartree, and therefore, stoichiometry is the dominating dimension. Among constitutional isomers, energy ordering depends on the relative positions of B and N: for monodoped mutants, the least stable is the one with B and N in the meta position, followed by the trans and ortho isomers. This ordering is analogous to the reactivity in electrophilic–nucleophilic aromatic substitutions.⁹¹ In fact, after a B–N doping of benzene, the electronic charge shifts from the boron atom to the nitrogen; in this context, nitrogen atoms act as Lewis acids (electrophiles) and borons as Lewis bases (nucleophiles).^92,93 We note, however, that these trends change upon removal of the nuclear repulsion terms, which do not affect the electronics. The predicted energy ordering is the same, both for the vertical and relaxed energies; the relaxation calculated only with the first order alchemical force is in fact not able to discriminate alchemical enantiomers⁹⁴ and does not change the energy order (Fig. 8).

Relaxation reduces the average RMSD by a factor of 60%: from 0.31 to 0.12 bohr. For singly doped mutants, the RMSD is reduced from 0.23 bohr to only 0.07 bohr.

The highest deviation from benzene’s geometry was found for mutants with adjacent N–N atoms and adjacent B–B atoms (mutants 4, 5 7, 15), the lowest for the mutants with alternating B–N (mutants 3, 14, 17). This is a consequence of alchemical symmetry:⁹⁴ at the first order, the alchemical force on a B–N bond is zero, and consecutively, there will be no stretching.

The mean absolute error in energy predictions is reduced by 70% after geometrical relaxation, passing from 89 to 27 mHa. In particular, for singly B–N doped mutants, it is reduced from 42 mHa to just 2 mHa.

In this work, we have shown that alchemical perturbation can successfully be applied to geometry relaxation and the prediction of equilibrium energies. A key step is the computation of the alchemical force $(∂2E∂Z∂R)$⁠, for which we gave an analytical formulation in the restricted Hartree–Fock case; an extension of the formula to other self-consistent field methods, such as UHF and DFT, should be possible. As we reported previously for vertical APDFT predictions,¹⁸ the basis set choice is crucial for accuracy. We confirm this conclusion is also for the APDFT based geometry relaxation and equilibrium energies studied in this work. More specifically, using the pcX-2 basis set, gradients and Hessians can be predicted with a fourth order error smaller than 2%. The consecutive geometrical relaxation using the employed Morse potential interpolation can predict equilibrium energies and bond lengths with a respective accuracy of ∼10 mHa and ∼0.01 bohr for small molecules (BF, CO, and N₂ and CH₄, NH₃, H₂O, and HF), and ∼27 mHa and ∼0.12 bohr for all the BN doped benzene mutants. The Morse potential was found to be a good descriptor of the bond dissociation energy profile for molecules within the range of bond distances encountered. Nevertheless, for different systems (e.g., solid state materials) or different bond lengths, the choice of another potential (e.g., Buckingham and Lennard-Jones) might be more suitable, but further investigations are beyond the scope of this paper.

In comparison to vertical APDFT predictions, we have confirmed the expectation that the inclusion of gradient and Hessian information results in a considerable increase in accuracy. The computational overhead is modest since first order APDFT predictions of gradients do not require more CPHF calculations than the one required for APDFT3. As for all alchemical methods, APDFT based geometry relaxations become particularly attractive when symmetry reduces number derivatives and simultaneously increases the number of target molecules. In particular, we demonstrated reasonable equilibrium energy and geometry predictions for all 17 BN doped mutant systems based on a single APDFT3 reference calculation for benzene. As an outlook, we believe that our findings and discussion would suggest that APDFT based geometry relaxation can be meaningful for many compounds on regular lattices, e.g., doped fullerenes or graphitic systems.

See the supplementary material for additional data on the accuracy of alchemical force predictions and for a more detailed discussion on the Morse potential interpolation.

We acknowledge support from the European Research Council (ERC-CoG grant QML and H2020 project BIG-MAP). This project has received funding from the European Union’s Horizon 2020 Research and Innovation Program under Grant Agreement Nos. 772834 and 957189. All results only reflect the authors’ views and the EU is not responsible for any use that may be made of the information. This work was partly supported by the NCCR MARVEL, funded by the Swiss National Science Foundation. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). The authors would also like to thank Dr. Max Schwilk for the helpful discussions about the formulation of analytic derivatives and Professor Frank Jensen from Aarhus University for the suggestion of using the pc-X basis set in the context of quantum alchemy.

The authors have no conflicts to disclose.

The data that support the findings of this study are openly available in Alchemical CPHF perturbator at http://doi.org/10.5281/zenodo.5606918.⁷⁰