Ammm:chgcorr

Corrections in free energy calculations when the simulation box/ligand has a net charge

Pengyu Ren

Hydration free energy

HFE is the free energy of solute transfer from gas to solution. Experimentally only salt property can be determined reliably and single ion HFE requires additional procedure to decompose. Thus there has been a great deal of controversy on what is experimental HFE (see discussion in (1)).  Experimental hydration free energy may contain two corrections that need to be removed before compared to PBC simulation results using Ewald:

1)     Standard state correction. The standard state of gas phase should be ideal gas (24.46 L/mol from PV=nRT). In simulation, we decouple the solute from water to “vacuum” but still in liquid concentration. This correction is the free energy needed to expand from 1 L/M to 24.46 LM. dG=RTln(24.46/1)=1.9 kcal/mol. Certain Expt HFE contains this correction already (e.g. using 1 M as standard). If not the expt HFE should be correction by -1.9 kcal/mol (more negative). Or you can add +1.9 to simulated HFE (make it less negative) to compare with expt if its standard state is ideal gas.

2)     Net charge related correction. There are multiple corrections related to the net charge of the solute in a PBC box.

a.      The first (important) is related to surface potential (Galvani potential) between gas-liquid interface (water at the interface has specific orientation). In decoupling simulation our solute never really crosses an interface. The surface potential has been calculated using force fields by integrating charge density from gas to solution (z direction). See Dang(2) (eq 2-4) and Beck paper (3) for how to calculate surface potential. The equation in Dang’s paper (based on field) is more useful for AMOEBA polarizable force field since we have both dipole and quadrupoles instead of just atomic charges; Or you can compute using ZJ’s method(4) . For fixed charge FF SPC/E and TIP3P etc, the SP has been reported. For SPC value is -0.55 or -12.7 kcal/mol-e (Tom’s paper (5)), For SPC/E -0.6 to -0.65V or -13.8 to -14.9 kcal/mol-e (3). For AMOEBA, ZJ calculated it to be -0.352 V or -8.1 kcal/mol.(4) Simulated HFE would be more negative when added this SP contribution.

b.      H+ HFE has been used as a reference. The expt value was reported to be -251.4 kcal/mol using 1 M as standard state and without SP. With SP (MAOEBA value above), H+ HFE would be -260.(4). Note that for ion HFE from Marcus 1991 paper, Marcus used the “wrong” reference so his value needs to shift by n*2.9 kcal/mol (n is charge). Marcus value also corrected for the SS but in wrong direction. So overall Marcus value needs two corrections (e.g. La3+ HFE= -3145/4.18 -2*1.84+3*2.9 =-747 kcal/mo where -3145 was from Marcus 1991 paper; the 1.84 is SS correction and 2.9 is H+ reference).

c.      I Beck’s paper, another positive local potential (LP) was mentioned. This refers to the LP due to nearby water even when the solute is a vdw sphere/neutral particle. This was calculated to be 8-9 kcal/mol. I don’t think we need to consider this as in decoupling simulations local water is sampled explicitly.

d.      Tom’s paper (5) mentioned additional correction, potential shift due to PME neutralizing plasma and polarization artifacts in PBC (infinite lattice). These terms disappears when box L is large and you can check by compute HFE for small vs large box. 45-50 Ang box size is sufficiently large for Mg2+ (see Tom, ZJ and Brandon’s paper).

e.      Side note, TLTR: Tom also mentioned a charge sorting/summation order related quantity that is 7-8 kcal/mol. Numerically this is the difference between PBC PME and water cluster simulations. I believe this is a combination of SP above and the PME boundary effect. PME uses a tinfoil conducting boundary. The boundary correction physically is related to how you treat boundary (conducting or no). This is different from correction due to conditional convergence of Ewald (which is called self-correction in AMOEBA code). For atom-charge only FF, when standard tinfoil boundary is used, the PME energy is the same as using atom based infinite cutoff (huge cutoff without Ewald). But you could also use group or molecule-based cutoff, and the energy will converge to a different value as you increase the cutoff. This energy can however also be achieved by PME but adding a boundary correction term that is related to the total dipole moment of the box. We have this boundary correction term in the AMOEBA code. AMOEBA however is complicated because we use atom-based cutoff for charges but dipole etc is automatically “group” based.

In summary, when computing ion HFE and comparing to expt data, you need consider two corrections 1) standard state and 2) surface potential 3) and use large box (45 or above). The best way to calibrate “expt” is to compare with the cited experimental values in ZJ paper (4) which uses 1 M as SS and no SP (so you can compare simulated HFE directly) . Note the corrections cancel between ions of same charges or in ion pairs (salt).

Binding free energy

The corrections above do not cancel between solvent and complex since protein dielectric is very different from water. The best practice is then to add a counter ion (dummy) to keep the system always neutral even when lambda =/= 1 (decouple the counter ion with your ligand together). If your ligand is +2 then counter ion has -2. This is only for ele part of free energy so vdw does not matter.  Select a counter ion that is far away from complex (you may use couple distance restraints to keep it away). The contribution of this counterion cancels between the solution and complex. This won’t work for absolute HFE above of course.

References

1.           Grossfield A, Ren P, Ponder JW. Ion solvation thermodynamics from simulation with a polarizable force field. J Am Chem Soc. 2003;125(50):15671-82. Epub 2003/12/11. doi: 10.1021/ja037005r. PubMed PMID: 14664617.

2.           Wick CD, Kuo IFW, Mundy CJ, Dang LX. The Effect of Polarizability for Understanding the Molecular Structure of Aqueous Interfaces. Journal of Chemical Theory and Computation. 2007;3(6):2002-10. doi: 10.1021/ct700098z.

3.           Beck TL. The influence of water interfacial potentials on ion hydration in bulk water and near interfaces. Chemical Physics Letters. 2013;561-562:1-13. doi: https://doi.org/10.1016/j.cplett.2013.01.008.

4.           Jing Z, Qi R, Liu C, Ren P. Study of interactions between metal ions and protein model compounds by energy decomposition analyses and the AMOEBA force field. J Chem Phys. 2017;147(16):161733. Epub 2017/11/04. doi: 10.1063/1.4985921. PubMed PMID: 29096462; PMCID: PMC5645194.

5.           Lin Y-L, Aleksandrov A, Simonson T, Roux B. An Overview of Electrostatic Free Energy Computations for Solutions and Proteins. Journal of Chemical Theory and Computation. 2014;10(7):2690-709. doi: 10.1021/ct500195p.