[gmx-users] Determining density of a gas in a simulation

[Mark Abraham]

Darrell Koskinen wrote:
> Hi Mark,
> I took a look at a tutorial where they perform an NPT simulation 
> (http://davapc1.bioch.dundee.ac.uk/prodrg/gmx.pdf, page 11) and it 
> appears to me that they simply ran the simulation with pressure and 
> tempurature coupling "on" for a sufficient amount of time to reach 
> equilibrium.
> I also searched in the GROMACS manual for NPT and came across the 
> following:
> /"The generally more useful Gibbs free energy G is related to the 
> partition function
> 
 of an N, p, T ensemble, which is assumed to be the equilibrium 
> ensemble generated by a MD
> simulation at constant pressure and temperature"./
> This appears to indicate that in order to run an NPT simulation, we 
> simply need to maintain constant temperature and pressure, which I 
> assume simply means letting the system reach equilibrium, since, at 
> equilibrium, the temnperature and pressure will be reasonably constant.

Sure... that's what equilibrium means - a stable point to which a system 
will return if perturbed. The starting configuration is a perturbation. 
In principle, you could use the thermostat and barostat until you 
achieved equilibrium, and then either turn them off so the dynamics are 
unperturbed, or continue using such coupling algorithms that do not 
affect the ensemble. Berendsen, notably, does perturb the ensemble 
despite its widespread use. In practice, you want to keep the coupling 
algorithm in cases where the integration errors are appreciable over the 
length of the simulation, etc.

Mark

> Darrell
> Date: Fri, 06 Nov 2009 19:50:48 +1100
> From: Mark Abraham <Mark.Abraham at anu.edu.au>
> Subject: Re: [gmx-users] Determining density of a gas in a simulation
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <4AF3E368.3060304 at anu.edu.au>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> Darrell Koskinen wrote:
> 
>> > Dear GROMACS Gurus,
>> > I wanted to determine the density of the ammonia gas in my 
>> simulation > and compare this to the reference values and thought I 
>> could use > g_density for this purpose. I recall being told that I 
>> need to perform a > NPT simulation and achieve equilibrium. How do I 
>> perform an NPT > simulation? I assume this means that I fix the number 
>> of molecules, the > pressure, and the temperature. The number of 
>> molecules is definitely > fixed and I define the temperature in my 
>> .mdp file through > "ref_t           =240    240". Is this the way in 
>> which to set the > temperature in an NPT simulation? I do not 
>> currently specify any > pressure in my simulation. Do I need to set 
>> the pressure as well through > a ref_p setting? And does the ref_p 
>> setting simply need to satisfy the > ideal gas equation for the 
>> specified N, T, & V?
>>   
> 
> P and T can't be "set". You can have an ensemble whose average over a 
> long period has a value. You can use an algorithm to move them closer to 
> a desired equilibrium value if the instantaneous one differs. See the 
> manual for more information. Doing some tutorial material and some 
> background reading sounds like a good idea, too.
> 
> 
>> > And how do I determine when equilibrium is achieved? Should the > 
>> potential energy of the system or some other value stabilize to within 
>> > some deviation of the absolute value? If so, what would would be the 
>> > acceptable limit of this deviation ... 5% of the value? Other 
>> percentage?
>>   
> 
> It depends - there's no hard-and-fast rule and it will vary with which 
> observables you intend to measure.
> 
> Mark
> 
> 
>