A New Formulation for the Concerted Alchemical Calculation of van der Waals and Coulomb Components of Solvation Free Energies

Alchemical free energy calculation via molecular dynamics has been widely used to obtain thermodynamic properties related to protein-ligand binding and solute-solvent interaction. Although soft-core modeling is the most common approach, the Linear Basis Function (LBF) methodology [Naden et al., J. Chem. Theory Comput., 10(3), p. 1128, 2014&11(6), p. 2536, 2015] emerges as a suitable alternative. It overcomes the endpoint singularity of the scaling method while maintaining essential advantages such as ease of implementation and high flexibility for post-processing analysis. In the present work, we propose a simple LBF variant and formulate an efficient protocol for evaluating van der Waals and Coulomb components of an alchemical transformation in tandem, in contrast to the prevalent, sequential evaluation mode. To validate our proposal, which results from a careful optimization study, we performed solvation free energy calculations and obtained octanol-water partition coefficients of small organic molecules. Comparisons with results obtained via sequential mode using either another LBF approach or the soft-core model attest to the effectiveness and correctness of our method. In addition, we show that a reaction field model with an infinite dielectric constant can provide very accurate hydration free energies when used instead of a lattice-sum method to model solute-solvent electrostatics.


Introduction
In recent years, calculations of free energy differences have been performed in several investigations to provide values of important thermodynamic properties. 1 Binding free energies of protein-ligand systems, for example, are applied in the discovery of new drugs through the prediction of structural stability and enzymatic activity. [2][3][4] Solvation free energy calculations, in turn, provide information on solubility and partition coefficients. [5][6][7][8][9][10] These properties are widely studied via molecular simulation, [11][12][13] with which the microscopic behavior of a system can be investigated.
Using Molecular Dynamic (MD), we can simulate alchemical transformations, 1,2 meaning that the free energy difference between two states is computed using a non-physical pathway. This is possible because free energy is a state function, i.e., it depends only on the initial and final states of a process. Intermediate states often connect the two extremes of the process so that the alchemical transformation occurs gradually. Several such states are usually required to continuously cover the relevant region of phase space. The number of such states is essentially dependent on the choice of the alchemical pathway.
The states are defined by different values of a coupling parameter λ, which controls the strength of the potential energy of interaction between particles of interest. For example, in calculating the solvation free energy, the initial state is a solute in a vacuum, and the final state is the same solute in a solvent medium. The system potential energy can depend on λ in such a way that the solute-solvent interaction is nullified when λ = 0 and fully activated when λ = 1. 14 Different formulations for the potential energy as a function of λ have been created, resulting in distinct alchemical pathways.
The scaling method 15 consists in multiplying the original potential of the alchemically affected interactions by a λ-dependent switching function. However, this method can lead to imprecise free energies in systems involving Lennard-Jones (LJ) interactions due to an endpoint singularity. The LJ potential goes to infinity as the distance between atoms approaches zero, which becomes an issue when λ changes from zero to a small finite value. To solve this shortcoming, Beutler et al. 16 proposed their well-known soft-core approach.
The soft-core model is widely applied in alchemical calculations [17][18][19][20] and has been implemented in MD software packages such as GROMACS, 21 CHARMM, 22 and AMBER. 23 This approach provides a smooth free energy profile as a function of λ and yields a low-variance estimate. 24 However, an important disadvantage stems from the modification done in the mathematical expressions of the LJ and Coulomb potentials to include λ. It involves changes in the inner force loop of an MD code, 25,26 which is usually highly optimized in terms of specialized hardware utilization. Moreover, multiplestate free energy estimators, such as the Multistate Bennett Acceptance Ratio (MBAR) method, 27 require computing the potential energy of each sampled configuration at all λ states under consideration in the analysis. With a soft-core potential, we must do such a calculation on the fly, which also demands changes in the inner loop of the MD code. In addition, careful planning becomes necessary, given the difficulty of including unanticipated states in the MBAR analysis.
In the scaling method, on the other hand, the LJ and Coulomb potentials remain unchanged as we multiply them by their respective switching functions. Hence, despite the endpoint singularity issue, it has considerable advantages over the softcore method. For example, we can defer the computation of potential energies at the multiple λ states to a post-processing stage. 25,28 Then, we can apply the MBAR method to estimate the free energies and other properties at either sampled or unsampled λ states via reweighting. 27,29 Such unsampled states can even be defined using switching functions other than those employed in the sim-ulations. This opens the possibility of searching for optimal switching functions, according to some criterion of choice, without performing new simulations.
Recently, Naden et al. 25,26 proposed a general framework for the scaling approach, which they refer to as the Linear Basis Function (LBF) model. They managed to avoid the endpoint singularity by softening short-range interactions 28,30 in the first step of the alchemical transformation and then applying, in a later step, a correction to restore the original LJ potential.
Here, we propose an LBF model that is mathematically simpler and involves fewer parameters than the original one. We provide a set of switching function parameters that we recommend for solvation free energy calculations in general despite having been obtained by optimizing the hydration of selected solutes. In contrast to most current methods, 17,31 our approach aims at providing concerted alchemical transformations, in which the total free energy difference, including both the van der Waals and Coulomb contributions, is obtained after a single sweep of λ values from 0 to 1.
Concerted alchemical transformations are suitable for advanced sampling schemes such as λ dynamics 32,33 and adaptive biasing. 34,35 Furthermore, they are computationally less expensive than sequential evaluation approaches, where each step is associated with different sampling requirements and statistical errors. 18 Lee et al. 18 also evaluated transformations with simultaneous changes in LJ and Coulomb potentials but they applied a modified soft-core model. In contrast, our proposal relies on the LBF approach, thus bringing about all the benefits mentioned earlier.
To test the effectiveness and robustness of our LBF method, we performed MD simulations to calculate the solvation free energy of organic molecules in both water and 1-octanol. It allowed us to compute their octanol-water partition coefficients and compare results across different methods and with experimental and simulation data from the literature.
In addition, while computing hydration-free energies, we analyzed the use of different pairwise electrostatic approaches [36][37][38] to model solutesolvent interactions in substitution to the Particle Mesh Ewald (PME) 39 method. We observed that a reaction field model 36,37 could provide highly accurate results if the dielectric constant of the medium is considered to be infinite, like in a conductor.
The results deviate slightly from PME-based ones, and, if necessary, they can be corrected with high precision through a simple perturbation method. 40 Along with the proposed LBF formulation, this observation is an important contribution to the field of alchemical free energy calculation.
The remainder of this document is organized as follows. In Sec. 2, we present the applied methodology. That includes the proposed LBF model and the schemes we used for obtaining optimized parameters. Results are presented in Sec. 3. Finally, we summarize our conclusions in Sec. 4.

