Fitting quantum machine learning potentials to experimental free energy data: predicting tautomer ratios in solution

Abstract

The computation of tautomer ratios of druglike molecules is enormously important in computer-aided drug discovery, as over a quarter of all approved drugs can populate multiple tautomeric species in solution. Unfortunately, accurate calculations of aqueous tautomer ratios—the degree to which these species must be penalized in order to correctly account for tautomers in modeling binding for computer-aided drug discovery—is surprisingly difficult. While quantum chemical approaches to computing aqueous tautomer ratios using continuum solvent models and rigid-rotor harmonic-oscillator thermochemistry are currently state of the art, these methods are still surprisingly inaccurate despite their enormous computational expense. Here, we show that a major source of this inaccuracy lies in the breakdown of the standard approach to accounting for quantum chemical thermochemistry using rigid rotor harmonic oscillator (RRHO) approximations, which are frustrated by the complex conformational landscape introduced by the migration of double bonds, creation of stereocenters, and introduction of multiple conformations separated by low energetic barriers induced by migration of a single proton. Using quantum machine learning (QML) methods that allow us to compute potential energies with quantum chemical accuracy at a fraction of the cost, we show how rigorous relative alchemical free energy calculations can be used to compute tautomer ratios in vacuum free from the limitations introduced by RRHO approximations. Furthermore, since the parameters of QML methods are tunable, we show how we can train these models to correct limitations in the underlying learned quantum chemical potential energy surface using free energies, enabling these methods to learn to generalize tautomer free energies across a broader range of predictions.

---

Introduction

The most common form of tautomerism, prototropic tautomerism describes the reversible structural isomerism involving the sequential processes of bond cleavage, skeletal bond migration and bond reformation in which a hydrogen is transferred.¹ Numerous chemical groups can show prototropic tautomerism. Common examples include keto–enol (shown in Fig. 2), amide/imidic acid, lactam/lactim, and amine/imine tautomerism.²

Tautomerism influences many aspects of chemistry and biology

Tautomerism adds a level of mutability to the static picture of chemical compounds. The change in the chemical structure between different tautomeric forms of a molecule is accompanied by changes in physico-chemical properties. By virtue of the movement of a single proton and the rearrangement of double bonds, tautomerism can significantly alter a molecule's polarity, hydrogen bonding pattern, its role in nucleophilic/electrophilic reactions, and a wide variety of physical properties such as partition coefficients, solubilities, and pK_a.^3,4

Tautomerism can also alter molecular recognition, making it an important consideration for supramolecular chemistry. To optimizing hydrogen bond patterns between a ligand and a binding site tautomerism has to be considered. In a theoretical study of all synthetic, oral drugs approved and/or marketed since 1937, it has been found that 26% exist as an average of three tautomers.⁵ While tautomerism remains an important phenomena in organic chemistry, it has not gained much appreciation in other scientific fields.

The typical small free energy difference between tautomers poses additional challenges for protein-ligand recognition. Local charged or polar groups in the protein binding pocket can shift the tautomer ratio and result in dominant tautomers in complex otherwise not present in solution.⁶ For this reason the elucidation of the dominant tautomer in each environment is not enough—without the knowledge about the ratio (i.e. the free energy difference) between tautomeric forms in the corresponding phase a correct description of the experimental (i.e. macroscopic) binding affinity might not be possible. To illustrate this, one might think about two extreme cases: one in which the tautomeric free energy difference in solution is 10 kcal mol⁻¹ and one in which it is 1.0 kcal mol⁻¹. It seems unlikely that the free energy difference of 10 kcal mol⁻¹ will be compensated by the protein binding event (therefore one could ignore the unlikely other tautomer form), but a tautomeric free energy difference of 1 kcal mol⁻¹ could easily be compensated. Examples for this effect—changing tautomer ratios between environments—are numerous for vacuum and solvent phase. One example is the neutral 2-hydroxypyridine (2-HPY)/2-pyridone (2-PY) tautomer. The gas phase internal energy difference between the two tautomers is ∼0.7 kcal mol⁻¹ in favor of the enol form (2-HPY), while in water an internal energy difference of 2.8 kcal mol⁻¹ was reported in favour of 2-PY. In cyclohexane, both tautomer coexist in comparable concentration while the 2-PY is dominant in solid state and polar solvents.⁷

The example above assumed two tautomeric forms—if there were multiple tautomeric forms, each with small free energy differences, using a single dominant tautomer in solvent and complex may still represent only a minor fraction of the true equilibrium tautomer distribution.

The study of tautomer ratios requires sophisticated experimental and computational approaches

The experimental study of tautomer ratios is highly challenging.⁸ The small free energy difference and low reaction barrier between tautomers—as well as their short interconversion time—can make it nearly impossible to study specific tautomers in isolation. Recently, efforts have been made to collect some of the experimental data and curate these in publicly available databases, specifically Wahl and Sander⁹ and Dhaked et al.¹⁰

In the absence of a reliable, cheap, and fast experimental protocol to characterize tautomeric ratios for molecules, there is a need for computational methods to fill the gap. And—even if such a method exists—predicting tautomer ratios of molecules that are not yet synthesized motivates the development of theoretical models. Computational methods themselves have a great need for defined tautomer ratios. Most computational methods use data structures in which bond types and/or hydrogen positions have to be assigned a priori and remain static during the calculation, introducing significant errors if an incorrect dominant tautomer is chosen.¹¹

The third round of the Statistical Assessment of the Modeling of Proteins and Ligands (SAMPL2) challenge included the blind prediction of tautomer ratios and provided an interesting comparison between different computational methods.¹² Most of the 20 submissions were using implicit solvent models and ab initio or DFT methods in combination with a thermodynamic cycle (as shown in Fig. 1) to evaluate tautomeric free energy difference.¹² However, as stated in Martin,¹³ “In summary, although quantum chemical calculations provide much insight into the relative energies of tautomers, there appears to be no consensus on the optimal method”. For the 20 tautomer pairs investigated in the SAMPL2 challenge (this includes 13 tautomer pairs for which tautomer ratios were provided beforehand) the three best performing methods reported an root mean squared error (RMSE) of 2.0,¹⁴ 2.5 (ref. 15) and 2.9 kcal mol⁻¹.¹⁶ While these results are impressive it also shows that there is clearly room for improvement. And, maybe even more importantly, the SAMPL2 challenge showed the need for investigating a wider variety of tautomer transitions to draw general conclusions about best practices and methods for tautomer predictions.¹² The excellent review of Nagy concluded that tautomer relative free energies are sensitive to “the applied level of theory, the basis set used both in geometry optimization and […] single point calculations, consideration of the thermal corrections […] and the way of calculating the relative solvation free energy”.¹⁷ We'd like to add “quality of 3D coordinates” to these issues and discuss them in the following.

Fig. 1 Tautomeric free energy differences in solution are typically calculated using a thermodynamic cycle, independent of the level of theory that is used to obtain the individual terms. A typical quantum chemistry protocol calculates the free energy in aqueous phase as the sum of the free energy in gas phase (; calculated using the ideal gas RRHO approximation) and the standard state transfer free energy (; obtained using a continuum solvation model). The tautomeric free energy difference in solution is then calculated as the difference between for two tautomers.

Standard approaches to quantum chemical calculations of tautomer ratios introduce significant errors

The tautomeric free energy difference in solution can be calculated from the standard-state Gibbs free energy in aqueous phase of the product and educt of the corresponding tautomer reaction, which itself is calculated as the sum of the gas-phase standard-state free energy and the standard-state transfer free energy , shown as a thermodynamic cycle in Fig. 1.

The standard-state Gibbs free energy is calculated using a quantum chemistry estimate for the electronic energy and, based on its potential energy surface, a statistical mechanics approximation to estimate its thermodynamic properties. While the accuracy of the electronic energy and the transfer free energy ΔG_S is dependent on the chosen method and can introduce significant errors in the following we want to concentrate on the approximations made by the thermochemical correction to obtain the gas phase free energy: the rigid rotor harmonic oscillator (RRHO) approximation and the use of single or multiple minimum conformations to generate a discrete partition function (as written in eqn (1)). In the following the standard-state of all the thermodynamic properties should be implied and will not be added to the notation.

A commonly used approach for the zero point energy (ZPE) and thermal contributions is based on the ideal gas RRHO, assuming that the molecule is basically rigid and its internal motions comprise only vibrations of small amplitude where the potential energy surface can be approximated as harmonic around a local energy minimum. This assumption leads to an analytical approximation that describes the local configuration space around a single minimum based on the local curvature of the potential energy surface (Fig. 2).¹⁸

Fig. 2 The rigid-rotor harmonic oscillator (RRHO) approach constructs a partition function using the curvature of the potential energy landscape at a minimum modelling all bonded terms as harmonic. The enol and keto form of a molecules (in the enol form: 1-(cyclohex-3-en-1-yl)ethen-1-ol) is shown and the main conformational degrees of freedom highlighted. The middle panel shows how the use of a local and harmonic approximation to the partition function can approximate the overall potential energy landscape, if all relevant minimum conformations are enumerated. The lower panel shows the probability density resulting from the harmonic approximations, with solid colored regions the difference between the true and approximated probability density resulting from anharmonicity and/or bonded terms that would better be modeled using a hindered/free rotor.

Errors are introduced because the harmonic oscillator approximation assumes that, for each normal mode, the potential energy associated with the molecule's distortion from the equilibrium geometry has a harmonic form. Especially low-frequency torsional modes would be more appropriately treated as hindered internal rotations at higher temperatures. This can lead to a significant underestimate of the configurational entropy, ignoring the contributions of multiple energy wells.¹⁹ The correct treatment of such low-frequency modes in the analytical rotational entropy part of the RRHO partition function can add considerable computational cost, since the numerical solution requires the calculation of the full torsion potential (periodicity and barrier heights).²⁰

