Dummy Atoms in Alchemical Free Energy Calculations

In calculations of relative free energy differences, the number of atoms of the initial and final states is rarely the same. This necessitates the introduction of dummy atoms. These placeholders interact with the physical system only by bonded energy terms. We investigate the conditions necessary so that the presence of dummy atoms does not influence the result of a relative free energy calculation. On the one hand, one has to ensure that dummy atoms only give a multiplicative contribution to the partition function so that their contribution cancels from double-free energy differences. On the other hand, the bonded terms used to attach a dummy atom (or group of dummy atoms) to the physical system have to maintain it in a well-defined position and orientation relative to the physical system. A detailed theoretical analysis of both aspects is provided, illustrated by 24 calculations of relative solvation free energy differences, for which all four legs of the underlying thermodynamic cycle were computed. Cycle closure (or lack thereof) was used as a sensitive indicator to probing the effects of dummy atom treatment on the resulting free energy differences. We find that a naive (but often practiced) treatment of dummy atoms results in errors of up to kBT when calculating the relative solvation free energy difference between two small solutes, such as methane and ammonia. While our analysis focuses on the so-called single topology approach to set up alchemical transformations, similar considerations apply to dual topology, at least many widely used variants thereof.


General considerations
Since we are concerned with classical mechanical systems, the contribution of the kinetic energy to the partition function can always be separated, and its contribution to double free energy differences cancels. We start by exploring under what conditions the configurational partition function Z of a molecule to which dummy atoms are attached, e.g., solutes L1 D and L2 D in Fig. 1(b) of the main manuscript, can be written as the product Z(L D ) = Z(L) × Z(D). L denotes the physical molecule and the superscript D indicates the presence of dummy atoms. Clearly, if the partition function can be written as a product, any dummy atom contribution cancels from the relative free energy differences of interest (cf. Fig. 1 of the main manuscript). If, on the other hand, there is no such separability of the partition function, or, in other words, if there is coupling between the energy terms of the physical system and the (bonded) energy terms involving dummy atoms, the equilibrium geometry as well as the dynamics of the physical atoms may be different in the presence and absence of dummy atoms. The physical atoms interact with the remainder of the system (water, protein etc.); these interactions depend on the (average) geometry of these atoms, as well as their dynamics. In the context of a thermodynamic cycle, these interactions take place in two different environments, e.g., gas phase and aqueous solution (calculation of relative solvation free energies), or in aqueous solution and in the protein (relative binding free energy calculations). Thus, in the presence of coupling, the identity ∆A 2 − ∆A 1 = ∆A 2 − ∆A 1 (cf. Fig. 1 of the main manuscript) may not hold.
We start with the potential energy function of some molecule L D in which one or more dummy atoms are present. Depending on the application our molecule could be in gas phase, aqueous solution or bound to a protein, we will refer to this as environment (E). The potential energy function can be schematically written as In Eq. 1 the term U LE encompasses all bonded and non-bonded interactions within the physical molecule and the non-bonded interactions between L and E. The second term U D comprises all interactions in which dummy atoms participate. By employing the notation b LD (r L , r D ) and b D (r D ), we emphasize that dummy atoms interact only through bonded terms. Since r L appears in both terms on the right hand side, no separation of the configurational partition function is possible in Cartesian coordinates.
The desired factorization can be accomplished by a partial change of variables from Cartesian to suitable internal coordinates {b LD , b D }. Adapting steps outlined, e.g., in Refs.

