[gmx-developers] Velocity verlet version of Gromacs with trotter-decomposition NPT

[Berk Hess]

Hi,

Sorry for not answering earlier.
I have been very busy.
I hope to have time to look into the code in detail next week.

Thanks for all the effort!

Berk

Michael Shirts wrote:
> Hi, all-
>
> After much excruciating debugging, and 10 months or so after I thought
> it would be done, I have a velocity Verlet Trotter decomposition based
> implementation including both Nose-Hoover chains and isotropic
> pressure control, all of which work with constraints.  It was the
> pressure control + constraints that held back everything for most of
> that time, but it appears to be working perfectly now.   By working,
> it conserves the conserved quantity(NVE,NVT,NPH,NPT), with quadratic
> dependence on error in the timestep as well good as the NVT
> simulation.  This works best in double, of course -- in single, it
> conserves energy as well as current methods, though starts to have
> problems with constraints for small timesteps.
>
> It also preserves previous MD methods without any differences in the
> output caused by more than order of operations changes in the
> numerical precision (at least, all that I've checked -- there may be
> some combinations of options that still have bugs).  It has also been
> verified with checkpoints, and passes all the regression tests. It is
> synchronized with the main git distribution as of September 25, 2009.
> I have not played around with multiprocessor jobs yet, however.
> Rerunning does not quite work for the ekin and pressures in the
> trotter-decomposition NPT version because of the integration of the
> Nose-hoover chains; that might take some extra work to sort out, but
> potential energies are fine, as is rerunning without NPT.
>
> The code is in my public repository at:
> ssh://git.gromacs.org/~mrshirts/public_git/gromacs_mrs.git
>
> It would be very good to have Berk look over it and see how much can
> be ported over as is (hopefully most of it!)  There's a few code
> architecture decisions that he would likely have some better ideas on;
> and I'm certain there are a few bugs that will pop up.  I'm personally
> willing to do whatever needs to be done to get these features in and
> working now so that I don't have to spend a day or two synchronizing
> and debugging every time a big change is made to the code base!
>
> There are also a number of free energy features that I'd like to add
> in as soon as these integrator changes are incorporated.  These should
> be much simpler to add since I've got them going in my 3.1.4 version
> already and they are pretty portable --- unlike this batch of changes,
> which is pretty heavily woven throughout a number of different parts
> of the code.
>
> There are a number of additional features to implement that have not
> yet been added, but will be relatively straightforward to implement in
> a unified code framework instead of working on a code branch, and can
> be done in the order that people require them.   These include:
>
> -Documentation of the methods -- I've been working off of a Mark
> Tuckerman preprint; I did the iteration a bit differently than he did,
> and I have not entirely finished the .tex documentation files.  These
> can be done over a day or two as needed, certainly in time for 4.1.
>
> -Nonisotropic scaling: This will be pretty trivial to put in.  It's
> mostly a matter of replacing scalars with matrices that are simple to
> compute and set up.  Will take a couple of days, but is
> straightforward.
>
> -Semiisotropic scaling: This would be a bit harder:  I'd have to sit
> down and do the new derivations for the trotter decomposition.  After
> the exact equations are worked out, it can be implemented in the
> nonisotropic scaling framework, so that won't require any additional
> coding at that point.  Will come after nonisotropic scaling, might
> take some thought.
>
> -velocity verlet temperature control using  KE(t) = 0.5*[(KE(t+dt/2) +
> KE(t-dt/2)] instead of KE(t) = \sum ( v(t-dt/2) + v(t+dt/2))^2. Using
> average KE vs. average velocities gives worse fluctuations in the
> conserved quantity but an average temperature that depends much less
> on the timestep.  Without temperature control, it doesn't matter --
> you're getting the same trajectory whichever definition you use, but
> with temperature control it can mean that for large steps, you are not
> keeping the system at the temperature you think you are.  For small
> timesteps, this difference is not large, but if you're taking 4 fs
> timsteps, it could cause problems in terms of the temperature that you
> -think- your simulation is running at.  This is pretty easy to code in
> without constraints.  With constraints, however, it requires solving a
> self consistent equation, which can be done -- it just takes some
> extra effort.  Might take some more mucking around.  I'm not exactly
> certain how important this is in the near term, could definitely be
> considered in the longer term.
>
> -Nose-Hoover with chains with the leapfrog Nose-Hoover algorithm.
> Leapfrog uses a linearized version, whereas my implementation follows
> Martyna-Kline-Tuckerman and is not; so it would take a bit of work,
> but wouldn't be too bad.
>
> -Changing pressure / thermostat coupling to only every N steps, thus
> requiring virial and ekin calculations, by rewriting the order of the
> Trotter decomposition. Not yet implemented, but -should- be (unless
> I'm missing something) pretty simple to do.  This should mean greatly
> improved scaling with pressure and temperature control.  You would
> have to have very long tau_t and tau_p for this to work, even with
> something fairly low like N=10; but for long simulations, long tau_t
> and tau_p are fine in any case.
>
> -Converging the iterative equations; Right now these are controlled
> with some heuristic guidelines which may occasionally lead to lack of
> convergence -- these will likely need to be adjusted slightly.
>
> Some more detailed coding issues:
>
> -I had to significantly rearrange the code of the do_md loop and the
> update() function to get this to work.  I had to have iterative loops
> to solve some of the implicit equations for constraint + NPT, so I
> introduced a number of helping functions to do this.  I used a
> modified secant search instead of self-consistent iteration ion, which
> greatly improves the convergence -- in most cases, it converges in
> just 6-8 iterations for even short temperature and pressure time
> constants (for example, a tau_p=50 fs tau_t=50 fs combination).  For 1
> ps tau_t/tau_p, it takes 4 iterations pretty much every time; with
> some fiddling, this could probably be reduced to 3 or so.
>
> -Two major structural components were splitting update() in to several
> different components, since I need to iterate over several different
> components of the updating algorithm, and grouping all of the "gather
> information from multiple states" into a single compute_globals call,
> since this needed to be done multiple times per loop -- there's a
> bunch of messy flag passing that could probably be simplified here.
>
> -There's a number of smaller changes here and there -- these include
> including a few more variables in structures to avoid passing things
> back and forth as much, but these are the biggest ones.
>
> -The number of Nose-Hoover chains is hard coded for simplicity.  This
> should probably be switched eventually, just on principle.  Will
> require changing code in a number of places, but should be
> straightforward.
>
> Theoretical note:
>
> You can write either velocity verlet or leapfrog as a trotter
> decomposition; it's just a matter of where you start from.  However,
> in most cases leapfrog requires things before you calculate them, so
> doing pressure control is very difficult.  Velocity verlet requires
> them after you calculate them, which usually makes it cleaner.  I can
> explore various other combinations as time permits and people are
> interested.
>
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> (434)-243-1821
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.
>