Validate density and dielectric computations for TIP3P against Swope/Rice/Horn/Pitera "Little Water Code"

[jchodera]

I now have the "Little Water Code" running on the cluster and can generate data for TIP3P (and other water models) for densities and dielectric constants. We should compare our computed values to the computed values from this code to make sure things are free of bugs.

[jchodera]

I've more or less figured out how to run it. What would be the most useful experiment for me to run here?

• Run a predefined set of systems (e.g. TIP3P) at a set of temperatures/pressures (e.g. 1 atm, and several temps) and write a paired analysis code so we can check by hand?
• Write a Python wrapper for the whole thing so you can plug it into the existing framework?

[kyleabeauchamp]

My goal is to switch to using amber files first.

1. Switch to using standard amber files
2. Benchmark against sander for something like EtOH
3. Then benchmark against LWC

My reason for this ordering is that we're using GAFF forcefields here, so I think the first sanity check we want is to compare against AMBER proper.

[kyleabeauchamp]

Also, the current benchmark isn't actually using TIP3P, so it's a bit unclear how to best run the LWC comparison.

[jchodera]

You'd have to set up a special TIP3P prmtop/inpcrd pair for you to use to do that test, but I can do that pretty easily. This is to test the simulation/analysis leg of the calculation.

[jchodera]

Note that AMBER has made some choices regarding what algorithms are supported that may make quantitative comparison a bit tricky. It may be good to rope in Jason Swails if we need help there.

• Lennard-Jones uses hard cutoff without switch (though I think OpenMM may use this by default as well)
• Langevin is supported (ntt=3, use gamma_ln for collision rate in 1/ps, temp0 for temperature, ig=-1 to randomly seed random number generator)
• Pressure coupling is ntp=1 for isotropic scaling (which requires ntb=2) but only the Berendsen algorithm is supported, which does not generate the correct distribution. taup is the coupling time (in ps) with default of 1 ps.
• PME selects grid parameters (nfft1, nfft2, nfft3) by default. The defaults for order (4), dsum_tol (1e-5) and rsum_tol (1e-5) may not have direct equivalents in OpenMM. I recall something about an OpenMM feature request to be able to retrieve these values from an OpenMM Context for easy comparison, but I'm not sure if that was implemented.

[kyleabeauchamp]

Have you been able to get stable long trajectories with LWC? I'm finding that trajectories seem to be stable on the short timescales used in all the included examples, but unstable when I extend things to be longer.

I've looked at both the econs_3p/run_econs.sh and boxes/tip3p_1000.lwc with boxes/lwc.inp. The most telling sign is volumes (fort.8) seem to grow linearly after about ~1ps or so.

In the example below, I've run the four simulations in run_econs.sh by a factor of 100X more steps than were originally present. The volume (y axis) goes nuts after about 200 steps (x axis).

[kyleabeauchamp]

Also note: to make lwc.inp work, you will have to disable the RDF calculation, which runs into memory issues for longer trajectories.

[jchodera]

Here's what I have, which I believe is a 1 ns NPT simulation that prints out energies, box dipoles, and box sizes:

cd /cbio/jclab/home/chodera/lwc/lwc/experiments/tip3p/npt-dielectric/production
grep "Etot=" lwc.out | cut -c 20-30 > Etot.out

[jchodera]

Here's a plot from using a 2 fs timestep.

Here are the parameters I believe I used:

! main controls
NSTEPS    =       50000000 ! number of MD timesteps (I)
DT        =   0.0020000000 ! timestep size, picosec (F)
ISEED     =              1 ! random number seed (I)
DENS_INI  =            1.0 ! density to be used as initial density (NpT) or system density (NVT/E)

! pressure controlled ?
IFBOX     =              1 ! pressure control flag (1=p-control; 0=const V) (I)
PEXT      =   1.0000000000 ! External Pressure for pressure control, Atm (F)
BOXMAS    =        0.00040 ! piston mass, AMU/(Angstrom**4) (F)

! temperature controlled ?
TBATH     =          298.0 ! Bath Temperature for temperature control, Kelvin (F)
INEWVEL   =              1 ! velocity update on startup (I)
N_NEWVEL  =            500 ! velocity update period (1=everytime 0=never) (I)

! potential-related flags
IFSWITCH  =              1 ! flag for switched (1) vs. truncated (0) interactions (I)
RLOWER    =            9.0 ! lower cutoff/switch boundary in Angstrom (F)
RUPPER    =            9.5 ! upper cutoff/switch boundary in Angstrom (F)
IFCORR    =              1 ! continuum vdw correction flag (I)