Linear Basis Function Approach for Alchemical Transformations
In the Linear Basis Function (LBF) approach, formalized in a series of publications by Naden et al., 25,26 the total potential between two alchemically interacting particles at a distance r is written as where u k (r) is the potential basis function of type k and h k (λ) is its corresponding switching function. The choice of basis functions can be made so as to avoid the endpoint singularity described in the introductory section (Sec. 1). Additionally, a judicious choice of switching functions can result in an efficient alchemical protocol. The original LBF model relies on a strategy proposed by Buelens and Grubmüller 28 (BG) to eliminate the endpoint singularity. It consists in enforcing the Lennard-Jones (LJ) potential to detour from its path to infinity as r goes down to zero, but only at an energy value high enough for the change to have a negligible impact on the sampling probabilities at usual thermodynamic conditions. In this way, the modified and original models, which Buelens and Grubmüller 28 referred to as the capped and native pair potentials, respectively, yield free energy values indistinguishable from each other. However, Naden et al. 25 observed that enforcing the detour at lower potential energies produces free energy results with smaller variances but at the cost of introducing non-negligible deviations from the native-model results. As a solution, they added a new step to the free energy calculation so as to account for the residual (or decapping) contribution, defined as the difference between the native and capped potentials.
Naden et al. 25,26 employed a basis with four potential functions. Besides using the Coulomb interaction for charged atoms, they performed a WCA decomposition 41 of the LJ potential into repulsive and attractive contributions, followed by the BGtype splitting of the repulsive part. Their version of the capped repulsive potential requires some constants to be calculated for each individual set of LJ parameters. They defined the switching function corresponding to each contribution so as to reduce the total variance of the alchemical calculations. The authors took advantage of the reweighting capability of the LBF approach in order to carry out a very efficient optimization procedure. They ended up recommending the trivial, linear coupling for all contributions except the capped repulsive potential, for which they suggested a fourthdegree polynomial with tuned coefficients. They also tested different orderings for the coupling of van der Waals and electrostatic interactions. One of the best performances was achieved by switching on the attractive and the capped repulsive potentials simultaneously at first, followed by the residual potential and, finally, by the electrostatic contribution.
In the present work, we propose a simpler formulation. Instead of using the WCA decomposition, we applied a BG-type capping to the full LJ model. Our LBF potential for the alchemical interactions is written as u(r, λ) = h c (λ)u c (r) + h r (λ)u r (r) + h e (λ)u e (r), (2) where u c (r), u r (r), and u e (r) are the capped, residual, and electrostatic pair potentials, respectively.
The LJ potential is expressed as u lj (r) = ν lj (r/σ), where σ and are the usual LJ parameters, 42 and ν lj (x) is a reduced potential given by We developed a simple version of the BG potential, in which a reduced capping function ν cap (x) replaces ν lj (x) for all x < 1. Therefore, u bg (r) = ν bg (r/σ), where The capping function is defined as a polynomial of degree m satisfying the conditions ν  lj (1) for n = 0, . . . , m 1 , where m 0 + m 1 = m and f (n) (x) is the n-th order derivative of f (x). These conditions, which completely determine the polynomial coefficients, guarantee that ν bg (x) is continuous and smooth at x = 1 and singularity-free at x = 0. The polynomial that fulfills these conditions with m 0 = m 1 = 3 is Therefore, rather than diverging to infinity, the potential energy approaches 119.2 as r → 0. As an example, this value is close to 42 kJ/mol in the case of the LJ interaction between oxygen atoms in the TIP3P water model. 43 Fig. 1 shows a comparison between the native and capped potentials in their reduced forms, given by Eqs. (3) and (4), respectively. The main advantage of the capped potential defined above, in contrast to those used elsewhere, 25,26,28 is its simplicity. By changing m 0 , it is also possible to alter the potential energy value at r = 0. Nevertheless, the reduced potential in Eq. (5) succeeded in all studies performed here. In summary, with the definitions above for the capped and residual potentials, we can rewrite Eq. (2) as u(r, λ) = h c (λ)u bg (r)+ h r (λ)[u lj (r) − u bg (r)] + h e (λ)u e (r). (6) A sensible choice of switching functions h c , h r , and h e can reduce the number of λ states required to obtain free energy estimates with a desirable precision. Here, we propose a family of smooth switching functions designed to allow the coupling of all basis potentials in a single sweep of λ values from 0 to 1. They are defined as where λ s k and λ f k are, respectively, the starting and finishing points of the actual switching of the k-th basis potential, and S(x) is a twice-differentiable smoothstep function given by 18 Thus, multiple partially coupled alchemical potentials can coexist at the same λ value, such as in the example shown in Fig. 2(a). In other words, a particular coupling stage can start before the previous one has finished. An interesting reinterpretation of the coupling scheme just described can be observed in Fig. 2 (b). It stems from rewriting Eq. (6) as u(r, λ) = h bg (λ)u bg (r) + h lj (λ)u lj (r) + h e (λ)u e (r), where h bg = h c − h r and h lj = h r . Hence, at the smallest λ values, only the capped potential is coupled to the system, thus preventing the endpoint singularity. Then, when it is safe to do so, the native LJ potential progressively replaces its capped counterpart.
For comparison, we also apply the formulation of Naden et al. 25,26 to our particular splitting of the LJ potential. It consists in a sequential coupling of the capped, residual, and electrostatic potentials separately. For this, we resort to a quartic switching function 25 (10) In the initial step, we make h c = P c , while keeping h r = h e = 0. In the middle step, we make h c = 1, h r = P r , and h e = 0. Then, in the final step, we keep h c = h r = 1 and make h e = P e . This procedure is illustrated in Fig. 3, whose each panel depicts the non-constant switching function used in each step. Applying a linear scaling (i.e., A k = B k = C k = 0) for both the residual and electrostatic potentials, as seen in the figure, was the recommendation of Naden et al., 25,26 who have also applied it for the attractive part of the WCA potential. For the capped repulsive part, they recommended using Eq. (10) with A = 1.62, B = −0.889, and C = 0.0255. Here we employ these values as an initial guess for obtaining an optimal switching function (see Sec. 2.3) for our capped LJ potential.

