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

[Alex]

Also, I just looked at your topology. Where are your angles? Try to
implement them and then turn off constraints.

Alex


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

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

A> 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







-- 
Best regards,
 Alex                            mailto:nedomacho at gmail.com