Alchemical geometry relaxation

We propose to relax geometries throughout chemical compound space (CCS) using alchemical perturbation density functional theory (APDFT). APDFT refers to perturbation theory involving changes in nuclear charges within approximate solutions to Schr\"odinger's equation. We give an analytical formula to calculate the mixed second order energy derivatives with respect to both, nuclear charges and nuclear positions (named"alchemical force"), within the restricted Hartree-Fock case. We have implemented and studied the formula for its use in geometry relaxation of various reference and target molecules. We have also analysed the convergence of the alchemical force perturbation series, as well as basis set effects. Interpolating alchemically predicted energies, forces, and Hessian to a Morse potential yields more accurate geometries and equilibrium energies than when performing a standard Newton Raphson step. Our numerical predictions for small molecules including BF, CO, N2, CH$_4$, NH$_3$, H$_2$O, and HF yield mean absolute errors of of equilibrium energies and bond lengths smaller than 10 mHa and 0.01 Bohr for 4$^\text{th}$ order APDFT predictions, respectively. Our alchemical geometry relaxation still preserves the combinatorial efficiency of APDFT: Based on a single coupled perturbed Hartree Fock derivative for benzene we provide numerical predictions of equilibrium energies and relaxed structures of all the 17 iso-electronic charge-netural BN-doped mutants with averaged absolute deviations of $\sim$27 mHa and $\sim$0.12 Bohr, respectively.


I. INTRODUCTION
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 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 works [13][14][15] succeeded in machine learning the optimized molecular geometries, but most QML methods do not predict directly geometrical structures.
APDFT 16 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 suba) Electronic mail: anatole.vonlilienfeld@univie.ac.at space; in this regard, APDFT is more versatile because it requires the explicit QM 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 16 . APDFT accuracy depends on the perturbation order and on the basis set choice 17 . If a high-quality calculation (e.g. CCSD) is used as reference, the error in third-order energy prediction can be smaller than the error made using a cheaper QM method such as HF and MP2. 16,17 Alchemical changes and derivatives were at first introduced as a way to rationalize energy differences between molecules. [18][19][20] Successive works showed that APDFT is also capable of predicting many other molecular properties, such as electron densities, dipole moments 16 , ionization potentials 21,22 , HOMO energy eigenvalues 23 , band structures 24 , deprotonation energies 25,26 . 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. [27][28][29] Using pseudopotentials, quantum alchemy can also be applied to crystal structures 24,[30][31][32][33][34] , and can even suggest catalyst improvements [35][36][37][38][39] .
One of the major remaining challenges in APDFT is to include geometry relaxation as the equilibrium molecular structure changes from reference to target; this would constitute a huge improvement over fixed atoms "vertical" predictions. To the best of our knowledge, only two papers 23,40 have tried to predict equilibrium energies, applying APDFT to "non-vertical" energy predictions (predicting energies for a perturbation that includes both changes in nuclear charges and in geometry).
In this paper we propose an alternative approach: we predict the geometry energy derivatives (gradients and Hessians) for the target molecule, and subsequently, use them to relax the geometry (see Fig.1).
FIG. 1. APDFT predictions of gradients and Hessians enable efficient geometry relaxation of target molecules. The predictive accuracy improves with Morse potential interpolation.
Alchemical forces are defined as the second-order mixed derivatives of the energy with respect to nuclear positions and nuclear charges ( ∂ 2 E ∂Z∂R ). Alchemical forces appear as non-diagonal terms in the generalized Hessian 41 ( 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 [42][43][44][45] ( ∂ 2 E ∂Ne∂R ), alchemical Fukui functions 21,22,[46][47][48][49] ∂Ne∂Z ) Those functions have already been studied, while to the best of our knowledge no paper 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 UHF or DFT is possible, and some of the conclusions derived for RHF may also hold for other levels of theory.
In the first section of this paper, 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 the second section we show that it is possible to compute higherorder derivatives through numerical differentiation and those derivatives can be used to predict geometrical gra-dients 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.

