The interpretation of pressure in LAMMPS for a non-ideal bonded system

Dear LAMMPS community,

Firstly, many thanks for this valuable mailing list, and for all the important answers that I have received from so many of you over the years.

This time, I would like to ask a question about pressure that I came across in while trying to understand how the enthalpy of a system is calculated in LAMMPS.

It is common knowledge that in MD, one evaluates pressure P as an ensemble of the instantaneous/microscopic pressure P*, that for a system of N particles in a volume V is as follows:

P*=1/V(1/3\sum_i m_i v_i •v_i+1/3\sum_i r_i•f_i)
where r is the position, v- velocity and f- force.

The macroscopic pressure is P=<P*>, a statistical average of an ensemble. In case of pairwise interactions, the pressure is defined as
P=<(N/V) k_BT>+<(1/3V)\sum_i sum_{j<i} r_{ij}•f_{ij}

Where r_{ij} is the intermolecular vector between i and j and f the corresponding force.

When we use NVE+Langevin (mimicking Brownian Dynamics), a simulation of an ensemble of ideal gas particles produces a linear relation between pressure and temperature upon cooling, as expected.

However, what happens if we have the same NVE+Langevin dynamics whereby the monomers are connected through a bond, say harmonic or FENE? What if one has a collapsing polymer chain (gradually with a specific rate of cooling), where pressure is not the same as it is in the ideal gas case due to bonding interactions? If one wishes to study the behavior of enthalpy as a function of Temperature, for instance, how does one define the pressure that is part of the enthalpy calculation?

Many thanks in advance! Your input is highly valued & appreciated.

With best wishes
Anna

Dear LAMMPS community,

Firstly, many thanks for this valuable mailing list, and for all the
important answers that I have received from so many of you over the years.

This time, I would like to ask a question about pressure that I came
across in while trying to understand how the enthalpy of a system is
calculated in LAMMPS.

It is common knowledge that in MD, one evaluates pressure P as an ensemble
of the instantaneous/microscopic pressure P*, that for a system of N
particles in a volume V is as follows:

P*=1/V(1/3\sum_i m_i v_i •v_i+1/3\sum_i r_i•f_i)
where r is the position, v- velocity and f- force.

The macroscopic pressure is P=<P*>, a statistical average of an ensemble.
In case of pairwise interactions, the pressure is defined as
P=<(N/V) k_BT>+<(1/3V)\sum_i sum_{j<i} r_{ij}•f_{ij}

Where r_{ij} is the intermolecular vector between i and j and f the
corresponding force.

When we use NVE+Langevin (mimicking Brownian Dynamics), a simulation of an
ensemble of ideal gas particles produces a linear relation between pressure
and temperature upon cooling, as expected.

​this last statement is not correct. particles in an ideal gas do not
interact, so there are no forces. you only have the kinetic energy
contribution and if you look at it closely, you will see that you can
recover the ideal gas law from it.​ also the interactions implicitly
contained in fix langevin do not apply to an ideal gas.

​the F <dot> R term only applies to interacting particles.

However, what happens if we have the same NVE+Langevin dynamics whereby
the monomers are connected through a bond, say harmonic or FENE? What if
one has a collapsing polymer chain (gradually with a specific rate of
cooling), where pressure is not the same as it is in the ideal gas case due
to bonding interactions? If one wishes to study the behavior of enthalpy as
a function of Temperature, for instance, how does one define the pressure
that is part of the enthalpy calculation?

​the F <dot> R relation can be applied for bonded interactions just as
well.​ it is easy to set up tests for that.
i've done this for example last fall to validate two bugfixes:
https://github.com/lammps/lammps/pull/213

​axel.​

Thank you Axel,

Dear Axel

What remains a confusion, that all these terms are located where the polymer chain is - are different on “droplet” -the collapsed polymer- surface, and certainly zero outside of “droplet”. Whereas thermodynamic pressure must be constant across the system in equilibrium. Am I interpreting this correctly?

Many thanks
Anna

Dear Axel

What remains a confusion, that all these terms are located where the
polymer chain is - are different on "droplet" -the collapsed polymer-
surface, and certainly zero outside of "droplet". Whereas thermodynamic
pressure must be constant across the system in equilibrium. Am I
interpreting this correctly?

​i don't think so.

pressure is a property of the entire observed system and the physical
interpretation of (total) pressure is (the average) force per area on the
bounding surface(s). thus the pressure of an isolated system is by
definition zero, since​ the volume is infinite, i.e. unbounded.

now, you *can* compute a "local" pressure, by subdividing the volume and
computing a pressure for this, but there is nothing requiring a system to
have equipartitioning on this kind of property in an inhomogeneous system.
all that is required would be an (on average) zero net pressure on the
dividing surfaces between two such subsystems.

axel.

I should also point out that if you are using the Langevin thermostat to model the presence of an implicit solvent, then the pressure reported by LAMMPS (based on velocities and forces of the solute “atoms”) is not capturing the dominant contribution to the pressure from the solvent. As an example, you could simulate polymer chains dissolved in water compressed to very high pressure using this approach, but the pressure reported by LAMMPS would not reflect the high pressure of the water molecules. The LAMMPS pressure is probably related to the osmotic pressure of the polymer.

Aidan