1-4, we obtain
Here |Ĵ LD,D (b LD , b D )| is the Jacobian resulting from the coordinate transformation. It is independent of r E and r L because of the choice of {b LD , b D } defining the position of the dummy atoms relative to the coordinates of the physical atoms. Eq. 2 has the desired form However, the factorization of Eq. 2 into a term for the physical system and a term comprising the dummy atom contributions depends on several prerequisites. In general, the steps leading to Eq. 2 assume that one uses exactly three internal degrees of freedom and associated force field terms per (dummy) atom; 1-4 Herschbach, Johnston and Rapp refer to these as non-redundant degrees of freedom. 1 E.g., given four non-linear, unbranched atoms, for which the Cartesian coordinates of atoms 1-3 are known, the position of atom 4 can be specified in terms of the bond distance d 43 , the angle θ 432 and the dihedral angle ϕ 4321 .
Any additional force field term acting on 4 constitutes a redundant degree of freedom. As already mentioned, in modern force fields an energy term is assigned to each valence and dihedral angle formed by the covalent bonds; this means that there are usually more than three bonded energy terms acting on an atom. Thus, in most practical cases the three nonredundant internal coordinates {b LD , b D } employed in the coordinate transformation leading to Eq. 2 are a subset of the bonded energy terms present. Consequently, if widespread practice is followed and the bonded energy terms attaching a dummy atom to a physical atom are the same as in the corresponding native state, then the position and orientation of the dummy atom relative to the physical system depends on (many) more than three degrees of freedom (bonded energy terms).
When investigating how the presence of such redundant terms affects the separability of the partition function, one has to distinguish two cases. Any redundant degrees of freedom (bonded energy terms), which depend only on positions (coordinates) of dummy atoms, are of no concern; their contribution to the partition function can be factored even in Cartesian coordinates. In fact, within large groups of dummy atoms all bonded interactions should be kept to maintain its structural integrity. By contrast, care is required in the "junction region", i.e., for redundant bonded terms, which involve coordinates of physical atoms as well as dummy atoms. It is these cases we will analyze in the following subsections.
In biomolecular simulations bonds are frequently held rigid by holonomic constraints to increase the integration time-step. As will be shown in the following subsections, complications from non-redundant degrees of freedom result from angle bending or dihedral angle terms, which are not subject to constraints. Therefore, all considerations given below apply regardless whether bonds are flexible or held rigid by SHAKE, 5 RATTLE 6 or similar means.  to water is completed, no intramolecular non-bonded interactions are present. Since we assume the molecules to be perfectly planar, the only force field terms needed are the usual quadratic bond stretching and angle bending terms, k(x − x 0 ) 2 = k(∆x) 2 , where ∆x denotes the displacement of a bond length r or an angle θ from its equilibrium value.
In 2D there are 2N − 3 non-redundant intramolecular degrees of freedom, i.e., in our example (six atoms, three physical, three dummy atoms) there are 9 non-redundant degrees of freedom. We start with the potential energy functions of physical water and the dummy part, respectively, where the subscripts refer to the atom labels of Fig. 1. Having used 6 out of 9 non-redundant intramolecular degrees of freedom, we have available three energy terms to attach the dummy atoms to the physical water part. Specifically, we anchor D C 2 relative to water and O relative to the dummy part. In 2D, in addition to the bond r O-D C 2 , one angle per side is sufficient; no dihedral angle is needed: Thus, using Eqs. 4 and 5 the full energy function for the WAT-D endpoint is To illustrate the steps leading to Eq. 2 in the general treatment, we use Eq. 6 and explicitly write the partition function for the WAT-D end state where the square brackets indicated the two multiplicative factors. Cartesian coordinates are kept for the physical system (water), but for the dummy part of WAT-D we transform to the internal coordinates present in the bonded energy terms of Eq. 6. The Jacobian factors r A-B resulting from the transformation from Cartesian to planar polar coordinates positioning atom B relative to atom A are written explicitly; one sees that they do not depend on degrees of freedom of the physical water molecules.
At the corresponding ETH endpoint, two additional angle bending terms, Adding a redundant degree of freedom is equivalent to introducing a geometrical constraint, which in this specific case is given by Adding 2 to U W AT −D (Eq. 6) and utilizing Eq. 8, we obtain for the potential energy function of WAT-D On the right hand side of Eq. 9 the redundant energy term has been expressed in terms coupling with the physical degrees of freedom (water) is introduced and the partition function of the WAT-D endpoint following from Eq. 9 remains separable. While the additional 2 term will affect single free energy differences, it will cancel exactly from double free energy differences between two environments and can safely be included as a force-field term to keep the dummy atoms attached to the physical system. An important observation here is that this can be deduced directly from the constraint equation Eq. 8, which does not contain any terms involving degrees of freedom of the physical molecule.
Explicitly writing expressions for the potential energy function, such as Eq. 9 above, is not needed, which simplifies analyses considerably.
Let us repeat the above steps for the second redundant degree of freedom, the θ H 2 -O-D C 2 angle, one would potentially like to keep. In this case, the relevant geometrical constraint Eq. 10 involves θ To summarize, the addition of an energy term in a redundant internal degree of freedom introduces a geometric constraint. If the constraint involves only internal degrees of freedom pertaining to one or more dummy atoms, then the physical part of the system is not affected, and any contribution, which the additional force field term makes to the partition function, cancels from double free energy differences. If, on the other hand the constraint couples the redundant degree of freedom with one or more physical degrees of freedom, then the physical system is affected and the partition function cannot be factorized.
2 Details concerning coupled three angles