Another, related shortcoming of the RRHO approach is the focus on a single minimum conformation. Assuming that the correct global minimum has been identified this approach ignores the configurational entropy of all other conformations to the partition function. Methods like the ‘mining minima’ approach can help to mitigate this problem and construct a partition function using local configurational integrals from multiple minimum energy wells.¹⁹ The standard state free energy G° of a molecule can then be calculated for each of the minimum conformations k separately and combined as follows

to obtain a weighted free energy.^21,22 The free energy obtained in eqn (1) mitigates the shortcomings of the RRHO approach to model torsions with periodicity of 1 since—as long as the different minima of the torsion are generated as conformations—these are included in the sum over all conformations. But the configurational integral is still restricted to the minimum conformation and the contribution of its energy well, ignoring anharmonicity in the bonded terms and contributions from conformations outside the minimum energy wells.

Alchemical relative free energy calculations with machine learning potentials can compute true tautomeric free energy differences, including all classical statistical mechanical effects

Limitations of the ideal gas RRHO approximation to the free energy and ZPE correction, challenges in enumeration of local minimum conformations, and the consistent treatment of internal/external rotational symmetry, as well as approximations in the continuum solvent model can lead to errors in the free energy that are difficult to detect and correct. The use of molecular dynamics (MD) simulations, explicit solvent molecules, and a rigorous classical treatment of nuclear motion to sample independent conformations from an equilibrium distribution can overcome the above mentioned challenges. To describe tautomer free energies in solution this would require computationally prohibitive MM/QM calculations in which the molecule of interest is treated quantum mechanically (QM) and the solute treated using a molecular mechanics (MM) description.

One of the most exciting developments in recent years has been the introduction of fast, efficient, accurate and transferable quantum machine learning (QML) potentials (e.g. ANI,²³ PhysNet,²⁴ and SchNet²⁵) to model organic molecules. QML potentials can be used to compute QM-based Hamiltonians and—given sufficient and appropriate training data—are able to reproduce electronic energies and forces without loss of accuracy but with orders of magnitude less computational cost than the QM methods they aim to reproduce. QML potentials have been successfully applied to molecular dynamics and Monte Carlo (MC) simulations.²⁶ Here, we present the application of a machine learning potential for the calculation of alchemical free energies for tautomer pairs in the gas phase (Fig. 3).

Fig. 3 Alchemical relative free energy calculations with quantum machine learning (QML) potentials like ANI can rigorously compute free energy differences. Linearly interpolating between two potential energy functions using the alchemical coupling parameter λ enables the sampling of bridging distributions between the physical endstates representing the two tautomers. Here, noninteracting dummy atoms—indicated by the empty circle at each endstate—are included to facilitate this transformation. In this work, we present the application of this concept to calculate relative free energies of tautomer pairs in vacuum.

We begin by investigating the accuracy of a current, state of the art approach to calculate tautomer ratios on a set of 468 tautomer pairs selected from the publicly available Tautobase dataset⁹ spanning different tautomer reactions, number of atoms and functional groups. We are using a popular DFT functional (B3LYP) and basis sets (aug-cc-pVTZ and 6-31G(d))^27–31 for the calculation of the electronic energy and a continuum solvation model (SMD³²) to model the transfer free energy.

Furthermore we are investigating the effect of a more rigorous statistical mechanics treatment of the gas phase free energies. To investigate this effect we calculate tautomeric free energy difference in vacuum Δ_tG^calc_vac using (1) alchemical relative free energy calculations and (2) multiple minimum conformations in combination with the RRHO approximation (as used in the quantum chemistry approach). To enable a direct comparison between the two approaches we are using a QML potential in both calculations.

Since we are interested in tautomeric free energy differences in solution we investigate the possibility to optimize the QML parameters to include crucial solvent effects. The framework we have developed to perform alchemical relative free energy calculations enables us to obtain a relative free energy estimate that can be optimized with respect to the QML parameters. We are using a small training set of experimentally obtained tautomer free energies in solution Δ_tG^exp_solv and importance sampling to obtain reweighted with the optimized parameters.

Here we report the first large scale investigation of tautomer ratios using quantum chemical calculations and machine learning potentials in combination with alchemical relative free energy calculations spanning a large chemical space and Δ_tG^exp_solv values.

Standard quantum chemistry methods predict tautomer ratios with an RMSE of 3.1 kcal mol⁻¹ for a large set of tautomer pairs

Calculating tautomeric free energy differences in solution Δ_tG^calc_solv is a challenging task since its success depends on a highly accurate estimate of both the intrinsic free energy difference between tautomer pairs as well as the transfer free energy.

We calculated the tautomeric free energy difference Δ_tG^calc_solv for 468 tautomer pairs (460 if the basis set 6-31G(d) was used) and compared the results with the experimentally obtained tautomer ratios in solution (expressed as free energy difference in solution Δ_tG^exp_solv) in Fig. 4.

Fig. 4 State of the art quantum chemistry calculations are able to calculate tautomeric free energy differences Δ_tG^calc_solv with a RMSE of 3.1 kcal mol⁻¹. The direction of the tautomer reaction is chosen so that the experimentally obtained tautomeric free energy difference Δ_tG^exp_solv is always positive. Panel (A) shows Δ_tG^calc_solv as the difference between the sum of the gas phase free energy and transfer free energy for each tautomer pair plotted against the experimental tautomeric free energy difference in solution Δ_tG^exp_solv. B3LYP/aug-cc-pVTZ is used for the gas phase geometry optimization and single point energy calculation, the ideal gas RRHO approximation is used to calculate the thermal corrections. The transfer free energy is calculated on B3LYP/aug-cc-pVTZ/SMD optimized geometries using B3LYP/6-31G(d) and SMD. Values in quadrant II (x-axis entries positive, y-axis entries negative) indicate calculations that assigned the wrong dominant tautomer species (different sign of Δ_tG^calc_solv and Δ_tG^exp_solv). The dashed line indicates the ideal behavior of the calculated and experimental values, the grey lines mark the ±1 kcal mol⁻¹ interval. Red dots indicate tautomer pairs with more than 10 kcal mol⁻¹ absolute error between Δ_tG^calc_solv and Δ_tG^exp_solv. These tautomer pairs are separately shown in Table S.I.1.† In panel (B), the top panel shows the kernel density estimate (KDE) and histogram of Δ_tG^calc_solv and Δ_tG^exp_solv. The red line indicates zero free energy difference (equipopulated free energies). In the lower panel the KDE of the difference between Δ_tG^exp_solv and Δ_tG^calc_solv is shown. MAE and RMSE are reported in units of kcal mol⁻¹. Quantities in brackets [X;Y] denote 95% confidence intervals. The Kullback–Leibler divergence (KL) was calculated using KL(Δ_tG^exp_solv‖Δ_tG^calc_solv).

Three quantum chemistry approaches were used to calculate these values:

• B3LYP/aug-cc-pVTZ/B3LYP/aug-cc-pVTZ/SMD. Multiple conformations were generated for each tautomer, geometry optimization and energy calculations were performed with B3LYP/aug-cc-pVTZ in gas phase and solution (using the SMD continuum solvation model) and Δ_tG^calc_solv was calculated as the difference between the free energies in aqueous phase for the individual tautomers. Transfer free energy was calculated using B3LYP/aug-cc-pVTZ/SMD on the optimized geometries in solution. This approach will be abbreviated as B3LYP/aug-cc-pVTZ/B3LYP/aug-cc-pVTZ/SMD in the following. The RMSE between Δ_tG^exp_solv and Δ_tG^calc_solv for the dataset is 3.4 [3.0;3.7] kcal mol⁻¹ (the quantities [X;Y] denote a 95% confidence interval).

• B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD. Multiple conformations were generated for each tautomer, conformations were optimizing with B3LYP/aug-cc-pVTZ in gas phase and solution (using the SMD continuum solvation model). Transfer free energy was calculated using B3LYP/6-31G(d)/SMD on the optimized geometries in solution. Δ_tG^calc_solv was calculated as the difference between the free energy in aqueous phase for the individual tautomers and conformations. The RMSE between Δ_tG^exp_solv and Δ_tG^calc_solv for the QM dataset is 3.1 [2.7; 3.4] kcal mol⁻¹. This approach will be subsequently called B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD.

• B3LYP/aug-cc-pVTZ/SMD. Generating multiple conformations, optimizing with B3LYP/aug-cc-pVTZ in solution phase (using the SMD solvation model) and calculating relative solvation free energy Δ_tG^calc_solv directly as the difference between the free energy in aqueous phase (on the solution phase geometry). In this case, the free energy in aqueous phase is not obtained through a thermodynamic cycle, but the frequency calculation and thermochemistry calculations are performed with the continuum solvation model.³³ The individual free energy in aqueous phase G^calc_solv,k for conformation k are averaged to obtain the final G^calc_solv. The RMSE between Δ_tG^exp_solv and Δ_tG^calc_solv is 3.3 [3.0;3.7] kcal mol⁻¹. This approach will be subsequently called B3LYP/aug-cc-pVTZ/SMD.

The protocol used to obtain the results described above did not explicitly account for changes in internal rotors (only changes in the point group were considered). The results including changes in internal rotors are shown in Fig. S.I.5.† Including internal symmetry numbers did not improve the results.

The following discussion will concentrate on the results obtained with the best performing method (B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD). The results for all other methods are shown in the ESI Section in Fig. S.I.4.†