A. Alchemical Perturbation Density Functional Theory
Within the Born Oppenheimer approximation 50 , the non-degenerate electronic ground state wave function Ψ 0 (r; R, Z, N e ) depends parametrically on the positions of the nuclei R, the nuclear charges Z, and the number of electrons N e . The nuclear charge vector Z specifies the chemical composition of a given molecule, for real molecules Z can only have integer values, though from a theoretical point of view it is possible to extend the concept of atomic charge to any real number. It is thus possible to define a smooth path for Z coordinates that connects molecules of different chemical composition 21 . We call "alchemical transmutation" the process of transforming one element into another for one or more atoms in a molecule. In order to connect a reference molecule R to a target molecule T we can define the alchemical coordinate λ as a linear transformation of the nuclear charges from the reference Z R = Z(λ = 0) to the target Z T = Z(λ = 1).
The linear transformation of the nuclear charges corresponds also to a linear transformation of the electronic Hamiltonian and the nuclear-electron attraction opera-torsV ne : For isoelectronic transmutations, the total electronic ground state energy of the molecule is continuous and differentiable in λ. 23 From the energy of the reference molecule E R = E(λ = 0) and its derivatives, we can express the energy of the target molecule E T = E(λ = 1) through a Taylor series expansion, where ∆λ = 1 and where we assume convergence.
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.
For a variational wave function, the first alchemical derivative of the electronic energy can be evaluated via the Hellmann-Feynman theorem 23,51 : For a monodeterminantal wave function (e.g. the solution of an Hartree Fock or a Kohn-Sham DFT calculation) the first alchemical derivative (Equation 4) can be rewritten as a contraction between the one electron density matrix P of the reference molecule with the change in the nuclear-electron potential energy operator ∆V : Higher order derivatives can be obtained via numerical differentiation or analytically using either Møller Plesset perturbation theory 40 or a coupled perturbed approach (CPHF or CPKS-DFT) 52 . The contribution of nuclear nuclear repulsion to the alchemical derivatives can be obtained analytically as the derivative of a classical charge distribution.

B. Coupled Perturbed Hartree Fock
Coupled perturbed methods are a powerful tool to calculate derivatives properties of a mono determinant wave function (HF,DFT). They represent the standard way to calculate polarization [53][54][55] and can also be used for alchemical perturbations. 27,29,47,56 The mathematical treatment to calculate the molecular response to an external electric field or to the electric field generated by the transmutation of one atom is indeed identical; here we would like to recall some theoretical aspects of the coupled perturbed Hartree Fock method.
In CPHF the first derivative of the coefficient's matrix C is written as a the product of C and a unitary response matrix U : From the orthogonality condition C T SC = I it is possible to show that U must be skew symmetric.
Since U acts as a rotation matrix, and since orbital rotations inside the occupied-occupied and the virtualvirtual blocks do not affect the total energy, we can choose the o − o and the v − v blocks of U to be zero; keeping only the o − v blocks of U non zero. The derivative of the Roothaan equation F C = SC with respect to λ becomes: Multiplying from the left by C T and simplifying using the identities C T SC = I and C T F C = , leads to: Since is a diagonal matrix and U, U , U are zero on the diagonal blocks we can divide Equation 9 in two parts: 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 λ, µ, ν, σ as indexes of atomic orbitals and capital letters I, J, k as atoms' indexes.
The partial derivative of the Fock matrix in Equation  9 contain the term ∂V /∂λ, but also an implicit dependence on U since F depends on C. A detailed solution to the CPHF equations can be found in Refs. 52,53,57 After solving Eq. 9 derivatives of the one particle density matrix can be evaluated as:

C. Second and third energy derivatives from atomic contributions
In APDFT any molecular property can be expressed as a function of the alchemical coordinate λ. Explicating the dependence on the nuclear charges (A(λ) = A (Z(λ))), differentiation with respect to λ can be performed using the chain rule: The last equation comes directly from Eq.1 . 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 .
Consistently with the Wigner 2n + 1 rule, from these first order derivatives it is possible to evaluate the second and the third derivatives of the electronic energy: This formulation is particularly convenient for highly symmetrical systems, where many of the derivatives ∂A/∂Z I are symmetrically equivalent. Only few explicit CPHF calculations are needed, while the number of possible targets grows combinatorially with the size of the system.

