[gmx-users] Gromacs - Difference between target temperature and pressure in mdp and simulation

[Justin Lemkul]

On 11/25/19 2:15 PM, Fabio Bologna wrote:
> Dear Gromacs user-base and experts,
>
>       I need your help, as a Gromacs rookie. I've been trying to simulate the folding of a 20 residues mini-protein inside a dodecahedric explicit solvent box (30542 atoms) using Gromacs 2016.1. Because I need to start from the linear structure, my box is quite big ; I had to use the group cut-off scheme without explicit buffering instead of the default Verlet, in order to obtain acceptable performances. Also, I’m using the V-rescale thermostat (tau-t = 0.1) and Berendsen barostat (tau-p = 1.0). As in the title, my simulations don’t reach the intended target values of temperature and pressure. These are 300 K and 1 bar respectively.
>

Please be aware that the Berendsen thermostat and barostat produce 
incorrect ensembles and are generally never used for production 
simulations any more.

>
> Because I didn’t notice earlier these issues, I let 6 production md run for 69 ns. The averages temperatures ranged from 302 to 303 K and the average pressures ranged from 0.86 to 1.1 bar. In my humble opinion, I don’t think it is caused by a lack of equilibration, given the small size of the system and the long simulation time.
>
>
>
> I checked out the previous 1ns NVT and 1ns NPT equilibration phases (each using the very same thermostat and barostat, when applicable) of those 6 simulations. The average temperatures among them ranged from 300 to 302 K and average pressures from 0.018 bar to 1.4 bar. The temperature averages seems to have risen from equilibration to production MD (which completely puzzles me) while at least the pressure got closer to the target 1 bar (will it reach 1 bar later on, meaning that there are not enough points to make the mean value close to the target?)
>

Your pressure values are as good as any I've seen; it's a pressure that 
oscillates wildly (see Dallas' post for specifics on this). I've never 
seen a simulation with an exact 1.0 bar average.

>
> I ran again the NVT and NPT equilibration phases, experimenting with different values of tau-t and tau-p. The following results are the averages outputted by gmx energy (more accurate than manually averaging the values in the .xvg file you get from gmx energy itself, if I have understood correctly) at the end of the NPT equilibration.
>
>
> Regarding temperature, I got:
>
>    *   tau-t = 0.5 --> 310-315<tel:+39310315> K
>    *   tau-t = 0.1 ••• 300-302<tel:+39300302> K
>    *   tau-t = 0.01 --> basically only 300 K
>
> Shouldn’t tau-t just affect the amplitude of temperature oscillations? Clearly tau-t = 0.01 is the right value for my system, but as far as I know it is not used in productions MDs and only in equilibration phases. Is it still a

tau-t affects the response time of the thermostat, and I would argue 
that tau-t = 0.01 is *not* right for your system. The goal of the 
thermostat is not to *fix* the temperature at a given value, but to 
sample a correct, canonical kinetic energy distribution. If tau-t is too 
tight, you're getting the "right" temperature but the wrong ensemble.

> safe value for production or does it create artifacts of any kind? Also, why does the most used values of 0.5 and 0.1 not give reasonable results? I've test it with the Verlet scheme and the simulation reach the correct target temperature in a total of 2 ns. Also, why did the average temperature rise during production from 300-302 to 302-303 K? The thermostat setting are the same.
>
>
>
> Regarding the pressure, I got:
>
>    *   tau-p = 1.0 ••• 0.02-2 bar
>    *   tau-p = 0.5 ••• 0.6-3 bar
>
> It seems tau-p doesn’t have any effect…
>

The magnitude of the fluctuation is just as (if not more) important than 
the average. If you get 0.6 ± 5, then is your average in any way 
distinguishable from 1.0? What about 0.6 ± 50?

>
> Is there something I’m doing wrong or is it a normal behaviour? Especially regarding the pressure control, given that tau-t seems to work fine.
>
>
>
> Finally, on a different note, 1 of my 69 ns simulations crushed because of too many LINCS warnings. By checking the .err file, I found out that they were caused by a lysine and the N-terminus of the protein. In some frames the dihedrals of the protonated amine groups were a bit anomalous and the groups didn’t have a tetrahedral conformation. However the rest of the system was completeley fine. Also, they weren’t “misbehaving” at the same time. First the lysine had this issue but then the group returned to its correct conformation; then many ns later the N-terminus encountered the same problem and the total count of LINCS warning reached 1000. The crashed simulation was completely identical in parameters and chemical entities to the other 5 simulations and I didn’t see any spikes in the various energies of the system. Am I right to assume that it is only a numerical rounding error on the hardware part?
>
>
>
> I'm addition, I've attached the .mdp file I used in the production runs.

This mailing list does not accept attachments. If you wish to share a 
file, either copy and paste its contents or upload to a file-sharing 
service and provide a URL. I have never seen a stable simulation fail to 
achieve the correct temperature (and ensemble, when using a proper 
thermostat).

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

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