question about pressure in the reduce potential energy expresion

[matiasf86]

Hello, community,

I have a question regarding the correct value for the input pressure in the reduced potential expression when performing mbar or bar with the PyMBAR code for NPT trajectories. I believe my question will be better understood with an example.

Suppose I will be performing mbar for three k states (usually referred to as lambda in literature and molecular simulation codes). For each k state, I have frames from my molecular dynamics simulations.

k_1 x1 x2 x3 .....xn

k_2 y1 y2 y3 .....xn

k_3 z1 z2 z3 .....xn

In my example, x, y, and z correspond to my trajectory at each corresponding k or lambda state. Additionally, for each state, I have instantaneous volumes and pressures.

To perform mbar, I need to reevaluate the energy of my trajectories at the other states. This is usually done by changing the potential parameters to accommodate the other states. Consequently, I would obtain something like this:

k_1 E(x1) E(x2) E(x3) .....E(xn)
k_2 E(x1) E(x2) E(x3) .....E(xn)
k_3 E(x1) E(x2) E(x3) .....E(xn)

k_1 E(y1) E(y2) E(y3) .....E(yn)
k_2 E(y1) E(y2) E(y3) .....E(yn)
k_3 E(y1) E(y2) E(y3) .....E(yn)

k_1 E(z1) E(z2) E(z3).....E(zn)
k_2 E(z1) E(z2) E(z3) .....E(zn)
k_3 E(z1) E(z2) E(z3) .....E(zn)

Here E is the Energy potential or total, then it is called U

With this results I can construct the reduce portential matrix u_k(x) = beta_k [ U_k(x) + p_k V(x) + mu_k' n(x) ]

I understand that V(x) represents the instantaneous value of volume for each configuration obtained at the original different (k) states. However, I have a question regarding the value of pressure. Since pressure is computed from the virial, it will naturally change when evaluating a trajectory for a different k. For instance, I have my trajectory obtained at k_2

k_2 y1 y2 y3 .....xn
this has a pressure for each frame
k_2 P(y1) P(y2) P(y3) .....P(xn)

But then if I want to evaluate that trayectory at state k_1, the pressure will change as the Energy does.
k_1 P'(y1) P'(y2) P'(y3) .....P'(xn)

I find that this is not taken into account in MD programs, there is no easy way to retrive pressure from changing states. and ussually just the value PV(x) correspoding to the state from wich the trayectory was obtainded is informed.

Am I understanding this okay,?
When performing NPT simulations, a matrix of pressures should be input in the reduced potential energy matrix, in the same way we do with the energy ?.

Thanks for any help

Matias F.

[chapincavender]

In the reduced potential u_k(x) = beta_k [ U_k(x) + p_k V(x) + mu_k' n(x) ], quantities with a k subscript that are not explicitly functions of configuration x are part of the thermodynamic state and should not be estimated from the output of simulation trajectories. In particular, for NPT simulations, the pressure p_k should be set to the value used for the target pressure of the barostat and not computed from the instantaneous virial. In contrast, the volume V(x) should be computed from the simulation configurations. Similarly, the temperature used to calculate beta_k should be set to the value used for the target temperature of the thermostat/integrator and not computed from the instantaneous kinetic energies. For completeness, in a grand canonical simulation, the chemical potentials of each species mu_k should be set to the target chemical potentials that define the ensemble while the particle numbers n(x) should be calculated from the simulation configurations. Hope this clarifies things!

[mrshirts]

@chapincavender, thanks for jumping in! Yes, to be super clear, the quantities denoted "temperature" and "pressure" that are output by many simulations programs are NOT in fact the temperature and pressure of the NPT ensemble; they are stochastic quantities whose averages must equal the temperature and pressure. In the case of the temperature, it's the instantaneous kinetic energy scaled to be in units of temperature. Specifically < 2/3 K.E. / k_B N > = T, where K.E. is the kinetic energy, k_B is the Boltzmann constant, and N is the number of atoms. But it's only the AVERAGE that is equal to the temperature. NPT means in contact with a temperature bath - of constant temperature, and a pressure "bath" at constant pressure on the system.

[matiasf86]

Thank you for responding, now it's much clearer.