D. Basis set energy correction
We pointed out in a previous paper 17 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: In the same paper we proposed a correction to the basis set error that can be added to any order APDFT energy prediction. The correction decomposes the total basis set error into a sum of individual atom contributions calculated from isolated atoms.
It is worth to mention the existence of the universal Gaussian basis set (UGBS) 58 , an uncontracted basis set where all elements share the same basis sets exponents, thus there is no difference between reference and target basis sets.

E. Alchemical forces
Alchemical forces are the mixed derivatives of the energy with respect to both nuclear charges and nuclear positions 41 . Alchemical forces can be computed analytically either by differentiating the alchemical potential with respect to the nuclear coordinates or 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.

First formulation
Differentiating the first alchemical derivative of the electronic energy (Eq.5) with respect to the nuclear coordinates R I of the atom I leads to: where: 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 is more convenient following the other approach.

Second formulation
Differentiating the geometrical gradient expression with respect to the nuclear charge of an atom (Z I ) requires only one density derivative ∂P ∂Z I , which is also needed for APDFT3 energy predictions.
The first derivative of the energy with respect to one nuclear Cartesian coordinate R is: 57 Where H (1) µν is the mono electronic part of the Hamiltonian,(µλ||νσ) is the sum of the bielectronic Coulomb and Exchange integrals, and W is the energy weighted density matrix: Eq.17 is composed of three terms, and we need to differentiate all of them with respect to the nuclear charge Z I of a given atom I.
For the first term of Eq.17: The mono electronic Hamiltonian H (1) is composed of two parts: the kinetic energy operator T , which is independent of Z I , and the nuclear electron potential V .
The mixed derivative of the one electron Hamiltonian (Eq.20) is only non zero if the coordinate R is referred to the position of atom I.
For the second term of Eq.17: Differentiating the third term of Eq.17 yields: The derivative of W can be obtained from Equation 18 using the response matrix U and the derivatives of the MO energies ∂ /∂Z I from the solution of the CPHF equation (Eq.9): Collecting together all terms, the electronic part of the alchemical force at a RHF level of theory can than be expressed as follows:

Nuclear-nuclear contribution
The nuclear nuclear contribution to the alchemical force can be evaluated analytically as derivatives of an electrostatic charge distribution: III. NUMERICAL RESULTS

A. Computational details
All through our work we used a locally modified version of PySCF 59 in which we implemented subroutines for the analytical alchemical derivatives (Eqs. 5,12,13), as well as the alchemical force (Eqs. 24,25). If not specified differently will be used the pcX-2 basis set 60 for second row elements in conjunction with pc-2 basis set for hydrogen atoms. Basis sets' coefficients and exponents were obtained from Basis Set Exchange 61-63 . Geometrical optimization in redundant internal coordinates was performed using a locally modified version of PyBerny 64 , 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" 65 implementing the Kabsch algorithm 66 . Molecular sketch of benzene's B-N mutants were plotted using RDKit. 67-69 All Python code used in this work is made available, free of charge on a Zenodo repository. 70

B. Convergence of gradients and Hessian series
Geometrical gradients and geometrical Hessians can be expanded as an alchemical series resulting into alchemical predictions of gradients and Hessians of the target molecules. The first alchemical prediction of the gradient is the alchemical force (Eqs. 24,25).
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 7 points equispaced by ∆λ = 0.1.
As a test case we choose to analyse the bond stretching of some diatomic molecules and the stretching of a C-H bond in a methane molecule, the results are plotted in Figure 2.
The zeroth order gradient prediction is just the gradient of the reference, which is zero in its geometrical minimum, this leads to a 100% prediction error for all molecules. Due to symmetry the alchemical force for the transmutation N 2 →CO is zero, therefore the first order term does not improve the prediction. The inverse transmutation CO→N 2 has also an 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 5 th 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 2 →CO where the green line has an offset of approximately 1.3%.
For Hessians the 0 th order prediction, i.e. approximating the target Hessian with the reference's one, have still a certain degree of accuracy. Additional terms increase the precision and going above second order will lead to an error below 2%. The best predictions were obtained at the third and fourth and order, while after the fifth order term the series diverge due to numerical errors in the calculation of the Hessians.

