[gmx-users] Effect of bond constraints on equilibrium densities of equilibrated liquids

[Miguel Caro]

Hi all,

I am running some test calculations with acetonitrile/water mixtures
using the GAFF topologies as available from virtualchemistry.org (by the
way, thanks for this great resource). I obtain equilibrium densities
from averages over different independent configurations after 1 ns
equilibration periods (298 K, 1 bar, 400 water molecules and 200
acetonitrile molecules at the end points). The equilibrium density I
obtain for pure acetonitrile is larger than that reported on the website
for GAFF
(http://virtualchemistry.org/molecule.php?filename=acetonitrile.zmat)
and slightly closer to the experimental value of 0.78 g/cm^3: it is
about 0.8 g/cm^3 for me compared to 0.73 g/cm^3 from the reference
calculation. I hope the graph copy-pasted below shows up on my email, to
help illustrate the issue (each data point is one configuration
equilibrated for 1ns). Going to the JCTC paper referenced on the website
(Caleman et al. JCTC 8, 61 [2012]), it appears unless I'm mistaken that
the reference results were obtained with constrained bonds (Table 1 on
that paper, "LIQ" item). However, I am allowing all bonds to vibrate in
my simulation (with appropriate small time steps).

My question is whether I can/should expect the difference in predicted
density to originate from the lifted bond constraints, or is there
something else flawed about how I'm running this simulation. I.e.,
should I run longer simulations, calculate for each configuration the
average density over the simulated trajectory instead of taking the
final density value, etc.?

Many thanks for you help,
Miguel


-- 
*Dr. Miguel Caro*
/Postdoctoral researcher/
Department of Electrical Engineering and Automation,
and COMP Centre of Excellence in Computational Nanoscience
Aalto University, Finland
Personal email: *mcaroba at gmail.com*
Work: *miguel.caro at aalto.fi*
Website: http://mcaroba.dyndns.org