The tautomeric free energy difference in solution Δ_tG^calc_solv calculated with B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD are shown in Fig. 4. The RMSE for the calculated values is 3.1 [2.7;3.4] kcal mol⁻¹. Overall, the method tends to overestimate the tautomeric free energy difference in solution, as can be seen from Fig. 4B (lower panel). For 59 tautomer pairs (13% of the dataset) the method was not able to correctly determine the dominant tautomer species. For these 59 tautomer pairs the mean absolute error (MAE) between predicted and experimental value was 2.9 kcal mol⁻¹.

5 out of the 6 tautomer pairs with highest absolute error (above 10 kcal mol⁻¹) between experimental and calculated free energy difference have common scaffolds (values and molecules shown in the ESI Table S.I.1†). tp_1668, 1669, 1670 are all based on 5-iminopyrrolidin-2-one with either nitrogen, oxygen or sulfur on position 4 of the ring. tp_1559 and tp_853 both have methoxyethylpiperidine as a common substructure.

For the three analogues based on 5-iminopyrrolidin-2-one, we note that the experimental values are estimated results based on a three-way tautomeric reaction.³⁴ While this might point to the unreliability of Δ_tG^exp_solv, we can not draw any definitive conclusions here without repeating the underlying experiment.

The best performing quantum chemistry method shows very good retrospective performance on the SAMPL2 tautomer set

Some of the tautomer pairs deposited in the Tautobase dataset were also part of the SAMPL2 challenge, specifically 6 out of 8 tautomer pairs of the obscure dataset (using the notation of the SAMPL2 challenge: 1A_1B, 2A_2B, 3A_3B, 4A_4B, 5A_5B, 6A_6B) and 8 out of 12 tautomer pairs of the explanatory dataset (7A_7B, 10B_10C, 10D_10C, 12D_12C, 14D_14C, 15A_15B, 15A_15C, 15B_15C).¹² Comparing these molecules with the results of the SAMPL2 challenge—specifically with the results of the 4 participants with the best overall performance (ref. 16, 35, 15 and 14)—helps to assess the quality of the used approach. Arguably, the selection of the best performing method in the section above introduces a bias since the methods described subsequently have only seen the tautomer data described in the SAMPL2 challenge. This bias is difficult to quantify and correct, especially since some of the methods described below have been modified to reproduce the SAMPL2 explanatory tautomer set, but nevertheless we want to point it out to the reader at this point.

Since all four mentioned participants employed different methods, we will briefly describe the best performing methods of these publications for which results are shown in Table 1. The references to the methods used below are not cited explicitly, these can be found in the publications cited at the beginning of the following paragraphs.

Table 1 Comparison between the results of this work and the SAMPL2 challenge. All values are given in kcal mol⁻¹. The first part of the names in the ‘name’ column refers to the tautomer pair naming convention from the SAMPL2 challenge and in brackets is the name as used in this work. The calculations for Klamt and Diedenhofen¹⁶ were performed with MP2+vib-CT-BP-TZVP, for Ribeiro et al.¹⁵ with M06-2X/MG3S/M06-2X/6-31G(d)/SM8AD, Kast et al.¹⁴ with MP2/aug-cc-pVDZ/EC-RISM/PSE-3, and Soteras et al.³⁵ with MP2/CBS+[CCSD-MP2/6-31+G(d)](d)/IEF-MST/HF/6-31G). On this subset the presented approach was the second best method based on the total RMSE. Bracketed quantities [X,Y] denote 95% confidence intervals

Name	ΔtG^expsolv [kcal mol^−1]	ΔtG^calcsolv [kcal mol^−1]	ΔtG^calcsolv [kcal mol^−1]	ΔtG^calcsolv [kcal mol^−1]	ΔtG^calcsolv [kcal mol^−1]	ΔtG^calcsolv [kcal mol^−1]
		This work	16	15	14	35
1A_1B (tp_982)	−4.8	−4.7	−4.0	−3.0	−7.7	−4.6
2A_2B	−6.1	−6.8	−5.7	−5.7	−9.7	−6.3
3A_3B (tp_999)	−7.2	−8.4	−7.7	−6.7	−11.2	−7.7
4A_4B	−4.8	−0.4	0.5	0.8	−4.6	0.6
5A_5B (tp_1614)	−4.8	−4.7	−3.9	−4.4	−6.2	−5.6
6A_6B (tp_1005)	−9.3	−11.4	−7.6	−9.7	−11.2	−10.0
RMSE		1.3 [0.7,1.7]	1.4 [0.6,2.1]	1.5 [0.4,2.3]	2.8 [2.0,3.5]	1.3 [0.3,2.1]
7A_7B (tp_141)	6.6	4.9	5.3	6.5	5.1	5.5
10B_10C (tp_1058)	−2.9	−5.3	1.7	0.0	−2.8	2.2
10D_10C (tp_1059)	−1.2	−1.7	3.8	2.6	−0.6	5.0
12D_12C (tp_1514)	−1.8	−2.1	3.3	3.1	−0.8	3.0
14D_14C (tp_1072)	0.3	−1.6	1.9	0.8	0.2	4.0
15A_15B (tp_504)	0.9	6.1	−3.0	3.6	0.0	0.9
15A_15C (tp_505)	−1.2	0.7	−0.8	2.3	−1.9	1.4
15B_15C (tp_506)	−2.2	−5.0	1.8	−1.2	−1.9	0.5
RMSE		2.5 [1.5,3.6]	3.6 [2.5,4.5]	2.9 [1.8,3.8]	0.8 [0.5,1.1]	3.8 [2.5,4.9]
Total RMSE		2.2 [1.4,2.9]	2.9 [2.0,3.6]	2.4 [1.6,3.1]	1.9 [1.2,2.6]	3.0 [1.9,4.0]

---

Klamt and Diedenhofen¹⁶ used BP86/TZVP DFT geometry optimization in vacuum and the COSMO solvation model. Free energy in solution was calculated from COSMO-BP86/TZVP solvation energies and MP2/QZVPP gas-phase energies. Thermal corrections (including ZPE) were obtained using BP86/TZVP gas phase frequencies.

Ribeiro et al.¹⁵ calculated the free energy in solution as the sum of the gas-phase free energy and the transfer free energy. The gas phase free energy was calculated using M06-2X/MG3S level of theory and the molecular geometries optimized with the same method. The corresponding transfer free energy were computed at the M06-2X/6-31G(d) level of theory with the M06-2X/MG3S gas-phase geometries using the SM8, SM8AD, and SMD continuum solvation models (Table 1 shows only the results with SM8AD, which performed best).

Kast et al.¹⁴ optimized geometries in gas and solution phase (using the polarizable continuum solvation model PCM) using B3LYP/6-311++G(d,p). Energies were calculated with EC-RISM-MP2/aug-cc-pVDZ on the optimized geometries in the corresponding phase. The Lennard-Jones parameters of the general Amber force field (GAFF) were used.

Soteras et al.³⁵ used the IEF-MST solvation model parameterized for HF/6-31G(G) to obtain transfer free energy values. Gas phase free energy differences were obtained by MP2 basis set extrapolation using the aug-cc-pVTZ basis set at MP2/6-31+G(d) optimized geometries. Correlation effects were computed from the CCSD-MP2/6-31+G(d) energy difference.

The approach used in this work (B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD) performs well compared to the four approaches described above. For the total set of investigated tautomer pairs our approach has a RMSE of 2.2 [1.4,2.9] kcal mol⁻¹, making it the second best performing approach only outperformed by Kast et al.¹⁴

The difference in RMSE between the explanatory and blind dataset is noteworthy. Approaches that perform well on the blind data set perform worse on the explanatory set and vice versa. This is to a lesser extent also true for our chosen approach—B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD performs worse for the explanatory tautomer set (RMSE of 2.5 [1.5,3.6] kcal mol⁻¹) than on the blind tautomer set (RMSE of 1.3 [0.7,1.7]), but in comparison with the other four approaches, it is consistently the second best approach.

Interesting to note are the three tautomer pairs 15A_15B, 15A_15C and 15B_15C from the explanatory data set. The absolute error for these three pairs are 5.24, 1.9, and 2.8 kcal mol⁻¹, respectively. It appears that the used approach has difficulty to model 15B correctly, showing larger than average absolute errors whenever 15B is part of the tautomer reaction. Most likely the hydroxyl group in 15B is critically positioned and sensitive to partial solvent shielding by the phenyl ring, something that has been noted before.¹⁴

The discrepancy between the different approaches (ours included) for the tautomer set shows that it seems highly difficult to propose a single method for different tautomer pairs that performs consistently with a RMSE below 2.0 kcal mol⁻¹. This issue is made substantially worse by the many different ways methods can be used/combined and errors can be propagated/compensated during tautomeric free energy difference Δ_tG^calc_solv calculations. While we believe that a RMSE of 3.1 kcal mol⁻¹ is a good value for the chosen approach, especially when compared to the results of the SAMPL2 challenge, it is by far not a satisfying result. The accuracy, compared to the cost of the approach, is not justifiable and there is still a dire need for more accurate and cheaper methods to obtain relative solvation free energies for tautomer pairs.

With this we close the investigation of the state of the art QM free energy calculations and will investigate some of the shortcomings of the RRHO thermochemistry calculations. Some of the insight gained above could have been concluded from the work of Kast et al.,¹⁴ Ribeiro et al.,¹⁵ Klamt and Diedenhofen,¹⁶ Soteras et al.,³⁵ but the small number of chemical species investigated there and the often contradicting results reported elsewhere (e.g. ref. 36 and 37) made a closer investigation of QM free energy calculations for a large dataset necessary. The reported results will be of importance for further method development for tautomeric free energies.