C. Relaxed basis set errors
In a previous paper 17 we stated that one of the major sources of error in APDFT vertical energy predictions is the alchemical basis set error, explained in section II D.
For the sake of this paper is of great interest to analyse the alchemical basis set error error, not only for vertical energy predictions but also for the predictions of energies and geometries of the target in its energy minimum "relaxed errors" (∆E BS eq , ∆R BS eq ).
In the simple cases of diatomic molecules (BF↔CO, CO↔N 2 ) at the RHF level of theory, we tested performances of some of the most commonly used the triple and quadruple ζ basis set. We compared the family of Karlsruhe def2-nZ 71,72 , the Dunning and coworker's correlation consistent polarized valence cc-pVnZ 73 and polarized core and valence cc-pCVnZ 74 , the Jensen's polarization consistent pc-n 75-77 and their uncontracted variant optimized for X-ray spectroscopy pcX-n 60 .

Figure 3 shows that ∆E BS
eq decrease with the increase in the basis set size, 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 ∆E BS eq 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 Figure 3 we can compare ∆E BS eq for cc-pCVTZ and cc-pCV(sQ)TZ basis sets: passing from 6 to 8 s type orbitals per atoms is enough to reduce the error by about one order of magnitude.
The cc-pVTZ basis set has a large geometry error, that's because the systems, in order to compensate for the high energy error, increase the orbital overlap between the atoms shortening the bond lengths. The other basis sets of the cc-pVnZ family (cc-pVQZ and cc-pV5Z) don't suffer as much from this problem. ∆R BS eq is similar for cc-pCV(sQ)TZ and cc-pCVTZ, both basis set 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 in valence orbitals and the presence of polarization functions is required for accurate geometry predictions.
Best results in terms of accuracy per basis sets number were achieved both for 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 describe effectively a coreionization 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-rays 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 second and third row elements, for hydrogens we will use the pc-2 basis set.

IV. APPLICATIONS TO GEOMETRY RELAXATION
In section 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 it is possible and can be accurate depending on the APDFT order.

A. Relaxation techniques
The simplest relaxation technique is a single step in the Newton-Raphson optimization method. The potential energy surface is approximated with a paraboloid; given a a set of Cartesian coordinates x, a gradient g and a Hessian H, the optimization step and the relaxation energy are: A better representation for the monodimensional potential of a bond stretching was given by Morse 78 : 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. This is legit because the the interpolated curves do not depend sensibly on D e , and analogously to what is done in the UFF 79 we will approximate D e with the value of 100kCal/mol times the total bond order (for UFF this parameter is 80kCal/mol).
For molecules bigger than diatomics is still possible to introduce the Morse potential interpolation as a correction to a Newton-Raphson step.
Using internal coordinates q, we calculated the correction treating every bond if b if it were an independent coordinate: The mixed terms between different coordinates are still calculated in the paraboloid (Newton Raphson) approximation.
The Morse potential correction affect also the energy: A more in detail discussion of the Morse potential interpolation is given in the supplementary materials of this paper. Taking as an example the dissociation of CO, shown in figure 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 reference and target are usually substantially different, this condition occurs quite often in alchemical relaxation.
At the Morse potential describes the whole energy curve of bond dissociations better, and this can lead to a more accurate relaxation, if the gradient and the Hessian are predicted accurately.

B. Diatomics
In the simple case of diatomic molecules we can derive some considerations regarding the accuracy of the alchemical relaxation at different perturbation order, and compare the relaxations obtained using Newton Raphson step and the Morse potential interpolation.
Tab.I shows the results obtained with a NR step. In the predictions CO→BF and N 2 →BF the error in bond lengths is consistent, and a higher alchemical order does not give better results.
Furthermore in this method is neglected the dependence of Hessian on the geometry, the different geometries of reference and target lead to a huge error in the prediction of harmonic vibrational frequencies.
Tab.I shows that the Morse method performs better than Newton Raphson's. The biggest errors were made in the predictions N 2 ↔CO where is involved an alchemical change of ∆Z = ±2. All the other predictions which have ∆Z = ±1 have an error in the bond length prediction on average of 0.011 Bohr at third order which decreases to 0.007 Bohr at fourth order. The error in the energy is also very low, 4.5 mHa and 3.2 mHa for APDFT3 and APDFT4 predictions respectively. The error in the harmonic frequencies prediction is about 70 cm −1 ; is not accurate, but it is indicative of the value of the frequency.