Alchemical Coupling of Solute-Solvent Interaction
When using our LBF model for computing solvation free energies, we adopt the so-called coupling route, 14,31 in which only solute-solvent interactions undergo alchemical transformation. Hence, the total potential energy of a system as a function of the particle coordinates r and coupling parameter λ is given by where letters A and S refer to the solute and the solvent, respectively. The potential U b includes all bonded interactions (bond stretching, bending, torsion, and so on) in the system, while U aa+ss nb refers to the solute-solute and solvent-solvent nonbonded interactions, which are not affected by the alchemical transformation. We employed a switched and truncated LJ model and PME electrostatic calculations for evaluating U aa+ss nb . The alchemically affected potentials U as k are intermolecular interactions between solute and solvent. The parts associated with van der Waals interactions can be represented as where k is either LJ or BG, r c is the cutoff distance, r s is a switching distance, and S(x) is given by Eq. (8). For every pair formed by a solute atom i and a solvent atom j, respectively located at Cartesian positions r i and r j , we have that r ij = r i −r j while ij and σ ij are the pair-specific LJ parameters, usually obtained from their atomspecific counterparts via mixing rules. 42 It is important to note the presence of the reduced pair potentials ν lj (x) and ν bg (x), defined in Eqs. (3) and (4), respectively. The electrostatic part of the solute-solvent interaction is computed using a reaction field model 36,37 considering a perfectly conducting medium whose dielectric constant is infinite. In this case, the solute-solvent electrostatic potential energy is given by where CRF stands for conductor-like reaction field, q i and q j are the electric charges of atoms i and j, respectively, r is the dielectric constant of the medium, and 0 is the permittivity of free space. The reduced potential is which satisfies the conditions ν crf (1) = ν crf (1) = 0, thus ensuring continuity at r ij = r c of both the potential and its resulting force. Nevertheless, a smooth switching from r s to r c was applied so that the second-order derivative is continuous as well.
The distance between atoms throughout the simulations will never be zero, avoiding a singularity in the electrostatic potential. This necessarily occurs because the electrostatic potential is only turned on when the Lennard-Jones interactions are, at least, partially coupled, creating a repulsive layer on the particles and, thus, preventing opposite charges from collapsing by attraction. A cutoff, pairwise treatment is applied for the solute-solvent electrostatics due to the ease of implementation and computational efficiency when compared to doing separate PME calculations involving distinct groups of atoms. As we will show in Sec. 3.1 through empirical tests, simulations using U as crf provide free energy results, which are often very close to those obtained using U as pme instead. It is, in fact, the best choice for this purpose among a set of tested pairwise electrostatic models. Furthermore, one can readily apply a correction by sampling the difference U as pme − U as crf during the simulation at the final λ state and then performing a Free Energy Perturbation 40 (FEP) analysis. Hence, the required correction is estimated by means of the identity where k b is the Boltzmann constant, T is the system temperature, β = (k b T ) −1 , and · λ=1 denotes an ensemble average at the final λ state. We will also show in Sec. 3.1 that FEP calculations are sufficiently precise for doing this type of correction, thus dispensing with the need for more elaborate methods. 44

