[gmx-users] Negative average pressure value after NPT equilibration
Dario Akaberi
dario.aka17 at gmail.com
Thu Oct 27 14:30:28 CEST 2016
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).
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.
>
More information about the gromacs.org_gmx-users
mailing list