[gmx-users] Pressure Question

Thu Nov 6 06:16:46 CET 2014

Well, it is definitely *not* "totally irrelevant to talk about pressure" 
when doing simulations in the NVT (canonical) or NVE (microcanonical) 
ensembles. Pressure, like temperature, volume, energy, numbers of 
particles, etc, is a thermodynamic property which is *always* defined for 
*any* system in equilibrium. These parameters characterize the 
thermodynamic state, regardless of the ensemble that you choose to 
describe the system in terms of statistical mechanics. And the microscopic 
counterparts of these thermodynamic parameters are defined in such a way 
that their ensemble average must necessarily equal the thermodynamic 
value, regardless of ensemble (although their fluctuations are 
ensemble-dependent, becoming asymptotically identical only for macroscopic 
systems).

In particular, the virial-based "instantaneous pressure" (call it P') 
computed in simulations has its ensemble average equal to the 
thermodynamic pressure P (check any good book on molecular simulation). 
But, as others already pointed out, this P' is well-known to show 
extremelly large fluctuations, meaning that its average computed from the 
simulation has usually a very large statistical spread. In other words, 
although the ensemble average of P' is strictly equal to P, its simulation 
average is a random variable that often shows large deviations from P 
(especially for short simulations). To get an idea of what is an 
acceptable error for the average of P', you may look at its distribution 
histogram in the NPT simulation.

As for the equilibration of the system, the only thing that matters is 
which termodynamic state you are aiming at and what is the best way to get 
there. For example, if you choose NVT but happen to start with a volume 
which is a bit too large (eg, because the parameterized model acquires a 
higher density than the true experimental value that you assumed when 
setting the box size), you may get into trouble because the system may 
then want to separate into two phases but, being unable to do so in a 
small simulation box, ends up in a weird metastable state (eg, if you take 
an amount of water into a syringe, seal the tip, and then further pull the 
piston, you will get an "empty" region that is actually filled with water 
vapor, because having only liquid water filling that volume at that 
temperature is not thermodynamically stable). So, it is usually a good 
idea equilibrate in NPT, because the system finds its proper density at 
some temperature and pressure, whose regions of interest you usually know 
for the system you are studying. Once a "good" volume is found for some 
relevant P and T values, you can do the same for the energy: use NVT and 
let the system find a "good" energy for that T value, moving then to NVE. 
So, if for some reason you really want to do a simulation in the NVE 
ensemble, the suggested sequential procedure NPT > NVT > NVE sounds 
reasonable to me. Actually, I have a simpler suggestion: just run an NPT 
simulation, look at the 2D distribution histogram of V and E values, 
choose one representative snapshot that is in the central region of that 
distribution, and use that snapshot to run NVE. In any case, unless you 
have some experimental indication of the material density (N/V) *and* of 
the energy density (E/V) of your system, which would be extremelly 
unusual, you will have to follow some kind of approach similar to this. Of 
course, we may also ask why you think you need an NVE simulation, but that 
is an entirely different question...

Best,
Antonio

-- 
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica - EAN, 2780-157 Oeiras, Portugal
phone: +351-214469619         email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista
--------------------------------------------------------------------------

On Wed, 5 Nov 2014, Johnny Lu wrote:

> Well... I finally have to accept that I really need NVE. Some forum posts
> suggested to run NPT, then NVT, and finally NVE.
>
> For ideal gas law, if we fix the number of molecules and the temperature,
> then fixing the volume at different value means changing the pressure.
> So, I guess pressure is still meaningful, if I have to decide the correct
> volume to run the NVT simulation.
>
> For a large number of ideal gas molecules, averages of all thermodynamic
> variables should be defined, regardless of which three of them is fixed or
> which ensemble do I choose to represent the state.
> For a small number of molecules, the distribution of the thermodynamic
> variables should still have a definition (so, the average and the
> distribution of pressure is still meaningful).
>
> For a protein simulation, it still has a partition function, and equation
> of state.
>
> I mainly want to know how much error in pressure (deviation from 1 bar) is
> acceptable.
>
> On Wed, Nov 5, 2014 at 2:35 PM, Téletchéa Stéphane <
> stephane.teletchea at univ-nantes.fr> wrote:
>
>> Le 04/11/2014 18:00, Johnny Lu a écrit :
>>
>>  Hi.
>>>
>>> If my NVT simulation of a protein in 30k molecules of water has a pressure
>>> of 11 bar (error 0.5 bar from g_energy), will the dynamics (not
>>> distribution of conformations) change enough that the mechanism inferred
>>> from this simulation be significantly more unreliable than the mechanism
>>> inferred from a 1 bar simulation? (Will the reviewers cut my paper into
>>> ribbons?)
>>>
>>> Thanks again.
>>>
>>
>> Hi,
>>
>> Considering only your "NVT" parameters for your simulation,
>> I would consider it totally irrelevant to talk about "pressure" where your
>> constrain the volume.
>> This value or any other one has not really a meaning in this situation,
>> and I seen
>> many variations in the pressure value in this microcanonical ensemble
>> without paying too much
>> attention on it.
>>
>> In an "NPT" simulation, then you should be able to find back a normal 1
>> bar simulation I think.
>>
>> Do you have any reason to do first an NPT simulation, and "then" an NVT
>> one?
>> I would personally let the system equilibrate in NVT, then swith to the
>> more natural NPT,
>> provided actual code and force fields are now good enough in this ensemble.
>>
>> Best,
>>
>> Stéphane
>>
>> --
>> Team Protein Design In Silico
>> UFIP, UMR 6286 CNRS,
>> UFR Sciences et Techniques,
>> 2, rue de la Houssinière, Bât. 25,
>> 44322 Nantes cedex 03, France
>> Tél : +33 251 125 636
>> Fax : +33 251 125 632
>> http://www.ufip.univ-nantes.fr/ - http://www.steletch.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.