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

Alex nedomacho at gmail.com
Mon Apr 20 18:49:24 CEST 2015


Hi Hang,

Are the angles in the system conserved via artificial bonds? In other
words, when you run your diamond, does it lose crystal shape, or does it
stay diamond? I am just trying to understand what you call "soft". You can
always try SHAKE instead of LINCS, if you insist on constraints, i'm just
trying to see if these constraints are masking something really wrong with
the unconstrained setup. You see, noone uses constraints in solid-state
simulations, the structure must be conserved, unless something is seriously
wrong.

Alex
On Apr 20, 2015 10:34 AM, "hang yin" <yinhang503 at gmail.com> wrote:
>
> Hi Alex,
>
> Thanks for your reply.
>
> I used large time step because it is the default value(20 ~40 ps) for the
Martini Coarse Grained ff. Applied that large time step, the angles and
bonds have been proved that they could not keep ND rigid enough. And for a
very large system I want to build in future, the CG model have to be used,
meaning that I have to solve the rigidity problem in a large time step
setup. Hence I really want to use constraints here.
>
> Hang
>
>
> 2015-04-21 0:22 GMT+08:00 Alex <nedomacho at gmail.com>:
>>
>> 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
>>
>> --
>> 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.
>
>


More information about the gromacs.org_gmx-users mailing list