Pressure calculation from the stress tensor

hghcomphys commented 4 years ago

Hi,

I noticed that for an example of the bulk liquid water MD simulation at ambient conditions there is a discrepancy between the pressure that is evaluated directly from LAMMPS Press variable and the from diagonal components of the stress tensor.

For any given potential, such as SPC/E, both pressures are expected to be exactly the same at any given time. However, as you can see from the part of the 'log.lammps' file, the pressures ('Press' vs. v_press) are different:

Per MPI rank memory allocation (min/avg/max) = 7.872 | 7.873 | 7.873 Mbytes Step Temp TotEng Press v_press v_dens 10000 300 -24805.565 16173.004 4116.885 1.0088118 10010 309.72891 -24805.818 -6442.4868 4249.4215 1.0085808 10020 370.76634 -24806.738 13054.728 5085.9069 1.0083952 10030 368.9284 -24808.713 -10678.673 5057.6896 1.0077962 10040 368.35577 -24812.762 14130.233 5048.4766 1.0075243 10050 324.97034 -24817.324 -8696.6602 4450.5506 1.0067756 10060 344.73275 -24822.403 11243.846 4718.753 1.0062534 10070 350.02017 -24827.523 -5208.2374 4786.6142 1.0053054 10080 355.17204 -24832.361 5075.6256 4853.0034 1.0044642 10090 344.0271 -24836.953 135.99129 4695.8461 1.0034225
...

lammps version: 20 Apr 2020 n2p2: the most recent version nnp potential parameters: I used mine but it can be also this pressure calculation from stress tensor:

compute     Tw all temp
compute     peratom all stress/atom Tw
compute     sxx all   reduce sum c_peratom[1]
compute     syy all   reduce sum c_peratom[2]
compute     szz all   reduce sum c_peratom[3]
variable    press equal   -(c_sxx+c_syy+c_szz)/(3*vol)  # pressure

Best, Hossein

singraber commented 4 years ago

Hey,

hmm, thanks for notifying me, that is indeed confusing.. I could reproduce the problem and found that the correct results were returned when I comment out the line https://github.com/CompPhysVienna/n2p2/blob/c7b4407317fdf9bfcfd5a20f805f5e8948083243/src/interface/LAMMPS/src/USER-NNP/pair_nnp.cpp#L70 I am not sure why this happens, but I will investigate further.. stay tuned..

singraber commented 4 years ago

I just realized that this topic was brought up previously, see #10. The per-atom stress computation is not implemented and unfortunately I forgot to throw an error if it is requested by the user. The pressure on the other hand should be correct. Do you need the stress computation for your work? I agree that this would certainly be a good test to see if everything works correctly and I am beginning to think it may be worth the effort to implement it...

hghcomphys commented 4 years ago

Hi Andi,

Thanks for the swift reply.

Based on my experience, I can tell you that the per-atom stress tensor is quite useful when there is a need to calculated applied pressure on a specific group of atoms. For instance, as a confined system of water between membranes (e.g. 2D channels, bubbles, etc), it can be used to evaluate applied pressure on the water. Furthermore, stress tensor can be used for calculating several dynamical coefficients such as bulk and shear viscosity components within the framework of Green-Kubo formalism, see this for more details.

There is no strict need for my work at the moment. But, I will definitely need it for future works.

Due to the fact that NNP is a many-body potential and can not be express in the summations of few-body terms, it admits that may not be straight forward to implement it. However, it is important that like many other many-body potentials such as ReaxFF and Airebo, NNP can properly produce the per-atom stress tensor components.

Best, Hossein

singraber commented 4 years ago

Hi Hossein,

thanks for your reply and the explanations, I totally agree with you that this functionality would be useful for many users. I am putting this on the ToDo list... I will report back if there is something ready for testing.

Best, Andi

hghcomphys commented 4 years ago

Hi Andi,

These days I am free and quite interested in the development of the N2P2 base code concerning the per-atom stress-tensor calculation. If you think that I can contribute with this part, please let me know. We can have a further discussion through email "hgh.comphys@gmail.com".

Best, Hossein

singraber commented 4 years ago

Hi Hossein,

thanks again for offering your help with this issue! To get things started I will put some remarks here:

For correctly implementing this feature I think it is best to first have a look how LAMMPS computes these quantities. The essential part of the documentation says:

Details of how compute stress/atom obtains the virial for individual atoms for either pairwise or many-body potentials, and including the effects of periodic boundary conditions is discussed in (Thompson). The basic idea for many-body potentials is to treat each component of the force computation between a small cluster of atoms in the same manner as in the formula above for bond, angle, dihedral, etc interactions. Namely the quantity r * F is summed over the atoms in the interaction, with the vectors unwrapped by periodic boundaries so that the cluster of atoms is close together. The total contribution for the cluster interaction is divided evenly among those atoms.

where they cite this paper:

Thompson, A. P.; Plimpton, S. J.; Mattson, W. General Formulation of Pressure and Stress Tensor for Arbitrary Many-Body Interaction Potentials under Periodic Boundary Conditions. J. Chem. Phys. 2009, 131 (15), 154107. https://doi.org/10.1063/1.3245303.

There, in section III are basically the instructions of how to calculate the virial for many-body potentials under PBC. I found it is definitely not an easy read... I will have to go through it once again to fully grasp how this works.
On the other end in the LAMMPS and n2p2 code the relevant parts are to call correctly the function ev_tally() or maybe rather ev_tally_xyz() with the right arguments in PairNNP::compute(). This will probably be very similar to the force summation over local and ghost atoms in InterfaceLammps::getForces().

hghcomphys commented 4 years ago

Hi Andi,

Thanks for the reply and remarks. To start out, I first need to read the paper and have a look at LAMMPS and N2P2 relevant code bases. This paper suggested by Prof. Behler is also useful on how to calculate the per-atom stress tensor for NNP.

For the bulk water example, I still don't understand where LAMMPs calculates the right pressure values when we comment below the line in the pair_nnp.cpp.

if (vflag_fdotr) virial_fdotr_compute();

It may be clear when we go through the details.

I will let you know about my progress...

Best, Hossein

songzhigong commented 4 years ago

Hi,

 I reported this issue of nnp before, I have tested some ML potential package. Do you think it is feasible to integrate the system virial stress as into the fitting framework, just like GAP?

Best

singraber commented 3 years ago

Just a bit of information...

It seems that computing the per-atom stress for manybody potentials is not trivial, in particular in the context of thermal conductivity calculations. There were even examples of broken code in LAMMPS for three-/fourbody potentials, see these two publications:

(1) Surblys, D.; Matsubara, H.; Kikugawa, G.; Ohara, T. Application of Atomic Stress to Compute Heat Flux via Molecular Dynamics for Systems with Many-Body Interactions. Phys. Rev. E 2019, 99 (5), 051301. https://doi.org/10.1103/PhysRevE.99.051301.

(2) Boone, P.; Babaei, H.; Wilmer, C. E. Heat Flux for Many-Body Interactions: Corrections to LAMMPS. J. Chem. Theory Comput. 2019, 15 (10), 5579–5587. https://doi.org/10.1021/acs.jctc.9b00252.

The method presented in (1) is already implemented in LAMMPS, see here: https://github.com/lammps/lammps/pull/1704

CompPhysVienna / n2p2

Pressure calculation from the stress tensor #42