Additional theoretical details
We outlined in the main manuscript ("Coupled three angles" in Sect. 2.1.2) that unintended coupling between degrees of freedom of dummy atoms and physical atoms arises when there are three atoms bound to a central atom, and one of these is transformed into a dummy atom. The simplest possible model for this is the alchemical transformation of ammonia to water, cf. Fig. 5 of the main manuscript. We showed that the following constraint applies to the physical bond angle where we use the shorthand The upper and lower bounds in Eq. 11 both arise in the limiting case of all four atoms being coplanar. The lower bound can be deduced immediately from Figs with The potential energy function U consists of all bonded energy terms present in water with the dummy atom attached to it. Since our integration variables are exactly the internal coordinates of these energy term, which are additive, the potential energy function U itself in this non-redundant coordinate set is unproblematic. However, the partition function Eq. 12 does not factorize because of the limits of integration and the functional form of the Jacobian factor sin Ψ, which couples all three angles.
Thus, attaching a dummy atom by a bond and two angles does not lead to a multiplicative contribution to the partition function. Two limiting cases, however, are of special interest.
If one attempts to keep all atoms, including the dummy atom, in the same plane, the hard and sin(Ψ) = 1 (14) With these simplifications Eq. 12 can be rewritten so that the partition function of the physical water molecule becomes a multiplicative factor, i.e., the desired separability of the partition function is regained. The somewhat surprising conclusion from this analysis is that whenever one anchors a dummy atom by one bond and two angles with respect to a physical molecule, one should set the equilibrium values of the two angles involving the dummy atom to 90 degrees in order to maintain separability of the partition function (Eq. 3), at least in good approximation since the dummy angles still fluctuate around 90 degrees. As shown next, this choice for the bond angles involving the dummy atom also effectively prevents flapping. Symmetric pathway: Y flaps through the midpoint of X 1 , X 2 and D or, equivalently, D crosses the X 1 -Y-X 2 plane following a path distant to X 1 and X 2 . (c) Canonical asymmetric path: D takes the path between X 1 and X 2 to cross the X 1 -Y-X 2 plane. The dash in θ 3 indicates that it is chosen as 0 ≤ θ 3 ≤ π, so θ 1 + θ 2 + θ 3 = 2π. For usual geometries encountered in chemistry, θ 1 and θ 2 are acute angles (in contrast to case (b)), i.e., 0 ≤ θ 1 , θ 2 ≤ π/2. If unusual geometries are permitted, one (but only one) of them may be obtuse. (d) Non-canonical asymmetric path I: X 1 takes the path between D and X 2 to cross the D-Y-X 2 plane. (e) Non-canonical asymmetric path II: same as (c) and (d), but X 2 takes the path between D and X 1 to cross the D-Y-X 1 plane.

