[gmx-users] Negative average pressure value after NPT equilibration

Justin Lemkul jalemkul at vt.edu
Thu Oct 27 20:47:54 CEST 2016



On 10/27/16 8:30 AM, Dario Akaberi wrote:
> Thank you very much for the answer, I thought that if the average pressure
> doesn't converge exactly of very close to 1 bar then something isn't
> working properly.
>
> I have one more question. In the mentioned Lysozyme tutorial the
> trajectories are corrected using the commands -pbc mol and -ur compact
> before calculating RMSD and radius of gyration. Is it normal to use these
> options? In my case when I visualize the simulation run with VMD there is a
> loop of my protein extends elongating in one direction for a couple of
> frames, then go back to the normal position, then extend again in other
> frames. If I use the option -pbc mol and ur compact the problem disappears
> (in particular the command I use is gmx_d file.tpr -f md_0_1.xtc -dt 20 -b
> 100 -o filextc -pbc mol -ur compact).
>

Yes, correcting for periodicity effects is pretty standard.  If you don't, and 
your simulation runs long enough, your protein will diffuse "out" one side and 
"back in" the other, leading to a huge spike in RMSD when the protein didn't 
actually do anything.

-Justin

> Thanks again for the help,
>
> Dario
>
> On Mon, Oct 24, 2016 at 4:25 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 10/23/16 10:07 AM, Dario Akaberi wrote:
>>
>>> So In my system I had 1 protein and a second smaller protein  associated
>>> to
>>> it. Since the second protein was not even complete, I eliminated it from
>>> the original PDB file I used to create the topology with pdb2gmx. I then
>>> repeated the whole process and now my average pressure after 500 ps of NPT
>>> equilibration is -3.31 with an Err. Est of 2.5, RMSD of 166.905, total
>>> drift -6.63.
>>>
>>> Is there anything I can do to get to the proper reference value of 1 bar?
>>> I
>>> tried a 1ns NPT equilibration which resulted in an average pressure of
>>> -4.77 and Err. EST of 3.8. So a longer equilibration is not gonna work In
>>> my opinion.
>>>
>>>
>> The RMSD (fluctuation) of the value means you likely can't establish with
>> certainty that the average pressure is statistically different from the
>> target of 1 bar.  Pressure fluctuates wildly, as the resources Mark noted
>> below tell you.  Rarely, if ever, will you see a value of exactly 1 on the
>> nanosecond scale.
>>
>> -Justin
>>
>>
>> On Sat, Oct 22, 2016 at 2:05 PM, Mark Abraham <mark.j.abraham at gmail.com>
>>> wrote:
>>>
>>> Hi,
>>>>
>>>> Negative pressure means the system wants to contract. But your
>>>> simulations
>>>> are probably too short to get a reliable measurement of pressure, nor to
>>>> equilibrate it if you were not already very close to the right density.
>>>> See
>>>> background at e.g. http://www.gromacs.org/Documentation/Terminology/
>>>> Pressure
>>>>  and
>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/
>>>> gmx-tutorials/lysozyme/07_equil2.html
>>>>
>>>>
>>>> Mark
>>>>
>>>> On Sat, Oct 22, 2016 at 1:34 PM Dario Akaberi <dario.aka17 at gmail.com>
>>>> wrote:
>>>>
>>>> Dear users
>>>>>
>>>>> I am trying to perform the MD simulation of a protein in water using the
>>>>> amberFF99sb-ildn force field.
>>>>>
>>>>> To so so I created the protein topology specifying TIP3P water. then I
>>>>> created a dodecahedron box with -d 1.0, I filled it using spc216.gro.
>>>>> and
>>>>> added 7 NA ions. This resulted in a system of  8837 water molecules plus
>>>>> the 7 ions and the protein. Minimization were carried out using
>>>>> conjugate
>>>>> gradient method (gromacs is compiled with double precision) which
>>>>>
>>>> converged
>>>>
>>>>> to emtol 0.01 in 48000 steps.The following nvt ensemble worked well (I
>>>>> think) and the resulting average temperature was 299.702 (ref
>>>>> temperature
>>>>> was 300 K) with Err. Est. = 0.32  RMSD 4.39922 and tot-Drift 1.58106.
>>>>> Now the problems begin. After NVT equilibration I performed an NPT
>>>>> equilibration step of 100 ps which resulted in an average pressure of
>>>>>
>>>> -50.
>>>>
>>>>> I then tried to equilibrate for a longer time: 150 and 200 ps with the
>>>>>
>>>> same
>>>>
>>>>> result. I also tried a 500 ps NPT equilibration which resulted in an
>>>>> average  pressure of -60 and average density 998.935 with Err. Est =
>>>>>
>>>> 0.12,
>>>>
>>>>> RMSD = 2.38 and Tot-Drift = 0.83 /kg/m^3).
>>>>>
>>>>> This is the NPT.mdp file I used:
>>>>>
>>>>> define                    = -DPOSRES
>>>>> integrator               = md        ; leap-frog integrator
>>>>> nsteps                   = 50000             ; 2 * 50000 = 100 ps
>>>>> dt                          = 0.002        ; 2 fs
>>>>>
>>>>> nstxout                   = 500        ; save coordinates every 1.0 ps
>>>>> nstvout                   = 500        ; save velocities every 1.0 ps
>>>>> nstenergy               = 500        ; save energies every 1.0 ps
>>>>> nstlog                     = 500        ; update log file every 1.0 ps
>>>>>
>>>>> continuation                = yes        ; Restarting after NVT
>>>>> constraint_algorithm        = lincs            ; holonomic constraints
>>>>> constraints                = all-bonds            ; all bonds (even
>>>>> heavy
>>>>> atom-H bonds) constrained
>>>>> lincs_iter                = 1                ; accuracy of LINCS
>>>>> lincs_order                = 4                ; also related to accuracy
>>>>>
>>>>> cutoff-scheme              = Verlet
>>>>> ns_type                   = grid
>>>>> nstlist                   = 10
>>>>> rlist                       = 1.0
>>>>> rcoulomb                = 1.0
>>>>> rvdw                    = 1.0
>>>>>
>>>>> coulombtype                = PME
>>>>> pme_order                = 4
>>>>> fourierspacing                = 0.16
>>>>>
>>>>> ; Temperature coupling is on
>>>>> tcoupl                = V-rescale                    ; modified
>>>>> Berendsen
>>>>> thermostat
>>>>> tc-grps                    = Protein Non-Protein
>>>>> tau_t                    = 0.1      0.1
>>>>> ref_t                    = 300       300            ; reference
>>>>> temperature, one for each group, in K
>>>>>
>>>>> pcoupl                    = Parrinello-Rahman            ; Pressure
>>>>> coupling on in NPT
>>>>> pcoupltype                = isotropic                    ; uniform
>>>>>
>>>> scaling
>>>>
>>>>> of box vectors
>>>>> tau_p                    = 2.0                ; time constant, in ps
>>>>> ref_p                    = 1.0                ; reference pressure, in
>>>>>
>>>> bar
>>>>
>>>>> compressibility             = 4.5e-5                    ; isothermal
>>>>> compressibility of water, bar^-1
>>>>> refcoord_scaling            = com                          * I get an
>>>>>
>>>> error
>>>>
>>>>> without this
>>>>>
>>>>> ; Periodic boundary conditions
>>>>> pbc                    = xyz                ; 3-D PBC
>>>>> ; Dispersion correction
>>>>> DispCorr                = EnerPres                    ; account for
>>>>>
>>>> cut-off
>>>>
>>>>> vdW scheme
>>>>> ; Velocity generation
>>>>> gen_vel                    = no                ; Velocity generation is
>>>>>
>>>> off
>>>>
>>>>>
>>>>> I don't know what I'm doing wrong, and unfortunately, I am new to
>>>>>
>>>> molecular
>>>>
>>>>> dynamics simulation which doesn't help. I tried to find a solution and
>>>>> found some other cases where people had a negative pressure after NPT
>>>>> equilibration step but it was like an average pressure of -2 not -60!
>>>>> Any help would be very appreciated.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Dario
>>>>> --
>>>>> 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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> 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

==================================================


More information about the gromacs.org_gmx-users mailing list