Hello,<div><br><div> I'm attempting to calculate a relative free energy of hydration between Argon and Neon. However when perturbing from Argon(aq)->Neon(aq) I was surprised to compute a change in free energy of 1.2kJ/mol, where the literature values suggest that the hydration energy of Argon is ~-9kJ/mol and neon is ~-8kJ/mol, which means that I should have seen a change in free energy of ~-1kJ/mol.</div>
<div><br></div><div> Interestingly, I noticed that this 2kJ/mol difference could be resolved if I subtract the dEkin/dlambda value obtained from the edr file. A clumsy look into the source code seemed to confirm that there is a dEkin/dlambda term included in the dH/dlambda calculation which results from the change in mass during the mutation. However it seems in the literature that most people calculate the free energy between two states as the difference of the configurational partition functions, and thus should be a function of the positions and independent of the mass. </div>
<div><div><br></div><div>It seems tempting to me to resolve my problem by simply subtracting off the dEkin/dlambda term, since I am in the NPT ensemble (and at constant temperature the Ekin shouldn't change). Does this seem reasonable? Is there a good reason for including this dEkin/dlambda term in FEP calculations? If I'm interested only in the difference between the configurational partition functions in the NPT ensemble, should the dEkin/dlambda term be included in the first place?</div>
</div><div><br></div><div><br></div><div><br></div><div>In case it helps, my topology file and mdp file:</div><div><br></div><div><b>Topology File</b></div><div><div>; Argon parameters from Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A.</div>
<div>;Intermolecular Forces: Their Origin and Determination;</div><div><br></div><div>[ defaults ]</div><div>; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ</div><div>1 2 yes 1.0 1.0</div>
<div><br></div><div>[ atomtypes ]</div><div> Ar Ar 39.948 0.000 A 3.41000e-01 9.96000e-01</div><div> Ne Ne 20.180 0.000 A 2.72000e-01 3.91000e-01</div><div> DNe DNe 20.180 0.000 A 0.00000e-01 0.00000e-01</div>
<div> DAr DAr 39.948 0.000 A 0.00000e-01 0.00000e-01</div><div><br></div><div>; Include water topology</div><div>[ atomtypes ]</div><div>OW OW 16.00 0.0000 A 3.15061e-01 6.36386e-01</div>
<div>HW HW 1.008 0.0000 A 0.00000e+00 0.00000e+00</div><div>[ moleculetype ]</div><div>; molname nrexcl</div><div>SOL 1</div><div><br></div><div>[ atoms ]</div><div>; id at type res nr res name at name cg nr charge mass</div>
<div> 1 OW 1 SOL OW 1 -0.834 16.00000</div><div> 2 HW 1 SOL HW1 1 0.417 1.00800</div><div> 3 HW 1 SOL HW2 1 0.417 1.00800</div>
<div><br></div><div>[ settles ]</div><div>; i j funct length</div><div>1 1 0.09572 0.15139</div><div><br></div><div>[ exclusions ]</div><div>1 2 3</div><div>2 1 3</div><div>
3 1 2</div><div><br></div><div><br></div><div>[ moleculetype ]</div><div>; Name nrexcl</div><div> Argon 1</div><div><br></div><div>[ atoms ]</div><div>; residue 3 LEU rtp CLEU q -1.0</div>
<div>; nr type resnr residue atom cgnr charge mass typeB chargeB massB</div><div> 1 Ar 1 Ar Ar 1 0.0000 39.948 Ne 0.0000 20.180</div><div><br></div><div>[ system ]</div>
<div>; Name</div><div>Argon Mutate Neon</div><div><br></div><div>[ molecules ]</div><div>; Compound #mols</div><div>Argon 1</div><div>SOL 215</div></div><div><br></div><div><b>MDP File</b></div>
<div><div>; Run control</div><div>integrator = sd ; Langevin dynamics</div><div>dt = 0.002</div><div>nsteps = 2500000 ; 5 ns</div><div><br></div><div>; Neighborsearching and short-range nonbonded interactions</div>
<div>nstlist = 10</div><div>ns_type = grid</div><div>pbc = xyz</div><div>rlist = 0.9</div><div>; Electrostatics</div><div>coulombtype = PME</div>
<div>rcoulomb = 0.9</div><div><br></div><div>; van der Waals</div><div>vdw-type = switch</div><div>rvdw-switch = 0.7</div><div>rvdw = 0.8</div><div><br></div>
<div>; Apply long range dispersion corrections for Energy and Pressure</div><div>DispCorr = EnerPres</div><div><br></div><div>; Spacing for the PME/PPPM FFT grid</div><div>fourierspacing = 0.1</div>
<div>; EWALD/PME/PPPM parameters</div><div>pme_order = 6</div><div>ewald_rtol = 1e-06</div><div>epsilon_surface = 0</div><div>optimize_fft = no</div><div><br></div><div>; Temperature coupling</div>
<div>; tcoupl is implicitly handled by the sd integrator</div><div>tc_grps = system</div><div>tau_t = 1.0</div><div>ref_t = 298</div><div><br></div><div>; Pressure coupling is on for NPT</div>
<div>Pcoupl = Parrinello-Rahman</div><div>tau_p = 0.5</div><div>compressibility = 4.5e-05</div><div>ref_p = 1.0</div><div><br></div><div>; Free energy control stuff</div>
<div>free_energy = yes</div><div>init_lambda = 0.0</div><div>delta_lambda = 0</div><div>foreign_lambda = 0.05</div><div>sc-alpha = 0.5</div><div>sc-power = 1.0</div>
<div>sc-sigma = 0.3</div><div>; options for bonds</div><div>constraints = h-bonds ; we only have C-H bonds here</div><div>; Type of constraint algorithm</div><div>constraint-algorithm = lincs</div>
<div>lincs-order = 12</div></div><div><br></div><div>I used 10 windows for the perturbation and ran each window for 3ns. The results were then combined using the g_bar utility</div>
<div><br></div><div>Thanks!</div><div>-Connie</div>
</div>