Energy barrier of inversion at a pyramidal center
We just showed that attaching a dummy atom through a bond stretching and two bond angle terms to a physical system leads to a three angle constraint (two dummy angles, one physical angle), which potentially influences the physical system. Somewhat unexpectedly, the effect of this constraint can be (mostly) avoided by setting the equilibrium value of the two angle terms involving the dummy atom to 90 • . In Sect. 2.2 of the main manuscript we showed that positioning a dummy atom relative to the physical system through a bond stretching and two bond angle terms is prone to flapping, as this choice of internal coordinates leads to two sets of Cartesian coordinates, which are equivalent with respect to their internal coordinates for the dummy atom. In the main manuscript, our example was the alchemical transformation of ammonia to water, with one of the ammonia hydrogens becoming a dummy atom attached to the water oxygen. This dummy atom can be either "above" or "below" the plane formed by the water atoms. When the dummy atom switches between the two possible positions, it has to traverse this plane, i.e., all four atoms are in the same plane, exactly the configuration for which the effect of the three angle constraint is most pronounced.
While the ammonia to water transformation may be considered a fringe case, positioning a dummy atom or group relative to a physical system through one bond and two bond angles is an appealing option. It is, therefore, necessary to understand the pathway of flapping which may arise in this case. Special emphasis has to be paid to the case where the equilibrium angles involving dummy atoms have been set to 90 • to mitigate the effect of the three angle constraint. Further, at the ammonia end state, "flapping" is normal since this is the nitrogen inversion occurring for ammonia and amines. Conversely, at the water end state (assuming a flexible water molecule), we would like to keep the dummy atom in a fixed position relative to the water molecule to prevent the planar transition state, which through a three angle constraint would influence the angle bending motion of water. For the toluene to pyridine transformation one has (see Fig. 9) D C= D, N=Y , C 1= X 1 and C 5= X 2 . Similarly, for acetone to 2-propenol (see Fig. 12), the correspondence is D H 13= D, Figs. 2(b)-(e) illustrate the four possible pathways, along which "pyramidal flapping" may occur. All involve a configuration in which the four atoms are coplanar. Fig. 2(b) shows the symmetric path, where D passes through the X 1 -Y-X 2 plane while staying at the opposite side of atoms X 1 and X 2 . Alternatively, this path can be viewed as the central atom Y tunneling through the plane formed by X 1 , X 2 and D. This latter view, of course, corresponds to nitrogen inversion. While this is the expected pathway for normal molecular geometries, we have to keep two aspects in mind. First, we are discussing dummy atom dynamics here, which we want to keep in a well defined position with respect to a physical molecule. Second, because we want to mitigate the effect of the associated three angle constraint on the physical system, we also have to consider the unusual case where angles The first possible asymmetric path is illustrated in Fig. 2(c): D crosses the plane X 1 -Y-X 2 between X 1 and X 2 . We refer to it as canonical asymmetric, as it is the dummy atom which crosses the plane. The second asymmetric path, shown in Fig. 2(d), is obtained if X 1 crosses the D-Y-X 2 plane between D and X 2 . In the third asymmetric path (Fig. 2(e)), X 2 crosses the D-Y-X 1 plane between D and X 1 .
To discriminate between these four potential pathways, we need to calculate the minimum energy of the planar configurations depicted in Figs. 2(b)-(e), which correspond to the respective transition state. The height of the energy barrier a particle has to cross along the four pathways also indicates how likely it is that flapping will occur. We start with the symmetric planar configuration shown in Fig. 2(b) and carry out an energy minimization under the constraint of planarity, i.e., the sum of the three angles has to equal 2π. The only energy terms we need to consider are the three angle bending terms; changes in bond lengths do not affect the angles. Due to non-bonded exclusions, there are neither Lennard-Jones, nor electrostatic interactions present between the physical atoms X 1 , X 2 , Y, and obviously none for D. Assuming quadratic potentials, the potential energy term including the Lagrangian multiplier is: We want to find the angles θ i ; θ 0 i and k i are the equilibrium angles and force constant of the respective force field terms. Taking the derivatives and setting them to zero leads to the following system of equations with solutions To find the energy of the transition state, we insert the expressions for the three angles into the energy term of the Lagrangian multiplier equation and obtain Next, we turn to the asymmetric cases Figs. 2(c)-(e). Because of our choice of labeling, the same derivation can be applied to all three cases. Here, the Lagrangian multiplier equation is: Because of the energy function for angle bending terms, we need to introduce the auxiliary variables θ i . Specifically, the instantaneous value of an angle in an angle bending energy term is expected to be in the range 0 ≤ θ i ≤ π. However, to define the constraint for the coplanar configuration of interest, we also need to allow for angles π < θ i < 2π. To proceed, we assume that Importantly, θ 3 is always the angle formed by the two atoms (plus central atom Y), between which the respective third atom crosses the plane (θ 3 is its complement to 2π). E.g., in Fig. 2(c), θ 3 describes the angle formed by atoms X 1 and X 2 (with Y as center), and atom D crosses the X 1 -Y-X 2 plane between them. In Fig. 2(d), X 1 crosses the D-Y-X 2 plane between atoms D and X 2 etc. Using this notation, the system of equations to solve is: which leads to Having found the θ , we can insert them in the energy function of the Lagrangian multiplier equation and obtain   This follows from the observation that the minimum energy for all three asymmetric flaps can be obtained by permuting indices 1, 2 and 3 in Eq. 25, while keeping the association of the force-field equilibrium angles θ 0 and force-constants k fixed. E.g., index 3 always refers to the physical angle degree of freedom. Doing so is not in line with the labeling used in Figs. 2(d,e), the utility of which was mostly unification of the mathematical treatment. As expected, the symmetric flap, corresponding to nitrogen inversion in the case of ammonia, has a far lower energy barrier than the hypothetical asymmetric paths. One sees in Fig. 3 that the height of this barrier depends only weakly on the value of the force constants k 1 and k 2 , i.e., extremely high values of the force constants for the two angle bending terms involving the dummy atom would be necessary to prevent flapping.