Optimization Schemes
We solved different minimization problems to obtain optimal parameters for the switching functions defined in Eqs. (7) and (10). We considered alternative objective functions, namely (i) a measure of local variations in the ∂U/∂λ vs. λ curve and (ii) the asymptotic variance of the coupling free energy.
When the sampling effort of the MD simulations is uniformly spread over the whole range from λ = 0 to λ = 1, regions with too large or too small free energy variations will result in inefficient sampling due to, respectively, poor or excessive probability distribution overlaps between adjacent states. Supposing that the amount of overlap is directly proportional to the free energy difference between two states, 45 an ideal free energy profile would be a straight line from G(0) to G(1), hence with constant ∂U/∂λ = ∆G. Since most cases of interest involve a free energy barrier, another efficient profile would entail a smooth increase from G(0) up to the barrier top, followed by a smooth decrease from there down to G(1). Although these ideal profiles are most often unrealizable for a given family of switching functions, we can look for a profile that minimizes some measure of local variations in the mean derivative ∂U/∂λ = G (λ). Here we use a finite-difference approximation for the mean absolute second-order derivative, where the λ domain has been divided into N p equally-sized intervals.
Let us now discuss variance minimization, which Naden et al. 25,26 have also applied to obtain optimal coefficients for Eq. (10) when acting as a switching function for their repulsive WCA potential. As they did, we compute the free energy variance using the expression which is related to the thermodynamic integration method. 46 In this case, we evaluate the variance of the mean derivative ∂U/∂λ and numerically integrate it with the λ range. The mean derivative ∂U/∂λ can be obtained in a post-processing stage by We apply the MBAR method to estimate both the value and the variance of the mean derivative at either sampled or unsampled states through reweighting. 27,29 Most importantly, the LBF approach allows us to perform reweighting aiming at states defined for switching functions other than those employed in the simulations. This possibility makes the optimization process very fast because evaluating the objective function for new parameter values does not require new MD simulations.
We define three optimization problems with the objective functions indicated in Eqs. (16) and (17). They will be referred to as schemes S1, S2, and Q, where S and Q stand for smooth and quartic switching functions, respectively. They all require an initial guess of switching function parameters, which we use as a reference system to perform MD simulations. The number of simulated λ states, which will be presented for each case in Sec. 2.5, must be sufficiently high to allow reweighting with appropriate precision. The optimization routine is carried out using the Powell algorithm 47 and Eq. (17) is integrated through Gaussian quadrature. We use the SciPy library 48 to apply both methods.

The S1 and S2 Schemes
Schemes S1 and S2, which rely on the smooth switching function defined in Eq. (7), differ only in the employed objective function. The former entails using L from Eq. (16) as the objective function, while the latter entails using var(∆G) from Eq. (17).
In both schemes, the reference system for carrying out the simulations was specified with param- 6, and λ f e = 1.0. During the optimizations, the coupling of the capped potential was forced to start at λ = 0 and that of the electrostatic potential was forced to finish at λ = 1, meaning that λ s c and λ f e are excluded from the set of optimization variables.
Two types of optimizations were performed. The first one involves four optimization variables, namely λ f c , λ s r , λ f r , and λ s e , under the constraints λ s r ≤ λ f c and λ s e ≤ λ f r , so that each coupling stage can only finish if the subsequent one has already started. The second type is a special case in which the residual and electrostatic potentials are coupled together, meaning that λ f r = 1 and λ s e = λ s r . In this case, the only two optimization variables are λ f c and λ s r , under the constraint λ s r ≤ λ f c .

The Q Scheme
Optimization problems were applied independently to coupling the capped, residual, and electrostatic potentials. For each case, the objective function of Eq. (17) was employed with a quartic polynomial as switching function as defined in Eq. (10). Thus, optimal parameters A k , B k , and C k are searched for in each case. The reference system for carrying out the simulations was specified based on Ref. 25, with linear scaling for the residual and electrostatic potentials and parameters A c = 1.62, B c = −0.889, and C c = 0.0255 for the capped potential.

Octanol-Water Partition Coefficient
The octanol-water partition coefficient (P o/w ) is defined as the equilibrium concentration ratio of a neutral solute between solvents 1-octanol and water. It is a signature of the relative affinity of a solute for an organic medium compared to an aqueous phase, as 1-octanol can represent soils, living tissues, and microorganisms. 8 Thus, P o/w has been used mainly to study systems of interest to the pharmaceutical industry, as in evaluating drug transfer between biological media. 49,50 Its experimental measurement has well-defined protocols, making the results very reliable. 51 In the common logarithmic representation, the octanol-water partition coefficient is related to the hydration free energy of the solute (∆G hyd ) and its solvation free energy in 1-octanol (∆G oct ), both at temperature T and atmospheric pressure, by 7,8 log Bannan et al. 7 organized a database of simulation results for 36 small organic molecules, as well as a comparison with experimental data found in the literature. [52][53][54] This database will serve as a reference for the studies carried out here. Like them, we used the hypothesis of complete immiscibility between the solvents, although experimental data indicate non-negligible miscibility of water in 1-octanol. 55 For each solvent, we consider a single solute molecule in the simulation box with periodic boundary conditions. This is reasonable, assuming that the experimental data to be compared were obtained at low solute concentrations.