Including multiple minimum conformations seems to have little effect on the accuracy of tautomer free energies

The three approaches described above consider multiple conformations to obtain the tautomeric free energy difference Δ_tG^calc_solv. Obtaining the global minimum conformation is an important and well established task in quantum chemistry calculations. Fig. S.I.6† shows the difference between the highest and lowest free energy for an optimized minimum conformation for the individual tautomers of the dataset generated with B3LYP/aug-cc-pVTZ. While many of the tautomers have only a single minimum conformation (278 out of 936 tautomers), for molecules with multiple minimum conformations the difference in the free energy emphasizes the need and justifies the cost for a global minimum conformation search. Molecules with more than 10 kcal mol⁻¹ difference between highest and lowest energy minimum are shown in Fig. S.I.6.†

While the scientific community agrees on the importance of the global minimum for property calculations, the importance of considering multiple conformations for quantum chemistry free energy calculations has not been well established (e.g. ref. 38). Using multiple minimum energy conformations can add substantial computational cost to the free energy calculations. Often, the global minimum conformation search can be performed with a lower level of theory than the single point energy calculation, making only a single high level electronic energy calculation necessary. Frequency calculations can also add considerably to the computational cost of the free energy calculation. If a single, minimum energy conformation is sufficient to obtain a good estimate for the free energy in gas phase substantial amount of computational time could be saved.

Fig. 6 compares the tautomeric free energy difference obtained using the weighted average over multiple minimum conformations and the single global minimum conformation. The results indicate that using multiple minimum conformations does not significantly improve the final result compared to using a single, global minimum structure. Only 12 tautomer pairs show an absolute deviation larger than 1 kcal mol⁻¹ in the tautomeric free energy difference Δ_tG^calc_solv between the two approaches. Much more relevant than including multiple conformations is locating the global minimum conformation (as clearly shown by Fig. 5).

Fig. 5 Computed free energies are highly sensitive to the selected minimum conformation. For each molecule, the number of minimum conformations N_conf is plotted against the difference between the corresponding highest and lowest obtained free energy value for the minimum conformations. 278 out of 936 molecules have only a single minimum; molecules with multiple minima show substantial free energy differences between the minimum conformations, highlighting the need for a global minimum search.

Fig. 6 A single, global minimum conformation can be used to calculate tautomeric free energy differences in solution without loss of accuracy. These results are based on the calculations with B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD. (A) shows the tautomeric free energy difference Δ_tG^calc_solv obtained with multiple minima (m.m.) plotted against the global, single minimum (s.m.) tautomeric free energy difference Δ_tG^calc_solv. (B) shows the KDE and histogram of the difference between s.m. Δ_tG^calc_solv and m.m. Δ_tG^calc_solv. The comparison indicates that there is little benefit using multiple minimum structures over a single, global minimum considering the high costs of the former. Comparing two estimated properties we use root mean squared deviation (RMSD) and mean absolute deviation (MAD) instead of error.

Alchemical relative free energy calculations with quantum machine learning potentials can rigorously capture classical statistical mechanical effects

Previous work using QML potential for free energy calculations have focused on hybrid QML/MM simulations wherein ligand interactions are treated with ML and the environment and ligand–environment interactions with MM.^39,40 To our knowledge, this is the first time alchemical free energy calculations have been performed using purely a QML potential for drug-like molecules.

Alchemical relative free energy calculations were performed for 354 tautomer pairs using 11 alchemical λ states in vacuum. In the following, we will compare the tautomeric free energy difference obtained using alchemical relative free energy calculations to the multiple minima, RRHO approximation using the same potential energy function (ANI-1ccx) to assess potential errors in the thermochemistry corrections. We will also show how a small number of experimentally obtained tautomer ratios in solution can be used to incorporate crucial solvent effects and recover tautomer free energies in solution by QML parameter optimization and importance weighting.

While ANI-1x was trained on energies/forces calculated with ωB97x/6-31G*(d), ANI-1ccx was retrained from the same level of theory to include energies recalculated with coupled cluster considering single, double, and perturbative triple excitations (CCSD(T)) with an extrapolation to the complete basis set limit (CBS).⁴¹ CCSD(T)/CBS is the gold standard for electronic energy prediction and calculates thermochemical properties to within the limit of chemical accuracy (e.g. ref. 42).

RRHO ANI-1ccx calculations show significant deviations from the alchemical relative free energy calculations

In the limit of infinite sampling alchemical relative free energy calculations approach the exact free energy difference. Alchemical relative free energy calculations can be used to quantify the error introduced by a discrete partition function in the form of multiple minimum conformations and harmonic treatment of all bonded terms (including torsions and internal rotors)—if the same potential energy function is used for both calculations. Such a comparison is especially useful for quantum mechanics potentials which typically do not include explicitly defined harmonic terms in their functional form (as is the case with most classical force fields). It is important to emphasize that what is quantified here is not the error that the RRHO introduces for free energy calculations, but for relative free energy calculations: the difference in the partition function that results from the tautomerization. Typically, only a small region of the molecule is affected by the tautomerization and there is potential for error compensation.

Since the simulation time of the individual lambda states for the alchemical relative free energy calculations were relatively short (200 ps) we repeated the calculations multiple times (5) with randomly seeded starting conformations and velocities to detect systems for which the simulation time was clearly insufficient. In the following we will only use systems that had a standard deviation of less than 0.3 kcal mol⁻¹ for 5 independent alchemical free energy calculations. Applying this filter resulted in the removal of 65 tautomer pairs for which the free energy calculation had not converged.

Results shown in Fig. 7 indicate the average deviation that every relative free energy calculation based on the RRHO approximation introduces, regardless of the accuracy of the actual potential to model electronic energies. The mean absolute deviation of 0.9 kcal mol⁻¹ should not be underestimated. Results shown in Table 1 would be significantly improved if an error of 0.9 kcal mol⁻¹ could be compensated by running a protocol that samples the relevant conformational degrees of freedom to obtain an exact partition function.

Fig. 7 Independent of the level of theory, thermochemistry corrections introduce a mean absolute deviation (MAD) of 0.9 kcal mol⁻¹. ANI-1ccx is used for the calculation of tautomeric free energy differences in vacuum Δ_tG^calc_vac using alchemical relative free energy (AFE) calculations and single point energy calculations on multiple minimum conformations and thermochemistry corrections based on the RRHO approximation. (A) shows a scatter plot between the two approaches and (B) the KDE and histogram of the difference between the two approaches. Error bars are shown in red on the alchemical relative free energy estimates, these were obtained from the MBAR estimate of the relative free energy (see Methods section). Comparing two estimated properties we use root mean squared deviation (RMSD) and mean absolute deviation (MAD) instead of RMSE/MAE.

For 12 (out of 289) tautomer pairs the multiple minima RRHO approximation deviated by more than 3 kcal mol⁻¹ (molecules are shown in Fig. S.I.1†). Most of these molecules have high conformational degrees of freedom and it seems unlikely that a naive enumeration of relevant conformations (e.g. with a conformer generator) will detect all of them—this might contribute to the observed error. That is certainly true for tp_113, 116, 565, 554, 403, 1674.

QML potentials can be optimized to reproduce experimental tautomer ratios in solution

Since the alchemical free energy calculations were performed in vacuum, a comparison with the experimental tautomeric free energy difference Δ_tG^exp_solv showed a high RMSE of 6.4 [5.9,6.8] kcal mol⁻¹, as expected from the well-known impact of solvation effects on tautomer ratios.³⁷

In the following we want to investigate if thermodynamic observable—in this specific case tautomer ratios—can be used to retrain a neural net potential derived from QM calculations. Further, we want to test how such an optimized parameter set would perform on the original dataset used to train the neural net potential and if it is possible to use the difference between the energies of the original and optimized parameter set for regularization.

Defining a loss function L as the error between the calculated Δ_tG(θ)^calc_vac and experimental Δ_tG^exp_solv free energies it is possible to optimize the loss with respect to the neural net parameters θ defining the ANI potential. Using importance weighting, a new free energy estimate can be calculated with the optimized QML parameters (θ*) using the modified potential energy function without resampling the equilibrium distributions. To ensure that the quality of the neural net parameters does not deteriorate beyond a reasonable threshold a regularization term was added that acts on the energy difference for each snapshot calculated with the endstate potential energy functions with the original and optimized parameter set (ΔE(θ, θ*)). The regularization term was included in the molecular loss function acting on the snapshots used for the free energy calculation for the 212 tautomer pairs in the training set—in the course of a single epoch ΔE(θ, θ*) is evaluated on a total of 699 600 snapshots (212 tautomer pairs × 11 lambda states × 300 conformations).

Initial results led to the introduction of scaling factors that allows to slowly increase the contribution of the tautomeric free energy difference in the loss function (and/or the regularization term) during the training. The protocol is described in more details in the Methods section.

It was possible to overfit the parameters on the training set within 50 epochs if no regularization term was used (training/validation performance shown in Fig. S.I.7I and M†). Without a regularization term the RMSE for ΔE(θ, θ*) on the training/validation/test snapshots rises to 40–100 kcal mol⁻¹. While the training set performance improves the validation set performance reaches a plateau around 2.2–2.6 kcal mol⁻¹. The red dotted line in Fig. S.I.7† top panels indicates the performance of the reported results in Fig. 8, indicating that even without regularization this performance can not be improved significantly.

