[gmx-users] Constraints applied for keeping rigid in ND

[Alex]

Is there a particular reason you're applying LINCS constraints to your
diamond structure? Usually, LINCS gives convergence errors/warnings if the system is poorly
equilibrated and/or your time step is too large. The fact that your
system was "soft" without constraints actually suggests poor setup,
which once again could include the time step.

Diamond structure, if properly set up, should be quite rigid, even if
your bond/angle terms are somewhat off, because, well, it's the
diamond structure.

Alex


hy> Hi all,

hy> My research focuses on the dynamics of a nanodiamond(ND) in a biological
hy> environment. I am working with gromacs version 5.0.2 and using Martini CG
hy> ff. In order to keep the rigidity of ND, I applied constraints as the stiff
hy> bond in the .itp file. But my system crashed in the NVT run with LINCS
hy> WARNING. So I tried to use a small ND, an octahedron, in the water box for
hy> debugging. It did work well. However, after I increased its size(just add 3
hy> beads actually), it crashed with same input parameters.

hy> How it could happen? I am the beginner using gromacs and CG model. I think
hy> the problem is due to the constraints because the same system with defined
hy> bond interaction worked well(but the ND in this situation is soft). I am
hy> just wondering that is there any limitation for applying constraints for
hy> crystal packing system or anything wrong with my configuration?

hy> Seeking for help. Thanks a lot.

hy> Kevin
hy> ******************************************************************************
hy> My itp file:

hy> ;itp file of ND
hy> [ moleculetype ]
hy> ; molname nrexcl
hy>   CG-ND         1

hy> [ atoms ]
hy> ; id type resnr residu atom cgnr charge
hy>   1 P1   1     FO   CD1   1     0.000000
hy>   2 C1   1     FH   CD2   2     0.000000
hy>   3 P1   1     FO   CD3   3     0.000000
hy>   4 C1   1     FH   CD4   4     0.000000
hy>   5 C1   2     FH   CD1   5     0.000000
hy>   6 C1   2     FH   CD2   6     0.000000
hy>   7 C1   2     FH   CD3   7     0.000000
hy>   8 P1   2     FO   CD4   8     0.000000
hy>   9 P1   3     FO   CD1   9     0.000000
hy>  10 P1   3     FO   CD2   10    0.000000

hy> [ constraints ]
hy> ; i  j funct length force.c.

hy> #ifdef FLEXIBLE
hy> [ bonds ]
hy> #endif

hy>  1    2     1 0.317 1000000
hy>  1    3     1 0.317 1000000
hy>  1    4     1 0.317 1000000
hy>  2    3     1 0.317 1000000
hy>  2    4     1 0.317 1000000
hy>  3    4     1 0.317 1000000
hy>  3    5     1 0.317 1000000
hy>  3    6     1 0.317 1000000
hy>  4    5     1 0.317 1000000
hy>  4    6     1 0.317 1000000
hy>  5    6     1 0.317 1000000
hy>  5    7     1 0.317 1000000
hy>  5    8     1 0.317 1000000
hy>  6    7     1 0.317 1000000
hy>  6    8     1 0.317 1000000
hy>  7    10   1 0.317 1000000
hy>  7    8     1 0.317 1000000
hy>  7    9     1 0.317 1000000
hy>  8    10   1 0.317 1000000
hy>  8    9     1 0.317 1000000
hy>  9    10   1 0.317 1000000
hy> ******************************************************************************
hy> My mdp file:

hy> title = debugging for ND
hy> ; Run parameters
hy> integrator = md ; leap-frog integrator
hy> nsteps = 2500        ; 20 * 2500 = 500 ps (0.05 ns)
hy> dt    = 0.02 ; 20 fs
hy> ; Output control
hy> nstxout = 1 ; save coordinates
hy> nstvout = 1 ; save velocities
hy> nstenergy = 10 ; save energies
hy> nstlog = 10 ; update log file
hy> nstcomm                  = 10
hy> fcstep          = 0
hy> nstcalcenergy   = 10

hy> ; Bond parameters
hy> unconstrained_start      = no
hy> constraints              = none
hy> constraint_algorithm     = Lincs
hy> lincs_order              = 4
hy> lincs_warnangle          = 30
hy> lincs-iter               = 1
hy> morse                    = no

hy> ; Neighborsearching
hy> cutoff-scheme            = verlet
hy> verlet-buffer-drift      = 0.001
hy> ns_type = grid ; search neighboring grid cels
hy> nstlist = 10    ;
hy> rlist = 1.4 ; short-range neighborlist cutoff (in nm)
hy> rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
hy> rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
hy> vdw_type        = cutoff
hy> rvdw_switch     = 0.9
hy> vdw-modifier             = Potential-shift

hy> ; Electrostatics
hy> coulombtype = Reaction_field
hy> coulomb_modifier= Potential-shift
hy> rcoulomb_switch = 0.0
hy> epsilon_r = 15
hy> epsilon_rf      = 0

hy> ; Temperature coupling is on
hy> tcoupl = Berendsen
hy> tc-grps = system   ; three coupling groups - more accurate
hy> tau_t = 1.0        ; time constant, in ps
hy> ref_t = 310         ; reference temperature, one for each group, in K
hy> ; Pressure coupling is on
hy> pcoupl = no    ; Pressure coupling on in NPT
hy> pcoupltype = isotropic    ; uniform scaling of x-y-z box vectors
hy> tau_p = 4.0        ; time constant
hy> ref_p = 1.0        ; reference pressure
hy> compressibility = 1e-5 ; isothermal compressibility, bar^-1

hy> ; Periodic boundary conditions
hy> pbc    = xyz ; 3-D PBC

hy> ; Velocity generation
hy> gen_vel = no
hy> gen_temp                 = 310
hy> gen_seed                 = -1