Error analysis in thermodynamic integration
Free energy differences were calculated using Thermodynamic Integration (TI), 10,11 as well as Bennett's acceptance ratio (BAR). 12 The results obtained with the latter were computed as a control; for additional details see Sect. 5 below. Each of our 5 × 21 simulations (N = 5 repetitions, K = 21 λ-states) per system and transformation are statistically independent.
Using BAR, each of the individual ∆A λ i →λ i+1 , (i = 1, . . . K − 1) was, therefore, obtained from five independent simulations; similarly, in TI we have five values for each of the K derivatives ∂U ∂λ λ i (i = 1, . . . K). In the case of BAR, Gaussian error propagation based on the standard deviation of the individual ∆A λ i →λ i+1 is straightforward; in the following, we outline how we estimated the overall statistical error using TI.
In TI, the integral over λ has to be approximated by numerical quadrature. Numerical integration algorithms, applied to TI, can formally always be written as where ∂U ∂λ λ i ,r is the potential energy derivative of the ith λ-state averaged over the rth of our N = 5 independent simulations and K = 21 denotes the number of simulated λstates. The W λ i are constants depending on the quadrature scheme. If they are known (as, e.g., when using the trapezoidal rule, Simpson's rule, Gauss-Legendre or Clenchaw-Curtis integration), error propagation based on the standard deviation of the individual ∂U ∂λ λ i is straightforward.
In this work, however, we fitted the integrand using natural cubic splines, approximating it in a piece-wise continuous manner by cubic polynomials, which were then integrated analytically. 13,14 Denoting this spline-fit-then-integrate approach as splint[], one can show that 13 i.e., splint[] can be regarded as any other numerical quadrature method. However, obtaining the coefficients W λ i is not as obvious as in standard numerical quadrature methods. Because of the linearity of the weighted sum in Eq. 27, we find for the free energy difference ∆G, i.e., the average over the N = 5 independent results ∆G r : The corresponding standard deviation σ is The coefficients W λ i can be extracted by means of custom code from spline routine used; this is, e.g., done in the widely-used software package Alchemical Analysis. 15 We took a different, yet mathematically equivalent approach, which requires neither knowledge of implementation dependent representations of the coefficients W λ i , nor of the computational details of the spline routine used. Specifically, we used the Python package SciPy (www.scipy.org), which wraps the FITPACK 16 library, but any numerical package capable of fitting and integrating natural cubic splines could have been chosen.
In order to compute the standard deviation σ (Eq. 29) without explicit knowledge of the Because Eq. 27 is a linear sum of the values of the integrand ∂U ∂λ λ i ,r at the λ-states, most of the terms ∂U ∂λ λ i cancel from Eq. 30 and one obtains  First, the array dudl aver is filled with the averages over replicas ∂U ∂λ λ i for each state λ i from the data ∂U ∂λ λ i ,r stored in dudl. The estimate for ∆G, i.e. result average, is calculated as the spline integral over dudl aver; cf. Eq. 29. This is followed by the estimation of the standard deviation s λ i for each λ-state; cf.
Eq. 31. For each λ i , the respective entry in dudl aver[λ] is replaced by each of the individual values ∂U ∂λ λ i ,r from dudl[λ] [r]. The resulting array, labeled dudl tmp, is evaluated with splineintegrate() and the result stored in spl int [r]. The latter is used to compute the standard deviation s λ i for each of the λ i .
The overall error estimate is obtained by Gaussian error propagation. For further details, the reader is referred to the input scripts and source code available at Zenodo (https://zenodo.org/record/4381708). 17 5 Detailed specifications for each transformation studied Full technical details and additional results are summarized on one page per transformation for each of the 12 relative solvation free energy differences reported in the main manuscript.
For each alchemical pair, there is a figure showing the atom labels used. The accompanying table lists on the left hand side gas phase and solvation free energy differences obtained along the absolute and relative paths. If a relative free energy difference was computed with more than one approach, then individual results are given for each of them.
The right hand side of each table for the twelve alchemical transformations lists force field terms involving dummy atoms, which were either deleted and/or modified. Illustrative examples of such lists were given in Table 3 of the main manuscript. If more than one approach was used to carry out the alchemical transformation, this is indicated by footnotes.
Urey-Bradley terms involving dummy atoms which would lead to a constraint were always deleted; this is not explicitly noted.
We report free energy differences obtained with thermodynamic integration (TI), 10 as well as Bennett's acceptance ratio (BAR) method 12 (see below). Tables 4 and 5 of the main manuscript were compiled from the TI results reported here. All free energy differences were also calculated using BAR based on the same underlying MD simulations used for TI. Coordinates were saved to disk every 50 steps (cf. Methods of main manuscript); each trajectory saved at λ i was used to compute energies at λ i−1 , λ i and λ i+1 . Thus, the overall free energy difference of a transformation was calculated as the sum of 20 individual free energy differences; except for PRP2DIM-3, where we employed 39 λ-states and hence the total free energy difference is the sum of 38 contributions.
As described in the main manuscript, each simulation was repeated five times (by starting from different random initial velocities); hence, for each free energy difference between λ i and λ i+1 we computed five values ∆A j λ i →λ i+1 (j = 1, . . . 5). Rather than just averaging over these (i.e., ∆A λ i →λ i+1 = 1/5 5 j=1 ∆A j λ i →λ i+1 ), we combined the raw data from the respective simulations (forward and backward energy differences), and computed our best expectation value for the free energy difference ∆A E λ i →λ i+1 from these combined data. The statistical uncertainty for this step was estimated as the standard deviation of the five individual ∆A j λ i →λ i+1 values. Thus, using BAR, the total free energy difference was estimated as ∆A E 0→1 = 20 k=1 ∆A E λ i →λ i+1 . The corresponding error estimate was obtained by Gaussian error propagation from the variances σ 2 (∆A λ i →λ i+1 ) of the individual steps.  Deleted:  Deleted:  2.38 ± 0.01 -2.38 ± 0.01 0.00 ± 0.00 Water TI -6.89 ± 0.01 6.89 ± 0.01 0.00 ± 0.00 BAR -6.89 ± 0.01 6.89 ± 0.01 0.00 ± 0.00 Water a TI -6.88 ± 0.01 6.88 ± 0.01 0.00 ± 0.00 BAR -6.89 ± 0.01 6.89 ± 0.01 0.00 ± 0.00 Methane b → Water a ∆∆G abs solv TI -9.27 ± 0.02 BAR -9.27 ± 0.02 a instead of the 1 fs integration step used in all other simulations, due to the small physical molecule size, a 0.25 fs integrator step was used here; b with 567 water molecules, see methane-ammonia); c all to 90 • equilibrium angle and 100 kcal/rad 2 force constant. Deleted:

Toluene-Pyridine
Modified: Modified: Modified a,c :   c with 567 water molecules; d all to 111.749 • equilibrium angle and 3.55 kcal/rad 2 force constant.

Dual Topology junctions
6 Comparison of absolute solvation free energies from simulation with experiment Absolute solvation free energies for small molecules calculated with the cgenff force field family [7][8][9] are not readily available in the literature; therefore, we list them here for our model compounds.
We compare these to the the experimental values compiled in the Minnesota Solvation Database, 18 as well as to recently reported computational results. 19,20 For the four systems for which a comparison to Ref. 19 is possible, the agreement is excellent. By contrast, the numbers reported by Huenenberger and co-workers 20 differ systematically. However, one should keep in mind the following differences in simulation setup. Calculations were carried out with GROMACS rather than CHARMM. Lennard-Jones interactions were truncated with a force switching function rather than switching the potential energy used in this work and in Ref. 19. A long-range correction for Lennard-Jones interactions was applied, which was not done here since our primary focus was closing the thermodynamic cycle (Fig. 1b of the main manuscript).  . For the reasons discussed in the main text, the integrand varies rapidly near λ = 1, which leads to a systematic error during numerical quadrature. The plot is shown for the gas phase, but the behavior is almost identical in aqueous solution, leading to fortuitous error cancellation for PRP2DIM-2 using TI for ∆∆A solv ; cf. the main manuscript.