Fig. 8 Optimizing QML parameters on a set of experimentally obtained tautomer free energies in solution Δ_tG^exp_solv enables ANI-1ccx to include crucial solvation effects and improved estimates for tautomeric free energy differences can be obtained by importance weighting from vacuum simulations using the optimized QML parameters. (A) Top panel shows the training (green) and validation (purple) set performance as ΔΔ_tG_solv. Validation set performance was plotted with a bootstrapped 95% confidence interval. The performance of the optimized parameter set is also shown on the original ANI-1ccx dataset in blue. The best performing parameter set (evaluated on the validation set and indicated by the red dotted line) was selected to evaluate its performance on the test set. The bottom panel shows the MAE for the energy difference between each of the 400 parameter sets and the original parameter set on all the snapshots used for the free energy calculations (≈1,2 million snapshots) split in training/validation and test set as well as the original ANI1-ccx dataset. Figure (B) shows the distribution of Δ_tG^exp_solv − Δ_tG^calc_solv for a hold out test set (71 tautomer pairs) with the native ANI-1ccx (θ) and the optimization parameter set (θ*). The optimized parameter set was able to improve the prediction of tautomeric free energy differences from initial 6.7 kcal mol⁻¹ to 2.8 kcal mol⁻¹ (MAE improved from 5.3 to 2.0 kcal mol⁻¹). The difference in Kullback–Leibler divergence (KL) indicates that the tautomeric free energy differences obtained with the optimized parameter set can reproduce the distribution of the experimental tautomer ratios much better than the free energy differences obtained with the original parameter set.

The results shown in Fig. 8 and S.I.7† indicate that there is a trade off between the accuracy of the ANI-1ccx parameters on its original training set (either directly shown in Fig. S.I.7† or indirectly through ΔE(θ, θ*) in Fig. 8) and the improvement for the prediction of tautomeric free energy differences. This is best exemplified in the training/validation set performance shown in Fig. S.I.7F and X.† For both training runs the scaling factor of the regularization term was increased throughout the epochs, shifting the focus from optimizing the free energy to keeping ΔE(θ, θ*) as small as possible. In Fig. S.I.7X† the ΔE(θ, θ*) term in the loss was kept constant throughout the first 100 epochs and then raised approximately tenfold during the next 50 epochs. Looking at the training/validation set performance it becomes evident that starting with epoch 100 the initial trend towards improved performance starts to reverse while the MAE for ΔE(θ, θ*) decreases. After 400 epochs the MAE for ΔE(θ, θ*) is near zero and there is no improvement in the free energy estimate for the training/validation set. This indicates that the parameter set approaches the same minimum that was occupied before training started.

Whenever both weights and biases are trained for the QML potential it often takes a few hundred epochs before ΔE(θ, θ*) approaches reasonable values (∼2 kcal mol⁻¹), shown e.g. in Fig. S.I.7L, J, G, and A.† Similar behavior was observed even if the learning rate for the biases was reduced to 1 × 10⁻⁶ as shown in Fig. S.I.7T.† In some instances it was not possible to converge and successfully minimize both terms of the loss function in 400 epochs. One such case is shown in Fig. S.I.7E.†

In Fig. 8A the training/validation set performance of a training run with 400 epochs is shown. Only weights were optimized with a learning rate of 1 × 10⁻⁴, regularization was used (the scaling factors are plotted in Fig. S.I.9†). Fig. 8B shows the obtained free energy estimates with the best performing parameter set (calculated on the validation set) on an independent test set (71 tautomer pairs). The optimized parameter set (θ*) was able to improve the prediction of tautomeric free energies on the test set from initial 6.7 [5.7,7.7] kcal mol⁻¹ with the original parameter set (θ) to 2.8 [2.2,3.3] kcal mol⁻¹, comparable to the performance of B3LYP/aug-cc-pVTZ/B3LYP/6-31G(d)/SMD. Fig. 8B shows the distribution of . The Kullback–Leibler divergence (KL) value for indicates that the optimized parameter set is able to reproduce the distribution of the experimental free energies better than the initial parameter set.

Fig. 8A shows the performance of the parameter set on the ANI-1ccx training dataset.§ The RMSE of the parameter set for the ∼500 000 data points of the dataset increases from initial 1.7 kcal mol⁻¹ to 5.5 kcal mol⁻¹. This increase in RMSE can be partially explained by the target of the retraining: tautomeric free energies in solution. Including solvation effects necessarily decreases the performance on a gas phase dataset.

There is a limit to the reweighting workflow, as indicated in Fig. S.I.8,† that requires to re-sample the equilibrium distributions with the perturbed parameter set. The uncertainty of the perturbed free energy estimate increases as the weights become more concentrated on fewer observations and the effective sample size shrinks.⁴³ Using an arbitrary cutoff of 1 kT as acceptable for the perturbed free energy uncertainty Fig. S.I.8† indicates that conformations should be resampled after around 400 epochs with the perturbed parameter set.

Coming back to the question posed in the beginning of this section: yes, it appears to be possible to use experimental values of a thermodynamic observable to retrain an already optimised QML potential to improve its performance for predicting free energies. This is astonishing, considering the small set of experimental measurements used to retrain the potential (≈250 data points) that provides an unambiguous and substantial improvement in performance on the held-out test set. The learning curve shown in Fig. S.I.12† indicates that using 5% of the training data already improves the MAE by 2 kcal mol⁻¹.

The regularization term restraining the QML parameter set so that it is still able to reproduce single point energies with comparable quality than the original parameter set (with an increase in RMSE of ∼3 kcal mol⁻¹) offers an interesting trade off between high quality QM single point energies and highly improved tautomeric free energy differences. It appears as if the regularization term is not necessary to improve the tautomeric free energy difference estimate, but introducing such a regularization term in the loss function does also not significantly hinder the parameter set optimization (if the scaling factors for the two terms in the loss function are chosen reasonably). These results highlight the incredible potential of QML energy functions that are differentiable with respect to its parameters in addition to coordinates, enabling model tuning based on experimental and quantum chemical data to be both facile and incredibly powerful.

Discussion & conclusion

In this work we use a state of the art density functional theory protocol and continuum solvation model to calculate tautomer ratios for 460 tautomer pairs using three different approaches to model the solvent contributions. The best performing method uses B3LYP/aug-cc-pVTZ and the RRHO approximation for the gas phase free energies (calculated on B3LYP/aug-cc-pVTZ optimized geometries). The transfer free energy was calculated using B3LYP/6-31G(d)/SMD on geometries optimized in their respective phase (with B3LYP/aug-cc-pVTZ). This approach performs with an RMSE of 3.1 kcal mol⁻¹.

One possible source of error—independent of the method used to calculate the electronic energy and model the continuum electrostatics—are the thermochemical corrections used to obtain the standard state free energy. Typically, an analytic expression is used to approximate the partition function which is based on the rigid rotor harmonic oscillator approximation. To obtain the correct and unbiased partition function and compute rigorous free energy estimates we implemented an alchemical relative free energy workflow using the ANI family of quantum machine learning (QML) potentials.²³ The method was implemented as a python package and is available here: https://github.com/choderalab/neutromeratio.

Using the same potential for the calculations based on the RRHO and performing alchemical relative free energy calculations we are able to show that the RRHO approximation introduces a mean absolute error of ∼1 kcal mol⁻¹ in the calculation of the investigated relative free energies. These errors can be attributed to anharmonicity in bonded terms, difficulties to enumerate relevant minimum conformations and in combining shallow local energy wells as well as the inconsistent treatment of internal and external symmetry numbers.

The ANI family of QML and comparable QML potentials have opened the possibility to investigate tautomer ratios using relative free energy calculations without prohibitive expensive MM/QM schemes or ab initio simulations. The calculated alchemical free energies obtained using the methods implemented in the “Neutromeratio” package can be used to optimized the parameters of the QML potential. Using a small set of experimentally obtained tautomer ratios we were able to optimize the QML parameters on a training set and significantly improve the accuracy of the calculated free energies on an independent test set. The tautomeric free energies obtained with the optimized parameter set improved the RMSE from initial 6.7 kcal mol⁻¹ to 2.8 kcal mol⁻¹ on a hold-out test set of 71 tautomer pairs.

What should be noted here: the experimental values are relative solvation free energies, while we calculate relative gas phase free energies. To calculate correct tautomeric free energies in solution, solvent effects need to be included. In this work solvent effects are modeled by fitting the ANI QML parameters to experimental relative solvation free energies. The use of explicit solvent molecules is the preferable and recommended solution. The optimization of QML parameters on a small set of experimental free energies can be extended easily for explicit solvent simulations (in Table S.I.2† we show results for alchemical free energy calculations for 6 of the previously investigated tautomer systems in a water droplet).

Obtaining accurate free energy differences between tautomer pairs in solvent remains an elusive task. The subtle changes and typically small difference in internal energies between tautomer pairs require an accurate description of electronic structures. Furthermore, solvent effects have a substantial effect on tautomer ratios; consequently, a proper descriptor of solvation is essential. The change in double bond pattern typically also induce a change in the conformational degrees of freedom and, related, in the conformation and rotational entropy. But—despite all of these challenges—we remain optimistic that further developments in fast and accurate neural net potentials will enable improved and more robust protocols to use relative free energy calculations to address these issues.

Detailed methods

Experimental data

The full dataset considered for this study was obtained from the DataWarrior File deposited in https://github.com/WahlOya/Tautobase (commit of Jul 23, 2019), described in detail in ref. 9. The dataset was sourced from the tautomer codex authored by P. W. Kenny.³⁴

From the dataset a subset of tautomer pairs were considered that

(1) were measured/calculated/estimated in aqueous solution

(2) had a numeric log K value between ±10

