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

[hang yin]

Hi all,

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

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

Seeking for help. Thanks a lot.

Kevin
******************************************************************************
My itp file:

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

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

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

#ifdef FLEXIBLE
[ bonds ]
#endif

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

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

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

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

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

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

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

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