[gmx-users] CG Lincs errors

Peter Kroon p.c.kroon at rug.nl
Fri Dec 16 10:02:55 CET 2016


As a note to Alex (and the rest of the list), the coarse-grained Martini
forcefield is usually run with timesteps between 20-40 fs. 15fs is
already rather low. I do agree that longer equilibration at low timestep
(5 or 10) might help.

Alternatively, Do you think a semiisotropic pressure coupling might be
applicable in this case, since it's an infinite collagen polymer? 


Peter (PhD in the Martini group)


On 16-12-16 00:21, Nash, Anthony wrote:
> Alex and Mark, thanks for the information. I’ll drop dt down,
> significantly, drop the temperature, and run it for a long time.
>
> Thanks for the ideas.
> Anthony 
>
>
> On 15/12/2016 21:54, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
> behalf of Alex" <gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
> behalf of nedomacho at gmail.com> wrote:
>
>> Mark is right, no two ways about it. For initial equilibration and
>> assessing preexisting structural strains try vacuum, _much_ smaller
>> timesteps and possibly low temperatures in vacuum, only then transfer to
>> solvent, etc. Algorithmically, LINCS requires convergence and you already
>> are using a pretty high LINCS order... From what I see, dt = 15 fs at 310K
>> looks like a cowboy mode simulation in this case.
>>
>> Alex
>>
>> On Thu, Dec 15, 2016 at 2:32 PM, Mark Abraham <mark.j.abraham at gmail.com>
>> wrote:
>>
>>> Hi,
>>>
>>> If a simulation isn't stable with a small time step (as I think you are
>>> saying) then moving to a larger time step is guaranteed to make that
>>> worse.
>>> Try an even smaller time step, for a long time, and see what happens. Or
>>> take a subset of your protein and see what happens. Or simulate in vacuo
>>> for a while. Your topology could be unsuited to your starting structure,
>>> e.g. some part is under a lot of tension that gets released at some
>>> point
>>> and no finite time step can in practice deal with the velocity of the
>>> recoil...
>>>
>>> Mark
>>>
>>> On Thu, 15 Dec 2016 23:06 Nash, Anthony <a.nash at ucl.ac.uk> wrote:
>>>
>>>> Hi all,
>>>>
>>>> I¹m hoping for some help. I¹m very sorry, this is a bit of a long one.
>>>>
>>>> I¹ve been struggling for almost a month trying to run a CG
>>> representation
>>>> of our all-atom model of a collagen protein (3 polypeptide chains in a
>>>> protein). Our original AMBER all-atom model had been successful
>>> modelling
>>>> using MD. I went on to use the latest version of Martinize.py with the
>>>> latest version of the MARTINI forcefield fields.
>>>>
>>>> After a little tweaking (the way AMBER names histidine residues), I
>>>> successful converted the molecule (approx 3100 amino acids) into a CG
>>>> representation. I successfully energy min the protein in vacuum to a
>>>> threshold of 500, and in solvent to a threshold of 750 using steepest
>>>> descent. Looking for a system at an energy min of a threshold around
>>> 300
>>> I
>>>> begin to see LINCS warnings. Observing the initial structure, there is
>>>> nothing obviously wrong with the bond network (both protein and
>>> polarized
>>>> CG water).
>>>>
>>>> I take the system that energy mins at 750 (protein-water mix, with no
>>>> fault reported), and went straight to NPT, 20fs step. Blew up. After a
>>> bit
>>>> of chatting with the MARTINI community, I¹ve started with an NVT
>>> ensemble,
>>>> beginning at 5s then through 10fs, 15fs, and 20fs. I only run for 1000
>>>> steps before switching. Keeping any of the simulations running for
>>> longer
>>>> throws lincs warnings followed by a segmentation fault from the
>>> warning:
>>>> "3 particles communicated to PME rank 7 are more than 2/3 times the
>>>> cut-off out of the domain decomposition cell of their charge group in
>>>> dimension x."
>>>>
>>>> Observing the trajectories of any of the extended simulations shows
>>> the
>>>> protein snapping like a rope, and always at the same place. I have
>>> watched
>>>> every trajectory at this point, using numerous energy min start
>>> points,
>>> to
>>>> try and understand why it is blowing up. I can¹t see any obvious
>>> reason.
>>> I
>>>> was told to consider how the temperature is changing. Below is an
>>> example
>>>> of the temperature and pressure from an NPT of 20fs step continued
>>> from
>>>> the very short 20fs step NVT simulation (hoping that perhaps CG
>>> without
>>>> pressure just doesn¹t behave happily; I was wrong).
>>>>
>>>>
>>>> TEMP:
>>>> Š
>>>> 6.630000  311.000336
>>>>     6.645000  311.371643
>>>>     6.660000  311.724213
>>>>     6.675000  313.878693
>>>>     6.690000  681558.937500
>>>>
>>>>
>>>> PRESSURE:
>>>> Š
>>>> 6.630000    3.559879
>>>>     6.645000    3.901433
>>>>     6.660000    3.589078
>>>>     6.675000    4.158611
>>>>     6.690000  81762.437500
>>>>
>>>> The final LINCS warning from this same run:
>>>>
>>>> Step 300, time 4.5 (ps)  LINCS WARNING
>>>> relative constraint deviation after LINCS:
>>>> rms 0.000035, max 0.003386 (between atoms 2125 and 2126)
>>>> bonds that rotated more than 45 degrees:
>>>>  atom 1 atom 2  angle  previous, current, constraint length
>>>>    2125   2126   68.3    0.2781   0.2691      0.2700
>>>>    2125   2127   45.9    0.2789   0.2701      0.2700
>>>>
>>>>
>>>> At this stage the structure ruptures as described above.
>>>>
>>>>
>>>> My NVT settings (with NPT included to save space) are:
>>>>
>>>> -----------------
>>>> title                    = Martini
>>>>
>>>> integrator               = md
>>>> dt                       = 0.015
>>>> nsteps                   = 1000
>>>> nstcomm                  = 100
>>>> ;comm-grps                       =
>>>>
>>>> nstxout                  = 1000
>>>> nstvout                  = 1000
>>>> nstfout                  = 0
>>>> nstlog                   = 1
>>>> nstenergy                = 1
>>>> nstxout-compressed       = 0
>>>> compressed-x-precision   = 0
>>>> ;compressed-x-grps        =
>>>> energygrps               = collagen solvent
>>>>
>>>> cutoff-scheme            = Verlet
>>>> nstlist                  = 20
>>>> ns_type                  = grid
>>>>
>>>> pbc                      = xyz
>>>> verlet-buffer-tolerance  = 0.005
>>>>
>>>> coulombtype              = PME ;reaction-field
>>>> rcoulomb                 = 1.1
>>>> fourierspacing       = 0.16 ;0.2  ;0.12
>>>>
>>>> epsilon_r                = 2.5 ;15      ; 2.5 (with polarizable water)
>>>> epsilon_rf               = 0
>>>> vdw_type                 = cutoff
>>>> vdw-modifier             = Potential-shift-verlet
>>>> rvdw                     = 1.1
>>>>
>>>> tcoupl                   = v-rescale ;berendsen ;v-rescale
>>>> tc-grps                  = collagen solvent
>>>> tau_t                    = 0.5 0.5 ;1.0 1.0
>>>> ref_t                    = 310 310
>>>>
>>>> Pcoupl                   = berendsen   ;parrinello-rahman
>>>> Pcoupltype               = isotropic
>>>> tau_p                    = 12.0  ; parrinello-rahman is more stable
>>> with
>>>> larger tau-p, DdJ, 20130422
>>>> compressibility          = 10e-4
>>>> ref_p                    = 1.0
>>>> refcoord_scaling         = com
>>>>
>>>> gen_vel                  = no
>>>> gen_temp                 = 310
>>>> gen_seed                 = 473529
>>>>
>>>> continuation = yes
>>>> constraints              = none
>>>> constraint_algorithm     = lincs
>>>> lincs-warnangle = 45
>>>> lincs-order=8
>>>> lincs-iter=4
>>>>
>>>>
>>>> ‹‹‹‹‹‹‹‹‹‹
>>>>
>>>> Every setting bar the lincs iter, order, warnangle were supplied with
>>> the
>>>> latest version of MARTINI. During many NVT runs I have adjusted the
>>> tau-t
>>>> to try and keep the thermostat from oscillating its way into infinity.
>>>>
>>>> I¹m curious, will an out of control thermostat break a structure, or
>>> will
>>>> a structure breaking (for what ever reason this structure is breaking)
>>>> cause the thermostat to go out of control?
>>>>
>>>> My only thought thought is the initial .itp file that Martinize
>>> created.
>>> I
>>>> informed the script that this was collagen, therefor it sets ³F² and
>>> the
>>>> corresponding bonded parameters. A human collagen is not perfect in
>>> its
>>>> helical structure. Could there be underlying forces contributed from
>>> badly
>>>> bonded backbone-backbone arrangements?
>>>>
>>>> [ atoms ]
>>>>     1    N0     1   GLY    BB     1  0.0000 ; F
>>>>     2    N0     2   LEU    BB     2  0.0000 ; F
>>>>     3    C1     2   LEU   SC1     3  0.0000 ; F
>>>>     4    N0     3   SER    BB     4  0.0000 ; F
>>>>
>>>>
>>>> Many thanks
>>>> Anthony
>>>>
>>>> Thanks
>>>> Anthony
>>>>
>>>> --
>>>> 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.
>>>>
>>> --
>>> 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.
>> -- 
>> 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