(3) had no charged species

(4) did not contain iodine

(5) only a single hydrogen and double bond change its position.

476 of the 1680 deposited tautomer pairs had these properties. We added two tautomer pairs from the SAMPL2 challenge (tautomer pair 2A_2B and 4A_4B).¹² 478 unique tautomer pairs were considered for further analysis. The term ‘unique’ tautomer pair was defined in the following as containing a unique combination of two molecules. We only consider tautomer pairs, not tautomer equilibria with multiple tautomer forms. In the following we will use an identifier containing the row number entry from the DataWarrier file to identify tautomer pairs in the dataset, e.g. tp_200 describing the tautomer pair (tp) at row number 200 in the original DataWarrior file.

The log K value was converted to free energies with Δ_tG^exp_solv = −RT ln K. The free energy difference of the tautomer pairs obtained from the Tautobase are subsequently referenced as Δ_tG^exp_solv in contrast to the calculated values which are called Δ_tG^calc. While we call all values deposited in the Tautobase Δ_tG^exp_solv we want to point out that some of these values are estimated or calculated.

Unfortunately, experimental error estimates are not deposited in the Tautobase and subsequently also not modeled in this manuscript. This highlights the need for further experimental measurements of tautomer ratios.

A closer inspection of some of the outliers from preliminary calculations identified molecules with incorrect structures in the database (row entry: 1260, 1261, 1262, 1263, 1264, 514, 515, 516, 517, 1587)—these 10 tautomer pairs were subsequently removed from all QM calculations (simulations using ANI included the corrected tp_1260, tp_1261 and tp_1587). These structures are corrected in the current version of the Tautobase. Additionally, 8 tautomer pairs containing bromide (row entry: 989, 581, 582, 617, 618, 83, 952, 988) were removed from calculations performed with the basis set 6-31G(d) due to the lack of adequate parameters.

ANI-1 is parameterized on molecules containing only carbon, nitrogen, hydrogen and oxygen, making a further removal of tautomer pairs containing elements beside the aforementioned necessary.²³ This resulted in 369 tautomer pairs.

The tautomer pair set used for relative alchemical free energy calculations needed one additional filter. Molecules with a stereobond that changes its position between tautomers had to be removed (this affected 15 tautomer pairs: tp_1637, tp_510, tp_513, tp_515, tp_517, tp_518, tp_787, tp_788, tp_789, tp_810, tp_811, tp_812, tp_865, tp_866, tp_867). Such stereobonds introduce additional complexity since it would be necessary to introduce restraints to define the stereochemistry during the lambda protocol.

The subset of the Tautobase used for the QM calculations can be obtained here as list of SMILES (468 tautomer pairs): https://github.com/choderalab/neutromeratio/blob/master/data/b3lyp_tautobase_subset.txt. The subset of the Tautobase used for the QML calculations can be found here as list of SMILES (354 tautomer pairs): https://github.com/choderalab/neutromeratio/blob/master/data/ani_tautobase_subset.txt. The distribution of the experimental tautomer free energies for both datasets is shown in Fig. 9.

Fig. 9 The tautomer dataset shows a wide variety of solvation free energies ΔG^exp_solv. (A) shows the ANI-Tautobase subset that was used for the QML calculations and (B) shows the QM-Tautobase subset used for the QM calculations. The full dataset considered for this study was obtained from the DataWorrier File deposited at https://github.com/WahlOya/Tautobase (commit of Jul 23, 2019), described in detail in ref. 9. The selection criteria for both datasets are described in detail in Detailed methods section.

Generating molecular conformations

The input tautomer pairs were specified as SMILES strings. 3D conformations were generated with the chemoinformatics toolkit RDKit (version 2019.09.2) which uses the distance geometry approach to generate conformations while enforcing chirality/stereochemistry^44,45 For each molecule 20 conformations were initially generated. The number of conformations was reduced to 10 if the average root mean square deviation (RMSD) of atom coordinates of the 20 conformations was below 0.5 Å, and further reduced to 5 conformations if below 0.2 Å.

RMSD calculations and filtering of conformations

For each molecule, pairwise RMSD between conformations were calculated using RDKit. Starting with a random conformation, if the RMSD to any other conformation of the molecule was below 0.1 the conformation is discarded, otherwise added to the list of unique conformations. The RMSD was calculated between heavy atoms and the hydrogen of selected chemical moieties including primary alcohols, imines, primary/secondary amines, cyanamides and thiols.

Combining energies of conformations

Energies of different minimum conformations of the same molecule were weighted and combined using the following schemewith N as the number of minimum conformations and G the free energy at the minimum conformation k. The obtained free energies will be called weighted free energies in the following.

Quantum mechanical calculations

The following quantum mechanical (QM) calculations were performed using the quantum chemical software orca 4.0.1.2.⁴⁶ The universal solvation model based on solute electron density (SMD) was used as continuum solvation model.³²

Geometric optimization, single point energy and frequency calculation. Geometric optimization was performed with the standard options of orca, redundant internal coordinates and the BFGS optimizer.^27–29 Frequency calculations were performed with the numerical Hessian computed using the central differences approach. If a conformation had negative frequencies (imaginary modes) after the geometry optimization it was excluded from further analysis. Single point calculations were performed on the optimized geometries using B3LYP and the basis set aug-cc-pVTZ or 6-31G(d).^30,31 A damping dispersion correction was applied (orca keyword D3BJ).⁴⁷

Continuum solvation model. The continuum solvation model SMD was used to model the molecules in aqueous environment.³² Since there is a volume change in the standard state from the gas phase (1 atm) to the solvent phase the gas phase standard state is indicated by ‘*’ and the solvation standard state by ‘°’. The standard-state transfer free energy is then defined aswith ΔG^*→° as standard-state adjustments (specifically, the correction of changing the volume from the gas phase to the solute phase with a constant value of 1.89 kcal mol⁻¹), ΔG_ENP describes the electronic (E), nuclear (N) and polarization (P) components of the free energy and ΔG_CDS free energy changes associated with solvent cavitation (C), changes in dispersion (D) and changes in local solvent structure (S).³²

Free energy and free energy in solution calculations. Thermal corrections were computed at standard state (298.15 K and 1 atm pressure) using the ideal gas molecular partition function and the rigid-rotor harmonic oscillator (RRHO) approximation. Low-lying vibrational frequencies (below 15 cm⁻¹) were treated by a free-rotor approximation⁴⁸—this method is also sometimes called rigid-rotor quasi harmonic oscillator. The external (rotational) symmetry number was obtained from the point group of the tautomer using the point group module of Jmol and visual inspection and used to correct the rotational entropy calculated by orca.⁴⁹

The Gibbs free energy in gas phase for a given coordinate set (k) was obtained by adding thermal corrections for conformer k at temperature T and zero point energy (ZPE) contributions ε_ZPE,K to the electronic energy E_k.²² The degeneracy D describes the entropy contribution of internal rotors and is not added for the calculations shown in the main text (results including degeneracy D are shown in the ESI†).

The degeneracy D was estimated by calculating the graph automorphism of the molecule. The implementation of the VF2 algorithm for graph isomorphism of networkx was used.⁵⁰ Nodes were defined to match if element and hybridization matched, edges were identical if bond order matched.

The tautomeric free energy difference in solution Δ_tG^calc_solv can be calculated from the standard-state Gibbs free energy in aqueous phase of the product and educt of the corresponding tautomer reaction, which itself is calculated as the sum of the gas-phase standard-state free energy and the standard-state transfer free energy expressed as

for a given conformation k and shown as a thermodynamic cycle in Fig. 1.

An alternative way to calculate the free energy in aqueous phase used in the manuscript is

where ψ^sol is the polarized wave function in solution, H^g the gas phase Hamiltonian, and V the potential energy operator associated with the reaction field. The bracket term describes the electronic energy, while G_NES is associated with non-electrostatic contributions (dispersion–repulsion and solvent structural terms) to the solvation energy and G_T,K are the thermal correction calculated directly in the continuum solvation model.⁵¹

ASE thermochemistry corrections

ANI-1ccx was used to calculate the electronic energy. The conformations were minimized using the BFGS optimizer as implemented in scipy.⁵² Frequency and thermochemistry calculations were performed using the optimized geometry. Thermal corrections were calculated at 300 K and 1 atm using the IG-RRHO approximation as implemented in the atomic simulation environment (ase).⁵³ The tautomeric free energy difference in gas phase Δ_tG_vac was then calculated as the difference between of the tautomer pair.

Relative alchemical free energy calculations

Relative alchemical free energies were calculated using a single topology approach. The topology of each tautomer pair only differed in the position of a single hydrogen (atom types and bonds were not specified). The hybrid topology (the superset of the two topologies) differed therefore by one hydrogen from each of the physical endstates.

By default, the coordinates of the hybrid topology were generated by using the coordinates of ‘Tautomer1’ (as defined in the Tautobase database). If a tautomer isomerism created or removed a cis/trans stereobond, the initial coordinates were taken from the topology with the stereobond present (therefore sometimes changing the direction of the tautomer reaction).

The coordinates of the added, non-interacting hydrogen were obtained by randomly sampling 100 positions on the surface of a sphere (with radius of 1.02 Å) defined around its new bonded heavy atom and subsequently using the lowest energy position as the starting conformation for each lambda window for the free energy calculation. The physical endstates (representing the two tautomer states, each with an additional non-interacting (dummy) hydrogen) were connected via 11 equidistant intermediate (lambda) states.

Energy and forces were calculated using ANI1-ccx as implemented in torchani https://github.com/aiqm/torchani). The energy and force was linearly scaled along the alchemical path as a function of lambda with