! ewald parameters
IFEWALD   =              1 ! normal nonbonds (0) or Ewald summation (1) (I)
IEEXCL    =              1 ! calculate Ewald exclusions (I)
IEDIR     =              1 ! calculate Ewald direct-space term (I)
IEREC     =              1 ! calculate Ewald reciprocal-space term (I)
KMAX      =             12 ! maximal k-index in any direction (I)
KSQMAX    =            144 ! spherical cutoff**2 for k-vectors (I)
ALPHA     =           0.35 ! width of screening gaussian, in Angstrom (F)
DIELEC    =   0.0000000000 ! continuum dielectric for RF/Ewald dipole term (F)

! neighbor list ?
IFNLIST   =              1 ! neigborlist flag (1=build/use list; 0=dont) (I)
RLIST     =   1.0000000000 ! cutoff used to build verlet list (F)

! rigid or flexible ?
IFRIGID   =              1 ! shake/rattle flag (1=constraints; 0=vibrate) (I)
SHKTOL    =   0.0000000001 ! shake/rattle tolerance, Angstrom (F)

! force field flags
IWATMOD   =              3 ! Water model flag (4=TIP4P; 5=TIP5P; ...) (I)
SIGMA     =   3.1643500000 ! LJ Sigma(O-O) (4P: 3.15365) (ignored unless IWATMOD<0) (F)
EPSILON   =   0.1627500000 ! LJ Epsilon(O-O) (4P: 0.155) (ditto) (F)
CHARGE_M  =  -1.0484400000 ! Coulomb charge(M) (4P: -1.04) (ditto) (F)
DIST_OM   =   0.1250000000 ! Distance O-M(idpoint) (4P: 0.150) (ditto) (F)
DIST_OH   =   0.9572000000 ! Distance O-H (TIPnP: 0.9572) (ditto) (F)
ANGLE_HOH = 104.5200000000 ! Angle H-O-H (TIPnP: 104.52) (ditto) (F)

! sampling & output periods
N_WRT     =             50 ! output period for instantaneous values file (fort.8) (I)
N_TRJ     =              0 ! period for dumping Coords+Velocs to trajectory file (fort.17) (I)

! compute RDFs ?
N_RDF     =              0 ! period for sampling RDF (0=never) (I)

! dielectric ?
IFDIPMO   =             50 ! compute system dipole every IFDIPMOth timestep (I) (on fort.11; for dielectric)
EFIELDz   =   0.0000000000 ! finite field (E(z), in kV/cm) for numeric dielectric computation (F)

! kinetics flags
IFKIN     =              0 ! compute kinetics props (diffusion coeff, etc) every IFKINth timestep (I) (on fort.10)
RHBOND    =   4.0000000000 ! minimum O-O distance threshold for H bond, in Angstrom (F)
AHB       =  30.0000000000 ! minimum O-H-O angle threshold for H bond (max. dev. from linear), in degrees (I)
ICOMVEL   =              0 ! if 1, remove the velocity of the simulation box (implied by ifkin > 0)

[jchodera]

Volumes do look weirder:

[kyleabeauchamp]

Ok, maybe we should ask the authors? On Mar 19, 2015 7:33 PM, "John Chodera" notifications@github.com wrote:

Volumes do look weirder: [image: volumes] https://cloud.githubusercontent.com/assets/3656088/6743331/cee46186-ce6e-11e4-8663-e4700c59debc.png

— Reply to this email directly or view it on GitHub https://github.com/choderalab/LiquidBenchmark/issues/27#issuecomment-83805430 .

[jchodera]

Hans Horn is the man to ask! Will send you his email info.

[jchodera]

Email sent.

I have a longer simulation running now too.

[kyleabeauchamp]

At this point, I believe we're both still waiting on more detailed responses from the LWC folks. Hopefully this week they will be able to provide us with working example files that give stable NPT trajectories at long time scales

[kyleabeauchamp]

This is now finished, with all code checked in to the LWC analysis repo. Densities are identical to 0.3%, with dielectrics also comparable to within uncertainty (using same error approaches as in manuscript). Even better agreement might be possible, but we would have to exactly match PME parameters and run LWC for ~40 days or so.

# LWC
In [15]: density.mean(), density_sigma
Out[15]: (1.034951753687456, 0.00043367431911135002)

In [16]: static_dielectric, sigma_dielectric
Out[16]: (91.664784881404046, 3.4119206720802784)

# OpenMM
In [25]: mu, sigma_mu
Out[25]: (1.0381425181387824, 4.0022215019214846e-05)
In [26]: dielectric, dielectric_sigma
Out[26]: (91.1848339567202, 0.53559782308956261)