[gmx-users] Unexpected lincs errors for OH groups in protein

[Erik Marklund]

Hi all,

This appears to happen only with group cutoffs, and not with verlet. The (probable) reason for that is that SOL-SOL short-range coulomb energy is zero when I use group cutoffs, which can't be right. I filed a new redmine issue: http://redmine.gromacs.org/issues/1426.

Kind regards,
Erik Marklund

On 22 Jan 2014, at 15:00, Erik Marklund <erik.marklund at chem.ox.ac.uk> wrote:

> Dear GMX-users,
> 
> I am initiating a series of simulations of a dimeric protein in water, i.e. something that should be fairly trivial. I use gromacs 4.6.5, vsites for the hydrogens, amber99sb-ildn for the protein, and tip4p-Ew for the water. With a time step of 0.004 I keep getting occasional lincs errors for OH groups that swing around to violently, to my surprise. There are no clashes in the structure that I can see and the OH groups are mostly at the protein-water interface and as such quite unconstrained by the surroundings. I am using a similar setup to what I have successfully been using with gromacs 4.5.X on numerous occasions, but non-bonded mdp options have changed so I suspect that's where the problem is rooted, but so far I have been unable to overcome these errors. EM proceeded just fine, a 100 ps NVT run with position restraints and a 0.002 fs time step, and a 1 ns NVT run with time step 0.002 fs preceded this NVT run  Is there anything obviously wrong with the following combination of mdp parameters?
> 
> integrator               = md
> ; Start time and timestep in ps
> tinit                    = 0
> dt                       = 0.004
> nsteps                   = 250000
> …
> 
> cutoff-scheme            = group
> ; nblist update frequency
> nstlist                  = 10
> …
> 
> rlist                    = 1.0
> ; long-range cut-off for switched potentials
> rlistlong                = -1
> nstcalclr                = -1
> 
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype              = PME
> coulomb-modifier         = Potential-shift-Verlet ; <- moot
> rcoulomb-switch          = 0
> rcoulomb                 = 1.0
> ; Relative dielectric constant for the medium and the reaction field
> epsilon-r                = 1
> epsilon-rf               = 0
> ; Method for doing Van der Waals
> vdw-type                 = Cut-off
> vdw-modifier             = Potential-shift
> ; cut-off lengths       
> rvdw-switch              = 0
> rvdw                     = 0.8
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                 = EnerPres
> …
> 
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling  
> tcoupl                   = v-rescale
> nsttcouple               = -1
> nh-chain-length          = 10
> print-nose-hoover-chain-variables = no
> ; Groups to couple separately
> tc-grps                  = Protein      Water_and_ions
> ; Time constant (ps) and reference temperature (K)
> tau-t                    = 0.2  0.2
> ref-t                    = 300  300
> ; pressure coupling     
> Pcoupl                   = Berendsen
> pcoupltype               = Isotropic
> nstpcouple               = -1
> …
> 
> ; GENERATE VELOCITIES FOR STARTUP RUN
> gen-vel                  = no
> gen-temp                 = 300
> gen-seed                 = 173529
> 
> ; OPTIONS FOR BONDS    
> constraints              = all-bonds
> ; Type of constraint algorithm
> constraint-algorithm     = Lincs
> ; Do not constrain the start configuration
> continuation             = no
> ; Use successive overrelaxation to reduce the number of shake iterations
> Shake-SOR                = no
> ; Relative tolerance of shake
> shake-tol                = 0.0001
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 6
> ; Number of iterations in the final step of LINCS. 1 is fine for
> ; normal simulations, but use 2 to conserve energy in NVE runs.
> ; For energy minimization with constraints it should be 4 to 8.
> lincs-iter               = 1
> ; Lincs will write a warning to the stderr if in one step a bond
> ; rotates over more degrees than
> lincs-warnangle          = 30 
> ...
> 
> 
> I tried various alterations to this setup, e.g. rvdw = 1.0, but still no success. Any help coming to grips with the new non-bonded stuff is much welcome.
> 
> Kind regards,
> 
> Erik Marklund, PhD
> Postdoctoral Research Associate
> 
> Department of Chemistry
> Physical & Theoretical Chemistry Laboratory
> University of Oxford
> South Parks Road
> Oxford
> OX1 3QZ
> 
> -- 
> 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-request at gromacs.org.