[gmx-users] energy minimisation - LINCS WARNING
Nash, Anthony
a.nash at ucl.ac.uk
Mon Oct 3 00:11:35 CEST 2016
Hi all,
I had a homology/de-novo model .pdb converted into .gro, solvated,
neutralised and now I¹m going through a series of energy minimisation
steps. Unfortunately, during energy minimisation I get LINCS WARNINGS
(angle relative constraint deviation). The the naked eye, the atoms
involved don¹t look untoward.
When I used LINCS all-bond I get the following error:
-----------------------------------
Step -1, time -0.001 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.001657, max 0.033500 (between atoms 1362 and 1363)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
When I reduced this to only be the hydrogen bonds (h-bonds) I get:
-----------------------------------
Step 20, time 0.02 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000000, max 0.000000 (between atoms 1322 and 1324)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1360 1361 37.5 0.1010 0.1010 0.1010
[continues with more errors until Segmentation fault 11]
I have tried both constraints with a reduced and increased deltaStep
(0.2,0.02,0.002,0.0002) no change. I have also taken constraints off, this
does work but then fails again when constraints are returned. I have tried
both Steep and cg, both fail with the same result.
The .mdp file I used (taken from the latest release of Justin¹s tutorials)
is:
; energy.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent
minimization)
emtol = 700.0 ; Stop minimization when the maximum force
< 1000.0 kJ/mol/nm
emstep = 0.002 ; Energy step size
nsteps = 200000 ; Maximum number of (minimization)
steps to perform
; Parameters describing how to find the neighbors of each atom and how to
calculate the interactions
cutoff_scheme = verlet
nstlist = 20 ; Frequency to update the neighbor list
and long range forces
ns_type = grid ; Method to determine neighbor list
(simple, grid)
rlist = 1.4 ; Cut-off for making neighbor list (short
range forces)
coulombtype = PME ; Treatment of long range electrostatic
interactions
rcoulomb = 1.4 ; Short-range electrostatic cut-off
rvdw = 1.4 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions (yes/no)
;ADDED THIS BIT MYSELF
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds; all-bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 8 ; also related to accuracy
The .pdb was generated using I-TASSER. I have tried refining the hydrogen
placement via -ignh during pdb2gmx, and I also tried ³Protein Preparation²
(bond ordering - again, hydrogen bond placement) in Schrodinger. All the
same.
I was hoping that I could actually watch a rough-and-ready minimisation
using Avogadro, however, the release I have keeps crashing whilst trying
to open large structures.
Any suggestions on what I could try next?
Thanks
Anthony
Dr Anthony Nash
Department of Chemistry
University College London
Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp
On 02/10/2016 19:06, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
behalf of Karel de Vries"
<gromacs.org_gmx-users-bounces at maillist.sys.kth.se on behalf of
kjcdevries at gmail.com> wrote:
>Mark,
>
>I understand your concern.
>I need to have the hydrogens for my final answer, so an all-atom force
>field is definitely inevitable. Justification for choosing OPLS over the
>other all-atom force fields is a little more tricky, but I have found
>literature where they used OPLS-aa for similar purposes so I figured it
>should work. As for the presence of a quaternary carbon, I expect that to
>be the "alkane C", opls_139. Given that opls_135 through 138 are "alkane
>CHn", this sounds like the most reasonable interpretation of opls_139.
>
>I will try to find better documentation for the exact meaning of these
>atom
>types though, thanks.
>
>On Fri, Sep 30, 2016 at 1:49 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 9/30/16 7:14 AM, Karel de Vries wrote:
>>
>>> Hi Justin,
>>>
>>> Thanks for your answer; I assumed that TopolGen is not intentionally
>>> incorrect. What I meant to ask was whether it was intentionally
>>>different
>>> from what the comments in the atomtypes.atp file indicate. I thought
>>>that
>>> perhaps, based on past experience, you know that some choices work
>>>better
>>> than others and included this into your script. I would not want to
>>>undo
>>> this by cooking up my own assignments based on no past experience at
>>>all.
>>>
>>>
>> Nope, it's just based on the normal principal of trying to construct
>> additive force fields from building blocks. Take what's known, assign a
>> best match, otherwise default to some stock parameter, then clearly warn
>> the user that the script can be dangerous :)
>>
>> -Justin
>>
>>
>> Karel
>>>
>>> On Fri, Sep 30, 2016 at 12:00 PM, Justin Lemkul <jalemkul at vt.edu>
>>>wrote:
>>>
>>>
>>>>
>>>> On 9/30/16 3:35 AM, Karel de Vries wrote:
>>>>
>>>> Hi again,
>>>>>
>>>>> I think the problem is in the force field description that I put in
>>>>>the
>>>>> aminoacids.rtp file. Like I said, I just got the output from Justin's
>>>>> TopolGen (V1.1) script. Upon comparison of the opls_xxx choices to
>>>>>the
>>>>> comments in the atomtypes.atp file, a few strange assignments turn
>>>>>up.
>>>>> For
>>>>> example, the Os are taken to be from the C=O in an amide. However,
>>>>>there
>>>>> is
>>>>> no nitrogen in PMMA at all. There is a C=O in PMMA, but it's an ester
>>>>> group. There are a few more doubtful choices, e.g. a C bound to four
>>>>> other
>>>>> Cs is assigned the "alkane CH3" atom.
>>>>>
>>>>> I have two questions.
>>>>> One concerns the TopolGen script: are these seemingly incorrect
>>>>> assignments
>>>>> intentional? Or is this just a limitation in the script? There is
>>>>> another
>>>>> TopolGen topology for PMMA online, at
>>>>> http://gromacs-oplsaa-topologies.wikia.com/wiki/PMMA_3-mer, which has
>>>>> slightly different but similarly doubtful assignments.
>>>>>
>>>>>
>>>>> I would never write a program that is intentionally incorrect.
>>>>>TopolGen
>>>> tries to guess the best parameter based on connectivity, and
>>>>sometimes it
>>>> guesses wrong, because it is pretty much impossible to design a script
>>>> that
>>>> will account for the vast chemical space that OPLS-AA encompasses. I
>>>>try
>>>> to make this abundantly clear in the final message printed to the
>>>>screen
>>>> by
>>>> TopolGen:
>>>>
>>>> "Output topology has been written. An attempt has been made to assign
>>>> charges and atom types based on existing functional groups, but they
>>>>may
>>>> not be correct. No charge calculations or other parameterization
>>>> calculations have been done. Guesses have been made for charge groups.
>>>> Please inspect and correct the topology before using it in any
>>>> simulations. The author of the script does NOT guarantee accuracy or
>>>> usability of any of the content; TopolGen was written as a convenience
>>>> for
>>>> outputting a skeleton topology, and nothing more."
>>>>
>>>> -Justin
>>>>
>>>>
>>>> The second is about the opls-aa force field. Is there any other
>>>>
>>>>> documentation on which assignments to make than just the comments in
>>>>>the
>>>>> atomtypes.atp file?
>>>>>
>>>>>
>>>>> Regarding Chaofu Wu's reply to my original email; thank you very much
>>>>> for
>>>>> your comments. I accidentally sent a template file, not one used for
>>>>>an
>>>>> actual run. I run gromacs with lower nst* values; when I need
>>>>> statistics I
>>>>> make sure there are about 1000 frames in the energy file. Values I
>>>>> adjust
>>>>> from the parameters I sent are nsteps, nst*, ref_t, ref_p, and for
>>>>>the
>>>>> initial run also gen_vel / gen_temp. I'm sorry for the confusion.
>>>>>
>>>>> Karel de Vries
>>>>>
>>>>>
>>>>> On Wed, Sep 28, 2016 at 11:43 AM, Karel de Vries
>>>>><kjcdevries at gmail.com>
>>>>> wrote:
>>>>>
>>>>> Hello all,
>>>>>
>>>>>>
>>>>>> I'm trying to simulate bulk PMMA, poly(methyl methacrylate) in
>>>>>>Gromacs.
>>>>>> I'm using the opls-aa force field. I added a PMMA unit cell to the
>>>>>> aminoacids.rtp file (listing below). I'm running a simulation with
>>>>>>30
>>>>>> molecules of 50 unit cells each. I start by putting straight strands
>>>>>> into
>>>>>> the box with gmx insert-molecules.
>>>>>>
>>>>>> I think that my input parameters are incorrect. My first reason for
>>>>>> thinking this is that it is very difficult to equilibrate this
>>>>>>system:
>>>>>> it
>>>>>> takes a lot of computation time to get the internal pressure to
>>>>>>equal
>>>>>> the
>>>>>> external one during NPT even at high temperatures. The second
>>>>>>reason is
>>>>>> that I find unrealistic values for the glass transition temperature
>>>>>>and
>>>>>> density of PMMA at room temperature.
>>>>>>
>>>>>> I have tried several combinations of NVT and NPT simulations at
>>>>>> different
>>>>>> temperatures to get the system to equilibrate. The most successful
>>>>>>one
>>>>>> is
>>>>>> as follows. I start with an NPT simulation at 1000 bar of isotropic
>>>>>> pressure at 900 K using Berendsen barostat with tau_p=5 and
>>>>>> compressibility
>>>>>> 4.5e-5. Then I slowly decreased the pressure to 1 bar in 8 steps of
>>>>>>10
>>>>>> ns
>>>>>> each. Finally I ran a 40 ns NPT equilibration step at 1 bar at the
>>>>>>end.
>>>>>> The
>>>>>> average pressure tensor over this 40 ns run is as follows:
>>>>>> Pressure (bar)
>>>>>> 9.67017e-01 -4.24140e-01 -1.51426e+00
>>>>>> -4.24137e-01 -2.97605e+00 -3.38247e+00
>>>>>> -1.51425e+00 -3.38246e+00 2.07151e+00
>>>>>> I'm not sure how good this is, but it's much better than anything
>>>>>>else
>>>>>> I've done.
>>>>>>
>>>>>> I then cool down from 900K to 550K in steps of 50K, 5 ns per
>>>>>> simulation.
>>>>>> I
>>>>>> do not expect the average cooling rate (0.01 K/ps) to be
>>>>>> unrealistically
>>>>>> fast: 550K is well above the glass transition temperature of real
>>>>>>PMMA
>>>>>> and
>>>>>> I sometimes see people in literature use similar rates during
>>>>>>annealing
>>>>>> in
>>>>>> other polymers. I then cool down from 550K to 300K in steps of 10K,
>>>>>>5
>>>>>> ns
>>>>>> per simulation. The average pressure tensor contains values that
>>>>>>keep
>>>>>> getting further away from zero. Moreover, I find a glass transition
>>>>>> temperature of around 470K, while for real PMMA it should be around
>>>>>> 400K.
>>>>>> The final density is 1.01 g/cm3, while it should be around 1.18
>>>>>>g/cm3.
>>>>>>
>>>>>> I'm starting to think there is something wrong with my parameters
>>>>>>here,
>>>>>> either in the mdp file or in the force field parameters. Has anybody
>>>>>> here
>>>>>> got any experience with PMMA and can they tell me what I'm doing
>>>>>>wrong?
>>>>>>
>>>>>> Thanks a lot in advance.
>>>>>> Karel de Vries
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> This is what I added to the .rtp file in the force field. I have
>>>>>> similar
>>>>>> definitions for the end chains, which have a H connected to C1 and
>>>>>>C2,
>>>>>> respectively. This was based on output by Justin Lemkul's TopolGen
>>>>>> script,
>>>>>> except I took the charges from an external source (don't remember
>>>>>> where I
>>>>>> found them, unfortunately).
>>>>>> [ atoms ]
>>>>>> C1 opls_135 -0.376 1
>>>>>> C2 opls_135 0.009 1
>>>>>> H3 opls_140 0.212 1
>>>>>> H4 opls_140 0.155 1
>>>>>> C5 opls_145 0.650 2
>>>>>> C6 opls_135 -0.288 3
>>>>>> H7 opls_140 0.096 3
>>>>>> H8 opls_140 0.096 3
>>>>>> H9 opls_140 0.096 3
>>>>>> O10 opls_236 -0.286 2
>>>>>> C11 opls_166 -0.161 2
>>>>>> O12 opls_236 -0.563 2
>>>>>> H13 opls_140 0.120 2
>>>>>> H14 opls_140 0.120 2
>>>>>> H15 opls_140 0.120 2
>>>>>>
>>>>>> [ bonds ]
>>>>>> C1 C2
>>>>>> C1 H3
>>>>>> C1 H4
>>>>>> C11 H13
>>>>>> C11 H14
>>>>>> C11 H15
>>>>>> C11 O10
>>>>>> C6 H7
>>>>>> C6 H8
>>>>>> C6 H9
>>>>>> C2 C5
>>>>>> C2 C6
>>>>>> C5 O10
>>>>>> C5 O12
>>>>>> C2 +C1
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> These are typical parameters I've got in the mdp file:
>>>>>> integrator = md
>>>>>> tinit = 0
>>>>>> dt = 0.001
>>>>>> nsteps = 1000000 ; 1 ns
>>>>>>
>>>>>> emtol = 100
>>>>>> emstep = 0.001
>>>>>> niter = 20
>>>>>> fcstep = 0
>>>>>> nstcgsteep = 1000
>>>>>> nbfgscorr = 10
>>>>>>
>>>>>> cutoff-scheme = Verlet
>>>>>> nstxout = 1000000
>>>>>> nstvout = 1000000
>>>>>> nstfout = 1000000
>>>>>> nstlog = 2000000
>>>>>> nstenergy = 10000
>>>>>> nstxout-compressed = 10000
>>>>>> compressed-x-precision = 1000
>>>>>> rcoulomb = 1.00
>>>>>> nstlist = 45
>>>>>> ns_type = grid
>>>>>> pbc = xyz
>>>>>> rlist = 1.00
>>>>>> vdw-type = cut-off
>>>>>> vdw-modifier = Potential-switch
>>>>>> rvdw-switch = 0.8
>>>>>> rvdw = 1.00
>>>>>>
>>>>>> Tcoupl = v-rescale
>>>>>> tc-grps = System
>>>>>> tau_t = 0.2
>>>>>> ref_t = 900
>>>>>> Pcoupl = berendsen
>>>>>> Pcoupltype = isotropic
>>>>>> tau_p = 5.0
>>>>>> compressibility = 4.5e-5
>>>>>> ref_p = 1.0
>>>>>>
>>>>>> gen_vel = no
>>>>>> gen_temp = 900
>>>>>> gen_seed = 1000
>>>>>>
>>>>>> DispCorr = EnerPres
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>> ==================================================
>>>>
>>>> Justin A. Lemkul, Ph.D.
>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>
>>>> Department of Pharmaceutical Sciences
>>>> School of Pharmacy
>>>> Health Sciences Facility II, Room 629
>>>> University of Maryland, Baltimore
>>>> 20 Penn St.
>>>> Baltimore, MD 21201
>>>>
>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>
>>>> ==================================================
>>>> --
>>>> 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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>
>> Department of Pharmaceutical Sciences
>> School of Pharmacy
>> Health Sciences Facility II, Room 629
>> University of Maryland, Baltimore
>> 20 Penn St.
>> Baltimore, MD 21201
>>
>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>> http://mackerell.umaryland.edu/~jalemkul
>>
>> ==================================================
>> --
>> 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