Simulation Details
In the present study, we employed the highly efficient and customizable OpenMM 56 package (version 7.5.0) for all MD simulations. The optimizations were performed with three different solute molecules: phenol, ibuprofen, and cyanuric acid. Initial configurations, interaction parameters, and atomic charges were taken from the Free-Solv database, 57 which also provides experimental and calculated hydration-free energy values. These systems are modeled with GAFF force field 58 parameters and AM1-BCC 59 charges for the solutes, as well as TIP3P 43 parameters and charges for water.
We used the proposed LBF model with different switching functions and optimal parameters, as well as the soft-core model implemented in Ref. 31, to calculate octanol-water partition coefficients. In this case, we built initial configurations via Packmol 60 and determined interaction parameters via LEaP 61 by using the GAFF 58 and TIP3P 43 force fields again. With Antechamber, 61 we applied the AM1-BCC method 59 to determine atomic charges for 1-octanol as well as for the solutes vanillin, pentachlorophenol, nicotine, and trimethylamine.
In all alchemical protocols involving a single set of simulations to calculate capped, residual, and electrostatic potentials, we simulated states determined by 21 evenly spaced values of λ from 0 to 1 (∆λ = 0.05). In the protocols entailing separate simulations, we used 5 evenly spaced values for the residual and electrostatic potentials (∆λ = 0.2), while for the capped and soft-core potentials, we used 15 unevenly spaced values (∆λ = 0.1 for 0.4 ≤ λ ≤ 0.9 or ∆λ = 0.05 otherwise).
The simulations were executed under the isothermal-isobaric ensemble (NPT) at 298.15 K and 1 atm. Pressure control was done by a Monte Carlo barostat 62 with attempts to change the periodic box volume at every 50 fs, and the maximum allowable volume change was adjusted throughout the simulation to obtain an acceptance rate of around 50%. Temperature control was done via Langevin dynamics with a leapfrog, middletype integration scheme (LFMiddle) 63 and with a friction coefficient of 10 ps −1 . Each simulation consisted of 4 ns of equilibration, followed by 20 ns of production for an aqueous medium or 40 ns for a 1-octanol medium. A 2 fs time step was employed with samples collected at every 200 fs. The nonbonded cutoff distance was 1.2 nm and a switching distance of 1.1 nm was employed. The simulation boxes had initial edge lengths between 3.5 nm and 4.1 nm. Furthermore, we constrained the length of every chemical bond involving a hydrogen atom and considered rigid water molecules in accordance with the TIP3P model.
When using the PME for the long-range electrostatic interactions, the Ewald damping parameter was assigned as α = − ln (2δ)/r c , with an error tolerance δ = 5 × 10 −4 . Moreover, the number of nodes in the mesh was determined as n mesh = (2/3)αL box δ −1/5 , where L box is the length of the box along each dimension.
Coupling free energies and average properties were estimated using the MBAR method implemented in the PYMBAR library. 27,64 For the decorrelation of the sampled properties, we used the Overlapping Batch-Mean (OBM) estimator, 65 with blocks of size b = √ n , where n is the original sample size. 66 Note that we used this method instead of the integrated autocorrelation function method implemented in PYMBAR. 64 Uncertainties in free energies were computed using the MBAR estimator, 27,29 their propagation was handled via the delta method, 67 and their values are reported in the form of 95% confidence intervals. Table 1 shows results for solutes of diverse chemical classes, whose simulated hydration free energies reported in the FreeSolv database 57 (∆G fs hyd ) span a wide range, from −9 to −93 kJ/mol, as one can see in the second column of the table.

