On 12/21/16 8:57 AM, Michael K. Gilson wrote:
Dear All, Below are more details on our efforts to maximize energy conservation in NVE runs with gromacs, while maintaining performance. I'd appreciate advice on the following: 1) Is this level of drift about "as good as it gets"? 2) Or, are there ways to make it even lower? 3) Are there settings we can dial back to improve speed without substantially increasing drift? Many thanks, Mike Gilson System: The system comprises 68679 waters, one protein with 6810 atoms, 8 sodium ions; total 75497 atoms The input file is below. We have deliberately turned off COM restraints because we are interested in how the protein moves.
Removal of net center-of-mass motion is not a COM restraint; your protein will diffuse fine with comm-mode = Linear. The removal of net COM motion is to avoid the flying ice cube effect. I suspect that failing to remove net COM motion is actually detrimental to conservation of energy.
-Justin
As noted before, for lincs-iter=4, the energy slope is approximately -0.01976 kJ/mol/ps; for lincs-iter=5, the slope is approximately 0.00991 kJ/mol/ps. Graphs of the energy for these respective conditions are provided here: https://drive.google.com/open?id=0B0zEqb9pykUWRHlqSEJTa2lwYkk https://drive.google.com/open?id=0B0zEqb9pykUWMWl0OEFtNGdVbzg MDP: ;DON'T POSITION RESTRAIN THE PROTEIN: NO -DPOSRES ;run parameters integrator = md ; leap-frog integrator nsteps = 10000000 ; 0.001ps * 10000000 = 10 ns dt = 0.001 ; 0.001ps smaller time step for NVE comm-mode = None ; let center of mass move, calculate diffusion constant ; Output control nstxout = 20000 ; save coordinates every 0.5 ns nstvout = 20000 ; save velocities every 0.5 ns nstenergy = 20000 ; save energies every 0.5 ns nstlog = 20000 ; update log file every 0.5 ns energygrps = protein non-protein ; Bond parameters continuation = yes ; starting after nvt.mdp constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds are constrained lincs_iter = 4 ; NVE requires high lincs-iter lincs_order = 4 ; accuracy ; Neighborsearching cutoff-scheme = Verlet ns_type = grid ; search neighboring grid cells nstlist = 10 ; Frequency to update the neighbor list and long range forces rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm) rvdw = 1.2 ; short-range van der Waals cutoff (in nm) verlet-buffer-tolerance = 0.0000005 ; kJ/mol/ps per particle pbc = xyz ; periodic boundary conditions DispCorr = EnerPres ; account for cut-off vdW scheme ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT ; Pressure coupling is off pcoupl = no ; fixed box size gives constant V ; Temperature coupling is off + no simluated annealing tcoupl = no ; NVE conditions ; Velocity generation gen_vel = no ; don't assign velocities: keep from thermalization
-- ================================================== Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul ================================================== -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.