Pressure profile

Hi Lammps users,

I am trying to plot the pressure profile of the liquid (nonane) - vapor system at 400 K along the z-direction using "compute stress/atom " and “fix ave/chunk” commands with LAMMPS (17 Nov 2016). A part of the input file reads as,

compute T all temp
compute peratom all stress/atom T
compute ptensor all chunk/atom bin/1d z -93.5 1.0
fix 1 all ave/chunk 1 1 10 ptensor c_peratom[1] c_peratom[2] c_peratom[3] file press.dat

No. of atoms (ncount) and the per atom stress components (sxx, syy, szz) in the bins are recorded in the press.dat file.
In the post-processing, I have calculated the pressure components as Pii = - ncount*sii/(volume of bin)

The pressure profiles look reasonable in the vapor region. But in the liquid region, the pressure is relatively low ( pressure in the liquid is not equal to the pressure in vapor). The pressure in both liquid and vapor regions should be the same.

Can anyone please suggest me what I am doing wrong or missing during the pressure profile calculations? I could not figure it out. I have attached the plot.

I will appreciate any suggestions.

Pauf Neupane

Pressure.jpg

Hello Pauf, I believe it is reasonable to have regions of different partial pressure along a unit cell dimension. When you apply a barostat, only the corresponding total pressure (integral along that dimension) will be approach the target pressure of the barostat on average.

However, there is another more subtle problem with binning the results of compute stress/atom with fix ave/chunk, because it allocates atomic values of the stress only at the point positions of the atoms themselves. This means that there are discontinuities on length scales comparable to the interatomic potentials’ cutoffs: for a homogeneous system (i.e. bulk) these discontinuities will cancel out, but they will not in an interfacial system.

So this means that you can definitely use the liquid/vapor setup to compute e.g. the surface tension of the liquid, but it is quite difficult to analyze partial pressures profiles unless you use a large system, and bin the pressure at a resolution much coarser than the non-bonded cutoffs, which would still contain biases from how the stress is accumulated but much smaller in magnitude.

Giacomo