Modeling of Solute-Solvent Electrostatic Interactions
As described in Sec. 2.2, for each solute we carried out separate simulations with the electrostatic solute-solvent interaction forces derived from either the CRF (U as crf ) or the PME potential (U as pme ). The former simulation is computationally cheaper than the latter. During these simulations, the values of both U as crf and U as pme were sampled at regular intervals. Then, in a post-processing stage, the free energy difference G as pme − G as crf was estimated using both the BAR and the FEP methods. Such difference is the correction required for estimating standard PME-based hydration free energies from values obtained using the CRF model to evaluate solute-solvent interactions instead. The results are shown in the two rightmost columns of Table 1.
All studied substances whose simulated hydration-free energies lie above −46 kJ/mol require minimal CRF-to-PME corrections (< 0.4 kJ/mol). It is important to remark that more than 95% of the compounds in the current FreeSolv database (611 out of 642) have ∆G fs hyd > −46 kJ/mol. For cyanuric acid, the compound with the most negative value of ∆G fs hyd , the required correction is only 1.1 kJ/mol. Moreover, the correction itself can be accurately obtained from unilateral FEP calculations since the values computed in this way are agreed with those from the bilateral BAR method. Therefore, we only need to run the cheaper simulation in which the CRF forces are computed at every step, but U as crf and U as pme are computed only from time to time for sampling. This analysis shows that the reaction field model with infinite dielectric constant can satisfactorily substitute PME calculations when we need to compute solute-solvent interactions separately.
We also compared the CRF model with other pairwise alternatives. Eq. (14) is the limit when rf → ∞ of a more general Reaction Field (RF) where rf is the dielectric constant of the medium outside the cutoff cavity. 36 In this model, a unitary dielectric constant (vacuum) is considered inside the cavity. We tested it with rf = 78.3, which is the dielectric constant of water at room temperature. Although the derivative ν rf (x) fails to vanish at x = 1 when rf is finite, 37 this is not a problem here due to the presence of a switching function in Eq. (13). We also tested the Shifted-Force (SF) version of the Coulomb potential, which can successfully replace lattice-sum methods in numerous situations. 38 It consists in replacing ν crf (x) in Eq. (13) by ν sf (x) = 1/x + x − 2, which satisfies ν sf (1) = ν sf (1) = 0 as well. In fact, this is a special case of the Damped-Shifted-Force (DSF) model 38 ν dsf (x) = erfc(αr c x)/x − erfc(αr c ) + (x − 1)[erfc(αr c ) + 2αr c erfc(−α 2 r 2 c )/ √ π], which we also tested here with a damping parameter α = 0.2 A −1 .
In columns 3, 4, and 5 of Table 1, we can see BAR-estimated values of the corrections needed to obtain PME-based solvation free energies from simulations applying the DSF, SF, and RF models, respectively, for computing the solute-solvent interactions. Only the SF model exhibits performance comparable to that of the CRF model. The needed corrections for the remaining ones are systematically much larger, making them unsuitable for this type of calculation. Table 1 suggests that the CRF model is still slightly better than the SF model, but a more detailed comparison is needed to attest to its superiority. Our analysis consists in plotting coordinate pairs (U as pme , U as pme − U as crf ) and (U as pme , U as pme −U as sf ) sampled during a simulation employing forces derived from U as pme . If these pairwise models are suitable substitutes for PME calcula-tions, those coordinate pairs are expected to lie close to the horizontal lines U as pme − U as crf = 0 and U as pme −U as sf = 0, respectively. These plots are shown in Fig. 4 for three selected solutes, namely phenol, ibuprofen, and cyanuric acid, covering a wide range of solvation free energies. The color of each plotted pair is picked from a dark-to-light palette according to the concentration of other pairs around it. This makes the plots resemble two-dimensional histograms. In all cases, the CRF model exhibits a high probability that U as pme − U as crf ≈ 0, with symmetric distribution around U as crf = U as pme . For the SF model, however, we observe a marked tendency of U as sf to overestimate the value of U as pme , which makes it a less reliable substitute than the CRF model.
We used the CRF model, and the subsequent correction to the PME model, in all simulations involving solvation free energy calculations in this work, even for 1-octanol solvent instead of water. This includes the simulations with both the LBF and the soft-core approaches to validate our new proposal.

Optimization Analysis
This section presents and discusses the results of the optimization procedures described in Sec. 2.3. They are applied independently for the solutes phenol, ibuprofen, and cyanuric acid in an aqueous phase. Nonetheless, we intend to find sets of switching function parameters that perform well in all three cases. We expect these consensus parameters to provide reasonable alchemical pathways for a wide range of systems besides the analyzed ones. As shown in Sec. 3.3, the same parameters are suitable for calculations having 1-octanol as solvent.

