I am finding that use of the VV2 integrator and TPCONTROL makes molecular dynamics run about 50% slower than the same system and options simulated with the leapfrog integrator and Berendsen thermostat. Use of VV2 without TPCONTROL results in about a 30% reduction in speed. The parallel timings show that significant amounts of time are going towards integration when VV2 is used. Why is this? What part of VV2 integration is taking up so much time?
I am using a system of DNA, an attached dye, water and ions with about 19200 atoms with particle mesh Ewald.
The attached zip file contains 3 files: test.out is the system with VV2 and TPCONTROL, test2.out is leapfrog and Berendsen, and test3a.out is VV2 only.
With leapfrog, you can also use the extended system CPT code (pressure.doc) for NVT, NPT, and other ensembles.
There are some problems with the weak coupling scheme employed by the Berendsen thermostat, and the (archaic) CHARMM implementation only supports a single heat bath. I believe the GROM* programs allow multiple heat baths, which avoids a common problem of T differences between phases. With a single heat bath, protein/water or octane/water systems will have T differences of 30-40 deg between the phases, and only their average matches the bath T. I've been told those programs now include an extended system for T and P regulation, in addition to the original weak coupling methods.
Why do you need a thermostat? NVE should be perfectly fine for DNA (CHARMM is sufficiently bugfree to conserve energy in NVE). Use the lookup command to increase speed (nbonds.doc).
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
I have tried CPT; the performance is only a little worse than berendsen. (The Berendsen thermostat was just the first thing that came to mind when I was thinking of something else to try.) What I'm really interested in is why VV2/TPCONTROL is so inefficient.
I would like to be simulating the canonical ensemble so that I can use enhanced sampling techniques like adaptive umbrella sampling, which assume constant temperature.
The use of lookup provides a modest speedup, about 8-9%, but the hint is gratefully appreciated anyway.
I looked at your outputs - the real space nonbonded force caclculation is a relatively small part of your calculation so speeding it up makes little difference for the total time. I have been running long DNA simulations this summer, and for a system with 15000 atoms I get 2ns/day with PME (2.5 ns/day w/o PME) on 8 cores of a dual quad core 2.5GHz intel E5420. No thermostat.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
If you intend to use or compare to the Drude polarizable FF in CHARMM, you must use the VV2 integrator. Otherwise, LEAPfrog is generally preferred, for NVE, NVT, NPT and for anisotropic ensembles such as NPAT or NPgammaT in most cases. I haven't used VV2 enough to understand it's weaknesses.
(PME is also preferred; force errors are typically 2 orders of magnitude less than spherical truncation methods. J.Phys.Chem.100:17011-17020 (1996).)
For what it's worth, I agree with Rick. The main advantage of the VV2 integrator is that it supports the Drude polarizable FF. If you're not using the Drude polarizable FF, it is wise to stick with the default LEAPfrog integrator, which is computationally efficient and thoroughly validated.
Same thing for the spherical truncation and the Berendsen thermostat. Generally spoken, everyone running MD simulations should bear in mind that every 2x increase in simulation time time allows you to you to access events only ~0.4kcal/mol higher in free energy, and only decreases errors in statistical properties by 1.4-fold. This implies that potentially compromising the accuracy of your results for a ~25% increase in speed is a clear example of a false economy! You can do these things if you're absolutely sure that you know what you're doing AND can explain it to your reviewers. Otherwise, stay away from Berendsen, PBC simulations without PME, and I would say even NVE (although that's my personal opinion).
Forget that I mentioned the Berendsen thermostat; the CPT thermostat is probably better and approximately equally efficient. I was just wondering why the VV2 integrator is so slow when I can't think of anything expensive that it's doing as part of the integration. Also I never mentioned using anything other than PME for the electrostatics on this system; other people on this forum brought that up.
I am using the CPT with the leapfrog integrator; it seems to be working well.
By the way, is there any way of recentering several segments together (other than joining them into a single segment)?
Not for image centering per se, but MERGE COOR (dynamc.doc) with RECENTER and ORIENT keywords and the appropriate second atom selection can be used to process trajectory files and make new ones.
A word from the original developer of VV2, here :-)
The VV2 code is slower because it was written using an "operator splitting" procedure, which is more concerned about getting the order of each operation right than having as few of them as possible. (For complicated integration algorithms, this is the only sane way to go...)
For instance, whereas the LEAP code may update velocities using: v = v + dt*F/m the VV2 code will often do something like: v = v + 0.5*dt*F/m do something else v = v + 0.5*dt*F/m
The code doesn't actually test to see if you are actually doing this "something else". That mean it will be slower if you are *not* -- which is often the case if you don't use that many of the features of VV2, such as Drude oscillators, multiple thermostats, barostat, SHAKE, etc.
An easy (but potentially messy) way to bring VV2 up to speed for "plain vanilla" dynamics would be to always check if the "something else" is done or not:
IF (doing something else) THEN v = v + 0.5*dt*F/m do something else v = v + 0.5*dt*F/m ELSE v = v + dt*F/m ENDIF