C. Alchemical protonation and deprotonation
Deprotonation and protonation energies are important quantities, because describe the changes in enthalpy within an acid-base reactions, and are needed in the calculation of the equilibrium constants pK a and pK b . [80][81][82] In this section, we demonstrate how APDFT can be used to predict protonation and deprotonation energies accurately. APDFT3 error for vertical deprotonation energies can be as small as 1.4 kcal/mol 25,26,40,83 , 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 10 electron series of second row hydrides CH 4 →NH 3 →H 2 O→ HF,(see figure 5 for illustration). Moving from one molecule to a neighbouring 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 chose 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 points finite difference scheme (∆λ = 0.1) with analytical differentiation (Eqs.12, 13, 24), 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. 30,31). 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; results show a mean absolute error of 4 milli Bohr and in all cases, the error is lower than 10 milli Bohr, 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 10 times higher than the error for hydrogens already present in the molecule, 56 and 4 milli Bohr respectively.
Angles are more flexible coordinates than bonds, therefore the prediction error of about 1 • is quite substantial, an exception is the prediction NH 3 → CH 4 where the smaller error might be due to a stiffer energy potential.

D. B-N Doping of Benzene
APDFT is efficient in predicting properties for multiple targets from a few derivatives when symmetric systems are chosen as reference. An appealing application is the B-N doping of aromatic hydrocarbons such as coronenes, fullerenes, or graphene. [27][28][29]41 FIG. 7. A single APDFT3 reference calculation for benzene affords equilibrium energy predictions for all seventeen possible iso-electronic charge-neutral B-N doped mutants at Hartree-Fock level of theory.
Here, we show that 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 84 . Energy prediction up to third-order can be obtained using Eqs.5, 12, 13, the same CPHF derivatives can be reused to calculate the alchemical force in Eq.24. 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 approximately 20% for the isolectronic 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 sensible to errors in gradient and Hessian. We measure geometry prediction errors through the root mean square deviation from the true target molecules' min-ima (RMSD 66 ). It is important to remark on the energy ordering of the mutants. Every B-N substitution corresponds to a decrease in energy of about 3 Hartree, therefore stoichiometry is the dominating dimension. Among constitutional isomers, energy ordering depends on the relative positions of B and N; for mono doped mutants, the least stable is the one with B and N in meta position, followed by the trans and the ortho isomers. This ordering is analogous to the reactivity in electrophilicnucleophilic aromatic substitutions 85 . 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). 86,87 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 the relaxed energies, the relaxation calculated only with the first order alchemical force is in fact not able to discriminate alchemical enantiomers 88 and does not change the energy order. Relaxation reduces the average RMSD by a factor of 60% : from 0.31 Bohr 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 88 : at 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 mHa to 27 mHa. In particular, for singly B-N doped mutants it is reduced from 42 mHa to just 2 mHa.

V. CONCLUSION
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 ( ∂ 2 E ∂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 17 , the basis set choice is crucial for accuracy. We confirm this conclusion also for the APDFT based geometry relaxation and equilibrium energies studied in this work. More specifically, using the pcX-2 basis set, gradients, and Hessian can be predicted with a 4 th 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, N 2 , and CH 4 , NH 3 , H 2 O, and HF), and ∼27 mHa and ∼0.12 Bohr for all the BN doped benzene mutants.
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 number of target molecules. In particular, we demonstrated reasonable equilibrium energy and geometry predictions for all the seventeen 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.

VI. ACKNOWLEDGEMENT
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 #772834 and #957189. All results only reflect the authors' view 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 Prof. Frank Jensen from Aarhus University for the suggestion of using the pc-X basis set in the context of quantum alchemy.