Smooth Switching Function
We present the optimization results obtained with the S1 and S2 schemes described in Sec. 2.3. Both depend on the same MD simulations performed with the reference switching functions whose parameters were set as the initial guess for each scheme. In the first column of Fig. 5, one can see the coupling free energy profile obtained in these simulations, as well as the mean derivative ∂U/∂λ and the free energy variance as functions of λ.
For all three solutes, virtually constant free energies are observed between λ = 0.5 and λ = 0.75, creating a local maximum close to zero in the  Figure 4: Comparison between solute-solvent electrostatic potential energies computed for the same configurations by the PME method and a pairwise alternative: the Conductor-like Reaction Field model (top panel) or the Shifted-Force potential (bottom panel). The configurations were sampled during MD simulations based on the PME method to compute forces. The warmer the color tone corresponding to a potential pair, the greater the probability of finding that pair during the simulation.
∂U/∂λ curve. In contrast, a sharp decay in ∆G, followed by a slower decay and final stabilization, occurs at 0.75 < λ ≤ 1, leading to a deep well in the ∂U/∂λ curve. This behavior indicates that the reference alchemical pathway is inefficient because too many intermediate states are allocated to a λ range with little variation in free energy, and too few states are allocated to a range with high variation. To overcome this issue, minimizing sudden changes in the derivative ∂U/∂λ , as done in Scheme S1, has been a suitable approach. The dashed lines in the second column of Fig. 5 result from the optimal switching functions ob-tained by applying Scheme S1. As done in every optimization step, including the one responsible for finding the optimal solution, these curves were inferred via reweighting. As desired, they show a more uniform free energy profile and a ∂U/∂λ curve without unneeded local maxima and narrow wells. We have checked the effectiveness of reweighting by performing new simulations with the optimal parameters. Their results, which we omitted to avoid redundancy, perfectly agreed with the reweighted curves.
Returning to the first column of Fig. 5, we now analyze the free energy variance for the reference Figure 5: Property profiles for the hydration of multiple solutes obtained by a concerted pathway with smooth switching functions. Circles and solid lines indicate, respectively, simulation and reweighting results for the reference system. The dashed lines represent reweighting results at states defined by the optimal switching functions obtained through S1 and S2 schemes.
system. Considerably higher variance values occur between λ = 0 and λ = 0.25, where a free energy barrier is present. Thus, another way to improve the alchemical pathway is by reducing the variance in this region. We expect this to occur when the area under the variance curve is minimized in Scheme S2.
The third column of Fig. 5 contains curves obtained by reweighting to a state defined by the optimal switching functions found with Scheme S2. Again, we validated these curves through new simulations. The free energy barrier between λ = 0 and λ = 0.25 exhibits a more gradual variation, and the area under the variance curve diminished considerably. The gentler rise of the free energy curve toward its maximum relates to the shrinking of the hump in a mean derivative curve, whose peak is aligned to that in the variance curve.
The optimization involving a space with fewer dimensions was preferable for both schemes due to the similarity of its results compared to the case with more parameters. Thus, all analyses in this section refer to the optimizations in which λ s r = λ s e . The optimal parameters are presented in Table 2, in which f init obj and f opt obj refer to objective function values for the reference system and the optimal parameters, respectively. Table 2: Optimal parameters, as well as the initial and minimum objective function values, obtained through Schemes S1 and S2 for the concerted alchemical pathway with smooth switching functions. The objective functions of these optimization schemes are measured in kJ/mol and (kJ/mol) 2 , respectively.  Table 2 are indicated as red dots. The darker blue region contains the points considered near-optimal, whose f obj values are not larger than 10% above the minimum. This region has a considerable extension for all cases, leading to the possibility of a wide family of efficient alchemical pathways.
The optimal pairs (λ f c , λ s r ) obtained through Scheme S1 for phenol and ibuprofen were close to each other but considerably distant from that obtained for cyanuric acid. On the other hand, Scheme S2 led to very similar optimal pairs for all three solutes. Notice that such a geometric center of these pairs lies at (1.0, 0.35) and would belong to the dark blue region in all six plots of Fig. 6. For this reason, we chose the values λ f c = 1.0 and λ s r = λ s e = 0.35 as the consensus parameters for the smooth switching functions used in our concerted LBF formulation.
We expect the consensus parameters to provide efficient alchemical pathways for solvation-free energy calculations with diverse solutes and solvents. If this is not true for a specific case, using one of the optimization procedures discussed here can provide parameters for the desired system. Although the consensus parameters ultimately stemmed from Scheme S2, it is worth noting that the more straightforward objective function of Scheme S1 was also effective in finding efficient parameter sets.

Quartic Switching Function
We performed independent simulations for the capped, residual, and electrostatic alchemical potentials to apply the Scheme Q optimization procedure, which involves the quartic switching function. Fig. 7 exhibits results for the capped potential for phenol. We used this solute as an example, but the others behave similarly, as shown in analogous figures made available in the Supporting Information. In the region between λ = 0 and λ = 0.25, the free energy increases more gradually with the optimal switching functions than with those of the reference system. In the region where λ > 0.3, we observe an increase in the local variance, which was not as significant as the reduction in the preceding region. Table 3 shows the optimal parameters and the variance values before (f init obj ) and after (f opt obj ) optimization. In general, the residual and electrostatic potentials have a small influence on the global variance of the alchemical pathway. Thus, following Naden et al., 25,26 we recommend the use of linear switching functions in these cases for the sake of simplicity. For the capped potential, we will adopt as a consensus the optimal parameters obtained for ibuprofen, A c = 2.568, B c = −2.104, and C c = 0.419, which shows good agreement with those of the other solutes, as noted in Fig. 8. When contrasted to phenol and cyanuric acid, ibuprofen is a larger molecule with an intermediate hydrationfree energy value.

