Stress tensor output for explicit volume dependent potentials

Hello everyone,

I am currently working with a system containing one atom type with the interactions between the atoms divided into bonded and non-bonded contributions. In the case of both contributions, the underlying potentials are not analytical but rather table potentials. In my Hamiltonian, I also happen to have an explicit volume dependent term, U_V(V). So my potential would look something as shown below, where U_R(R) denotes the sum of all the usually-existent bonded and non-bonded contributions that depend solely on relative spacing between my atoms (which in my case are tables) and U_V(V) is this “unusual” explicit volume dependent term.

image

It’s not a problem to address the dynamics in the NPT ensemble in such special condition (where I have an explicit volume dependent term) given that LAMMPS has a package (BOCS) that allows me to do so in the context of a modified version of the Nose-Hoover approach.

Naturally, like for many properties, the default compute (and output) pressure of LAMMPS would be missing the contribution related to delta(U_V(V))/deltaV. For example, below there is the version of pressure calculation that would be the “complete” one in a scenario where there is this explicit volume dependent addend in the Hamiltonian. The phi(r,V) would be my U_V(V).

image

Computing and outputting the “due pressure” is not so complicated given that I know the function U_V(V) and could manually add the missing term to the default pressure value by defining variables. Even easier, there is an option to already output the “due pressure” with a thermo_modify command specific of the BOCS package. Anyways, this is not the problem.

My problem is that now I would very much like to output the stress tensor. I have been trying to find (in a reliable yet effective (low time cost) way) a ready-to-go, complete formula for P_ij. My plan was to do the calculation of the the stress tensor by means of me defining variables as mentioned previously. However, what I am finding seems a little bit out of my league at the moment to be honest as the formulas are “a bit not straightforward” and so are their derivation :"D

Before grabbing my hat and coat and trying to hardcore’ly understand the paper and book I have found in order to do things properly, I thought of asking if there is some command within LAMMPS I could use to compute this in an easier/faster way where I can remain ignorant on the formula for P_ij under such condition but be sure that I am computing and outputting it correctly. It doesnt seem that there is one within the BOCS package… but maybe in some other package or some fancy thing that I dont yet know…

I will leave here how I solved the issue in case someone in the future finds himself/herself in the same situation. So, basically I didnt solve it within LAMMPS. Rather I figured out what was the proper term I needed to add to the pressure tensor LAMMPS standard’ly computes and added it to get the due pressure tensor in my case. Basically I used was a book by Mark Tuckerman (Statistical Mechanics: Theory and Molecular Simulation), mostly chapter 5 (and sometimes 4). The precise section in which the “complete” formula appears is 5.7. The derivation of the formula as a whole as shown in the book is quite fancy, requiring a background in Lagrangian mechanics (which I personally dont have) and there are some other equations that are not so trivial (to me at least). But even if you cannot follow the derivation of the final expression for the components of the stress tensor step-by-step and thus have to accept some things, by reading some different sections of the book, you can for sure be able to use the formula correctly and therefore safely compute the correct thing. It is not presented in a super complicated way, as I found to be when reading some other sources.

If you have some ability to code in c++, it is possible that you can maybe get LAMMPS to output this for you by modifying some source codes of the BOCS package. It seems that it even calculates P_xx, P_yy and P_zz already for secondary purposes.

1 Like