[gmx-users] Constraints applied for keeping rigid in ND
hang yin
yinhang503 at gmail.com
Wed Apr 22 09:21:37 CEST 2015
Hi Alex,
Sorry for late reply.
Because I could not find any reference of crystal CG mapping strategy, I
parameterized my ND based on 4-1 mapping of the Martini CG model.
Previously, I used default parameters and got an "oscillatory crystal" so I
thought maybe I should try constraints. After talking with you, I tried the
CG model with different values of angular force constant and managed to run
it at 20 fs time step. As long as its rmsd plot converged with the plot of
the atomic model, I could say that I got a rigid enough CG crystal.
Finally, I got angular constant 40 and bond constant 4500.
And I found that I am wrong previously. The results of rmsd of the same
structure in different time step basically converged with each other.
Because the system with large time step moved quicker in the VMD, I thought
incorrectly that the time step somehow related to the oscillatory behavior.
After I put out every frame of the trajectories and watched them in the
different speed corresponding to their time step, they were same actually.
Thanks for your advises. You are really kind and warmhearted. Sorry for my
bad English.
Hang
2015-04-21 1:42 GMT+08:00 Alex <nedomacho at gmail.com>:
> Hang,
>
>
> This is exactly the problem. Your system isn't set up properly, and you
> are trying to mask this fact behind constraints. By the way, oscillatory
> behavior does not necessarily mean the system is too soft. It may actually
> mean that your set of connstants is causing resonant or near-resonant
> oscillations, given the time steps you've got. For instance, it is good
> practice to check if the coarse-grain time step corresponds to the ratio
> between the bead mass and the average per-particle mass of all atoms
> assigned to a bead. In the case of homogeneous crystals, it is merely the
> number of atoms in one bead. If the number of atoms in one bead is five and
> you're hiking up the time step from 1 ps to 20 ps, there's likely going to
> be a problem.
>
>
> It is good that LINCS is screaming, it would be far more dangerous if it
> didn't.
>
>
> Diamond (as a crystal type) is one of the most stable ones, and wild
> oscillations are usually indicative of the setup being wrong. There are of
> course ways to try making LINCS/SHAKE work in this situation, but that's
> just bad science, in my opinion.
>
>
> Good luck!
>
>
> Alex
>
>
>
> >
>
> Hi Alex,
>
> When I ran my diamond with bonds and angles interactions, the beads in the
> diamond oscillated significantly, larger time step larger amplitude. The
> crystal shape could not conserved.
>
> Tomorrow I'll try SHAKE.
>
> Really thanks for your help! It is deep night here, so I may reply you on
> later morning. :)
>
> Hang
>
> 2015年4月21日 上午12:49于 "Alex" <nedomacho at gmail.com>写道：
>
>
> 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
> <http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List> before
> posting!
>
> >>
>
> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> <http://www.gromacs.org/Support/Mailing_Lists>
>
> >>
>
> >> * For (un)subscribe requests visit
>
> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
> <https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users> or
> send a mail to gmx-users-request at gromacs.org
> <gmx-users-request at gromacs.org>.
>
> >
>
> >
>
>
>
>
>
> --
>
> Best regards,
>
> Alex mailto:nedomacho at gmail.com
> <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