E = (1 − λ)E₁ + λE₂ (7)

with E₁ and E₂ representing the potential energy at the physcial endstates. Before each simulation, initial coordinates were minimized using a BFGS optimizer as implemented in scipy.⁵² Coordinates were sampled using Langevin dynamics at 300 K with a collision rate of 10 ps⁻¹ and a 0.5 fs time step using the BAOAB integrator.⁵⁴ Initial velocities were obtained from a Maxwell–Boltzmann distribution at the simulation temperature.

Because nuclei can rearrange to form distinct chemical species in highly perturbed simulations, we applied a flat-bottom harmonic restraint to covalent bonds to ensure we sampled the desired chemical species in initial and final states of the free energy calculation. The restraint was defined as

with H as the Heaviside step function, Δ(r_i,j,r⁰_i,j) as the difference between the reference bond length r⁰_i,j and the current bond length r_i,j and r_fb as half of the well radius. For all heavy atom pairs r⁰_i,j was set to 1.3 Å and r_fb to 0.3 Å with k_i,j set to k_BT/0.1 Å. For C–H/O–H/N–H bond pairs r⁰_i,j was set to 1.02 Å (the average of the three different equilibrium bond length of C–H/O–H/N–H bond pairs) and r_fb to 0.4 Å with k_i,j set to k_BT/0.2 Å. The restraint well was chosen so that the restraint does not interfere with normal bond stretching but will activate once a bond is broken.

Samples were obtained from 200 ps simulations for each lambda state. For each tautomer pair calculations were repeated 5 times with randomly seeded initial velocities (and coordinates). Relative alchemical free energies ΔG^calc_vac were calculated using MBAR as implemented in the pymbar package.⁴³ 300 uncorrelated snapshots were considered for the MBAR analysis from each lambda state.

Neural net parameter optimization based on experimental relative solvation free energies

The tautomer data set was randomly split (20 : 20 : 60) into a test set (71), validation set (71) and training set (212 tautomer pairs). Neural net parameters were optimized using a routine modified from the TorchANI tutorial.¶ To limit the capacity to overfit, only the weights and biases of the final layer of each of the 8 pretrained ANI-1ccx models were optimized for each of the atom nets (one net per element), resulting in roughly 8 × 4 × 97 tunable parameters (8 neural nets, each with 4 atom nets containing 97 weights and bias) for ANI-1ccx. As in the TorchANI tutorial, the weight matrices were updated using the Adam optimizer with decoupled weight decay (AdamW), and the bias vectors were updated using Stochastic Gradient Descent (SGD). The training data was randomly partitioned in each epoch in mini-batches of 10 tautomer pairs and gradient updates were performed for each mini-batch. Training was performed for 400 epochs. The best model was chosen based on the RMSE on the validation set and model performance reported on the test set.

The model was trained by minimizing the mean squared error (MSE) loss between calculated and experimental relative free energies. If regularization was used, the mean absolute error (MAE) between the energy for the parameter set θ* and the original parameter set θ was calculated on the individual snapshots used for the free energy calculation (11 × 300 snapshots) with the potential energy function at λ = 0 and λ = 1. The regularization term was normalized using the number of atoms of the tautomer system. The per molecule pair (m) loss function l is defined as

with ΔG^exp_solv,m as the experimental and as the calculated tautomeric free energy differences for tautomer pair m using parameters θ*. The two scaling factors f(epoch) and g(epoch) were used to control the contribution of the two terms as a function of training epoch. For the final results the values for g (labeled ‘scaling dG’) and f (labeled ‘scaling dE’) is shown in Fig. S.I.9.†

The overall loss is then

The perturbed relative free energy ΔG(θ*)^calc_vac,m required for computing l(m; θ*) is calculated by importance sampling, or “reweighting” the original MBAR estimate from the original pre-trained parameters θ to the current parameters θ*. To compute this estimate efficiently at arbitrary θ*, we first collect configuration samples at a reference value θ (corresponding to the original parameters of the pretrained ANI-1ccx model) for each intermediate value λ. For each configuration sample x, we compute the reduced potential u(x, λ; θ), to form the N × M matrix of inputs to MBAR (where N is the total number of snapshots, M is number of λ windows). MBAR equations are solved to yield a vector of reduced free energies , where m indexes the intermediate λ values. The relative free energy prediction implied by the model parameters θ* is then .

To compute relative free energies at a new value of the parameters θ*, we need to compute u(x, λ = 1; θ*) and u(x, λ = 0; θ*) for all configurations x. The optimization routine evaluates the gradient of the loss function L w.r.t. θ* by automatic differentiation and updates the parameters.

300 uncorrelated snapshots per λ state were used for the MBAR estimate and importance sampling for the vacuum simulations.

Disclosures

JDC is a current member of the Scientific Advisory Board of OpenEye Scientific Software, Redesign Science, and Interline Therapeutics, and holds equity interests in Redesign Science and Interline Therapeutics. The Chodera laboratory receives or has received funding from multiple sources, including the National Institutes of Health, the National Science Foundation, the Parker Institute for Cancer Immunotherapy, Relay Therapeutics, Entasis Therapeutics, Silicon Therapeutics, EMD Serono (Merck KGaA), AstraZeneca, Vir Biotechnology, Bayer, XtalPi, Foresite Laboratories, the Molecular Sciences Software Institute, the Starr Cancer Consortium, the Open Force Field Consortium, Cycle for Survival, a Louis V. Gerstner Young Investigator Award, and the Sloan Kettering Institute. A complete funding history for the Chodera lab can be found at http://choderalab.org/funding.

Funding

JF acknowledges support from NSF CHE-1738979 and the Sloan Kettering Institute. MW acknowledges support from a FWF Erwin Schrödinger Postdoctoral Fellowship J 4245-N28. JDC acknowledges support from NIH grant P30 CA008748, NIH grant R01 GM121505, NIH grant R01 GM132386, and the Sloan Kettering Institute.

Data availability

• Python package used in this work (release v0.2): https://github.com/choderalab/neutromeratio

• Data and notebooks to reproduce the plots/figures (release v0.2): https://zenodo.org/record/4562426

Author contributions

Conceptualization: JDC, JF, and MW; methodology: JDC, JF, and MW; software: JF, and MW; investigation: JF, and MW; writing – original draft: JDC, JF, and MW; writing – review & editing: JDC, JF, and MW; funding acquisition: JDC and MW; resources: JDC and MW; supervision: JDC and MW.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

MW, JF and JDC are grateful for discussions from the Tautomer Consortium, specifically Paul Czodrowski, Brian Radak, Woody Sherman, David Mobley, Christopher Bayly and Stefan Kast as well as input from Adrian Roitberg and Olexandr Isayev. This work was only possible through the time and effort Thomas Sander and Oya Wahl invested in curating the open tautomer database (Tautobase). Further, MW is grateful for the support of Thierry Langer, who donated significant computational resources to perform much of the calculations/simulations and Oliver Wieder, for helpful discussions.

References