Octanol-Water Partition Coefficient
We used Eq. (19) to calculate octanol-water partition coefficients from solvation free energies determined in both water and 1-octanol media for the solutes trimethylamine (TMA), nicotine, pentachlorophenol (PCP), and vanillin. We obtained such free energy values through two different alchemical pathways, formulated with the switching functions and their consensus parameters discussed in Sec. 3.2. We have also used the ubiquitous softcore model for the LJ potential to validate the newly proposed methodology. Table 4 shows the solvation free energies in water and 1-octanol for the four analyzed solutes and each different alchemical pathway. We can see a good agreement among the alchemical pathways, with uncertainties of similar magnitudes. The table also displays simulation results by Bannan et al. 7 . The observed deviations are probably due to differences in details of the adopted MD procedures, but mainly in the atomic charges used for the solutes and 1-octanol. We note that the authors deliberately used OpenEye's QUACPAC application 68 instead of Antechamber 61 to apply the AM1-BCC method. Since the partial charges of the water model do not come from this method, the influence of charge dissimilarity would be better quantified if we could compare the LJ contributions to the hydration-free energies. However, this specific information is not available in Ref. 7. We remark that although the optimal quartic switching functions yield partial coupling free energies with low variances, as seen in Table 3, error propagation leads to overall free-energy uncertainties only slightly smaller than those resulting from the optimal concerted coupling. Such achievement of Scheme S is outstanding, considering that it involves only two optimization parameters, while Scheme Q involves nine parameters for a full solutesolvent coupling. Table 5 presents more details for comparing the Figure 6: Objective function contour plots for different solutes in water. The first and second rows refer to optimization schemes S1 and S2, respectively. Dark blue regions contain objective function values smaller than 10% above the minimum. Red dots correspond to the optimal parameter values.    concerted LBF approach using smooth switching functions to the sequential one using quartic functions. It contains the capped (∆G c ), residual (∆G r ), and electrostatic (∆G e ) contributions to the solvation free energies. In the case of the concerted approach, these contributions were com-puted by reweighting the simulation data to the milestone states of the sequential pathway, which are those defined by (h c , h r , h e ) equal to (1,0,0) and (1,1,0). We can see a good agreement between the free energy components for both LBF approaches. The uncertainty is systematically smaller in the sequential case for the residual term, which shows that the number of states used for this calculation was oversized. The Supporting Information also discriminates the LJ and electrostatic contributions to the solvation free energies obtained from soft-core-based simulations. A single electrostatic contribution calculation applies for both the sequential and soft-core pathways as a linear coupling has been performed for these cases. The soft-core strategy refers only to the coupling of the LJ potential. Fig. 9 allows an in-depth analysis of the concerted alchemical pathway for the different solutes in 1-octanol. The free energy presents smooth variations, the ∂U/∂λ curve is devoid of narrow valleys, and the variance is in a proper order of magnitude. The same conclusions are valid for the water medium results shown in the Supporting Information.
We can observe minor kinks in the meanderivative curves at λ = 0.35, where the coupling of the residual and electrostatic potentials starts (i.e., λ = λ s r = λ s e ). Hence, the second-order derivative of the free energy profile is discontinuous at such a point. This is not surprising because the employed switching functions, defined in Eq. (7) and (8), are only second-order differentiable at both λ s k and λ f k . Notwithstanding, the free energy profile is continuous everywhere, and the presence of kinks in its derivative does not influence the free energy calculation via independent simulations of intermediate states. However, in the context of λ dynamic simulations, where the free energy derivative is the mean force exerted on λ, it might be advisable to use an even smoother family of sigmoid-shaped switching functions.
In the Supporting Information, we can find figures aimed at studying simulations performed using the sequential LBF approach with quartic switching functions for the four solutes. They show the negligible contribution of the residual potential to the total variance of the alchemical transformation concerning the other potentials. It is interesting to note that, despite this, the residual potential contributes significantly to the free energy. We can also verify that the free energy variances of water  as solvent were smaller than those with 1-octanol instead. In the latter case, the electrostatic potential had a significant weight in the total variance, which was not observed with the water system. Table 6 shows the logarithms of partition coefficients calculated from the free energy data pre-sented in Table 4. Experimental data from the literature 52,54 are also provided for comparison. All values are positive, which denotes that all solutes have more affinity for 1-octanol than for water. For each chemical species, we can observe good agreement among the results of simulations using the soft-core model and the LBF schemes but significant deviations when comparing them to the experimental data. For example, the deviations are around 120% for nicotine and 37% for vanillin. This issue is not due to the calculation methods but to the inadequacy of the employed force fields.

Conclusions
We have proposed a new splitting of the Lennard-Jones (LJ) potential into capped and residual terms and employed the Linear Basis Function (LBF) formalism to perform alchemical free energy calculations. This splitting is more straightforward than previous proposals once it involves coefficients independent of the LJ parameters. Additionally, we applied a set of smooth switching functions to control the coupling of the multiple potential energy terms, allowing a concerted evaluation of LJ and Coulomb components of the solvation-free energy. The capped potential avoids numerical issues by replacing the native LJ potential at the onset of the alchemical transformation. The coupling of the residual and electrostatic potentials starts at a later stage when the risk of numerical instabilities has vanished.
The inherent reweighting-related features of the LBF model allowed us to carry out fast optimization procedures. As a result, we could find switching function parameters that led to safe and ef- Table 6: Logarithm of the octanol-water partition coefficient for different solutes using the free energy data previously obtained in ficient alchemical pathways for solutes spanning a wide range of molecular sizes and hydrationfree energies. Given the similarity of the optimal parameters, we determined consensus values that should be efficient in general solvation free energy calculations. Using our concerted LBF approach with consensus parameters, we performed solvation free energy simulations, both in water and 1-octanol media, of solutes other than those used for optimization. Their results allowed us to calculate the octanolwater partition coefficients and compare the obtained values with literature data. We also employed the sequential LBF pathway and the ubiquitous soft-core model to perform analogous simulations. Through a successful comparison with these established methodologies, we attested to the correctness and robustness of the concerted LBF pathway.
In addition, we observed that a reaction-field model with infinite dielectric constant was the best among several pairwise options for substituting lattice-sum methods in solute-solvent electrostatic calculations.  Free energy values in kJ/mol for the coupling of Lennard-Jones (∆G lj ) and electrostatic (∆G e ) potentials in the application of the concerted LBF, sequential LBF, and soft-core alchemical pathways. Note: a single electrostatic contribution calculation applies for both the sequential and soft-core pathways.

Solute
Pathway Water as solvent 1-Octanol as solvent ∆G lj ∆G e ∆G lj ∆G e trimethylamine Figure SI.2: Property profiles regarding the coupling of the residual-LJ potential with a quartic switching function in the hydration of phenol, ibuprofen, and cyanuric acid. Circles and solid lines respectively indicate simulation and reweighting results for the reference system. Dashed lines represent reweighting results in the case of the optimal switching function derived from Scheme Q.