1. I. Kapetanović, Drug Discovery and Development: Present and Future, IntechOpen, 2011.
2. B. Bax, C. W. Chung and C. Edge, Getting the chemistry right: protonation, tautomers and the importance of H atoms in biological chemistry, Acta Crystallogr., Sect. D: Struct. Biol., 2017, 73(2), 131–140.
3. M. Sitzmann, W. D. Ihlenfeldt and M. C. Nicklaus, Tautomerism in large databases, J. Comput.-Aided Mol. Des., 2010, 24(6–7), 521–551.
4. A. R. Katritzky, C. Dennis Hall, B. E. D. M. El-Gendy and B. Draghici, Tautomerism in drug discovery, J. Comput.-Aided Mol. Des., 2010, 24(6–7), 475–484.
5. Y. C. Martin, Let's not forget tautomers, J. Comput.-Aided Mol. Des., 2009, 23(10), 693–704,  DOI:10.1007/s10822-009-9303-2.
6. R. A. Sayle, So you think you understand tautomerism?, J. Comput.-Aided Mol. Des., 2010, 24(6–7), 485–496.
7. S. Hejazi, O. Osman, A. Alyoubi, S. Aziz and R. Hilal, The Thermodynamic and Kinetic Properties of 2-Hydroxypyridine/2-Pyridone Tautomerization: A Theoretical and Computational Revisit, Int. J. Mol. Sci., 2016, 17(11), 1893,  DOI:10.3390/ijms17111893.
8. P. J. Taylor and L. Antonov, “Triage” for Tautomers: The Choice between Experiment and Computation, in Tautomerism, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2016, pp. 11–34, http://doi.wiley.com/10.1002/9783527695713.ch2,  DOI:10.1002/9783527695713.ch2.
9. O. Wahl and T. Sander, Tautobase: An Open Tautomer Database, J. Chem. Inf. Model., 2020, 60(3), 1085–1089,  DOI:10.1021/acs.jcim.0c00035.
10. D. K. Dhaked, L. Guasch and M. C. Nicklaus, Tautomer Database: A Comprehensive Resource for Tautomerism Analyses, J. Chem. Inf. Model., 2020, 60(3), 1090–1100,  DOI:10.1021/acs.jcim.9b01156.
11. P. Pospisil, P. Ballmer, L. Scapozza and G. Folkers, Tautomerism in Computer-Aided Drug Design, J. Recept. Signal Transduction, 2003, 23(4), 361–371,  DOI:10.1081/RRS-120026975.
12. M. T. Geballe, A. Geoffrey Skillman, A. Nicholls, J. Peter Guthrie and P. J. Taylor, The SAMPL2 blind prediction challenge: introduction and overview, J. Comput.-Aided Mol. Des., 2010, 24(4), 259–279.
13. Y. C. Martin, Overview of the perspectives devoted to tautomerism in molecular design, J. Comput.-Aided Mol. Des., 2010, 24(6–7), 473–474.
14. S. M. Kast, J. Heil, S. Güssregen and K. F. Schmidt, Prediction of tautomer ratios by embedded-cluster integral equation theory, J. Comput.-Aided Mol. Des., 2010, 24(4), 343–353,  DOI:10.1007/s10822-010-9340-x.
15. R. F. Ribeiro, A. V. Marenich, C. J. Cramer and D. G. Truhlar, Prediction of SAMPL2 aqueous solvation free energies and tautomeric ratios using the SM8, SM8AD, and SMD solvation models, J. Comput.-Aided Mol. Des., 2010, 24(4), 317–333.
16. A. Klamt and M. Diedenhofen, Some conclusions regarding the predictions of tautomeric equilibria in solution based on the SAMPL2 challenge, J. Comput.-Aided Mol. Des., 2010, 24(6–7), 621–625.
17. P. Nagy, Theoretical Calculations on the Conformational/Tautomeric Equilibria for Small Molecules in Solution, Biochem. Pharmacol.: Open Access, 2014, 03(02), 1–23,  DOI:10.4172/2167-0501.s4-001 , http://www.omicsgroup.org/journals/theoretical-calculations-on-the-conformationaltautomeric-equilibria-for-small-molecules-in-solution-2167-0501.S4-001.php?aid=12401.
18. K. K. Irikura. Computational Thermochemistry, ed. K. K. Irikura and D. J. Frurip, American Chemical Society, Washington, DC, 1998, vol. 677 of ACS Symposium Series, https://pubs.acs.org/doi/book/10.1021/bk-1998-0677,  DOI:10.1021/bk-1998-0677.
19. H. X. Zhou and M. K. Gilson, Theory of Free Energy and Entropy in Noncovalent Binding, Chem. Rev., 2009, 109(9), 4092–4107,  DOI:10.1021/cr800551w.
20. C. Y. Lin, E. I. Izgorodina and M. L. Coote, How Accurate Are Approximate Methods for Evaluating Partition Functions for Hindered Internal Rotations?, J. Phys. Chem. A, 2008, 112(9), 1956–1964.
21. C. E. A. Chang, W. Chen and M. K. Gilson, Ligand configurational entropy and protein binding, Proc. Natl. Acad. Sci. U. S. A., 2007, 104(5), 1534–1539.
22. A. V. Marenich, J. Ho, M. L. Coote, C. J. Cramer and D. G. Truhlar, Computational electrochemistry: prediction of liquid-phase reduction potentials, Phys. Chem. Chem. Phys., 2014, 16(29), 15068–15106.
23. J. S. Smith, O. Isayev and A. E. Roitberg, ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost, Chem. Sci., 2017, 8(4), 3192–3203,  10.1039/C6SC05720A.
24. O. T. Unke and M. Meuwly, PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges, J. Chem. Theory Comput., 2019, 15(6), 3678–3693.
25. K. T. Schütt, H. E. Sauceda, P. J. Kindermans, A. Tkatchenko and K. R. Müller, SchNet - A deep learning architecture for molecules and materials, J. Chem. Phys., 2018, 148(24), 241722,  DOI:10.1063/1.5019779.
26. P. O. Dral, Quantum Chemistry in the Age of Machine Learning, J. Phys. Chem. Lett., 2020 mar, 11(6), 2336–2347,  DOI:10.1021/acs.jpclett.9b03664.
27. A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A: At., Mol., Opt. Phys., 1988 sep, 38(6), 3098–3100,  DOI:10.1103/PhysRevA.38.3098.
28. C. Lee, W. Yang and R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37(2), 785–789,  DOI:10.1103/PhysRevB.37.785.
29. R. A. Kendall, T. H. Dunning and R. J. Harrison, Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions, J. Chem. Phys., 1992, 96(9), 6796–6806,  DOI:10.1063/1.462569.
30. W. J. Hehre, R. Ditchfield and J. A. Pople, Self—Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian—Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules Published by the AIP Publishing Articles you may be interested in Self - consistent molecular or, J. Chem. Phys., 1972, 56(5), 2257–2261.
31. M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. DeFrees and J. A. Pople, Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements, J. Chem. Phys., 1982, 77(7), 3654–3665,  DOI:10.1063/1.444267.
32. A. V. Marenich, C. J. Cramer and D. G. Truhlar, Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions, J. Phys. Chem. B, 2009, 113(18), 6378–6396.
33. J. Ho, Are thermodynamic cycles necessary for continuum solvent calculation of pK_as and reduction potentials?, Phys. Chem. Chem. Phys., 2015, 17(4), 2859–2868.
34. P. Kenny, The Prediction of Tautomer Preference in Aqueous Solution (Version 1.0), 2019, 7, https://figshare.com/articles/preprint/The_Prediction_of_Tautomer_Preference_in_Aqueous_Solution_Version_1_0_/8966276,  DOI:10.6084/m9.figshare.8966276.v1.
35. I. Soteras, M. Orozco and F. J. Luque, Performance of the IEF-MST solvation continuum model in the SAMPL2 blind test prediction of hydration and tautomerization free energies, J. Comput.-Aided Mol. Des., 2010, 24(4), 281–291,  DOI:10.1007/s10822-010-9331-y.
36. N. Tielker, L. Eberlein, G. Hessler, K. F. Schmidt, S. Güssregen and S. M. Kast, Quantum–mechanical property prediction of solvated drug molecules: what have we learned from a decade of SAMPL blind prediction challenges?, J. Comput.-Aided Mol. Des., 2020, 0123456789,  DOI:10.1007/s10822-020-00347-5.
37. S. Schlund, E. M. Basílio Janke, K. Weisz and B. Engels, Predicting the tautomeric equilibrium of acetylacetone in solution. I. The right answer for the wrong reason?, J. Comput. Chem., 2010, 31(4), 665–670,  DOI:10.1002/jcc.21354.
38. A. Nicholls, S. Wlodek and J. A. Grant, The SAMP1 Solvation Challenge: Further Lessons Regarding the Pitfalls of Parametrization, J. Phys. Chem. B, 2009 Apr, 113(14), 4521–4532.
39. D. A. Rufa, H. E. Bruce Macdonald, J. Fass, M. Wieder, P. B. Grinaway, A. E. Roitberg, O. Isayev and J. D. Chodera, Towards chemical accuracy for alchemical free energy calculations with hybrid physics-based machine learning/molecular mechanics potentials, bioRxiv, 2020 DOI:10.1101/2020.07.29.227959.
40. S. L. J. Lahey and C. N. Rowley, Simulating protein-ligand binding with neural network potentials, Chem. Sci., 2020, 11(9), 2362–2368,  10.1039/c9sc06017k.
41. J. S. Smith, B. T. Nebgen, R. Zubatyuk, N. Lubbers, C. Devereux, K. Barros, S. Tretiak, O. Isayev and A. E. Roitberg, Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning, Nat. Commun., 2019, 10(1), 1–8,  DOI:10.1038/s41467-019-10827-4.
42. A. Karton, Thermochemistry of Guanine Tautomers Re-Examined by Means of High-Level CCSD(T) Composite Ab Initio Methods, Aust. J. Chem., 2019, 72(8), 607–613,  DOI:10.1071/CH19276.
43. M. R. Shirts and J. D. Chodera, Statistically optimal analysis of samples from multiple equilibrium states, J. Chem. Phys., 2008, 129(12), 124105.
44. J. P. Ebejer, G. M. Morris and C. M. Deane, Freely available conformer generation methods: How good are they?, J. Chem. Inf. Model., 2012, 52(5), 1146–1158,  DOI:10.1021/ci2004658.
45. T. F. Havel, Distance Geometry: Theory, Algorithms, and Chemical Applications, Encyclopedia of Computational Chemistry, 2002.
46. F. Neese, The ORCA program system, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2(1), 73–78.
47. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132(15) DOI:10.1063/1.3382344.
48. S. Grimme, Supramolecular binding thermodynamics by dispersion-corrected density functional theory, Chem.–Eur. J., 2012, 18(32), 9955–9964,  DOI:10.1002/chem.201200497.
49. Jmol: an open-source Java viewer for chemical structures in 3D, accessed 10-10-2020, http://www.jmol.org/.
50. A. Hagberg, P. Swart, and D. S. Chult, Exploring Network Structure, Dynamics, and Function using NetworkX, in Proceedings of the 7th Python in Science Conference, ed. G. Varoquaux, T. Vaught and J. Millman, Pasadena, CA, USA, 2008, pp. 11–15.
51. J. Ho and M. Z. Ertem, Calculating Free Energy Changes in Continuum Solvation Models, J. Phys. Chem. B, 2016, 120(7), 1319–1329,  DOI:10.1021/acs.jpcb.6b00164.
52. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson and C. J. Carey, et al., SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 2020, 17, 261–272,  DOI:10.1038/s41592-019-0686-2.
53. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard and J. B. Maronsson, et al., The atomic simulation environment—a Python library for working with atoms, J. Phys.: Condens. Matter, 2017, 29(27), 273002.
54. B. Leimkuhler and C. Matthews, Molecular Dynamics: With Deterministic and Stochastic Numerical Methods, Springer, 2015.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sc01185e

‡ Current address: Relay Therapeutics, Cambridge, MA 02139, USA.

§ https://qcarchive.molssi.org/apps/ml_datasets/

¶ https://aiqm.github.io/torchani/examples/nnp_training.html

---