Virial calculation with strain-based energy pair styles

Hello Lammps community!

I am currently working on my own Magnetoelastic pair potential within the SPIN package. I’m utilizing the a modified version of the magnetoelastic energy equation from Kittel (see here) and an atomistic strain tenson developed by Zimmerman (see here).

In my current implementation, an atom’s magnetoelastic energy is calculated through the atomic strain tensor via a displacement gradient of my neighboring atoms given an initial configuration at time t=0, and the orientation of the current atom’s spin vector direction.

In taking the derivative of this energy interaction with respect to displacement, I am left with only a single net force on my atom instead of the force between two atoms. I want this force to be added into the virial calculation so I can monitor its contribution to the stress arrays during energy minimization. What would be the best way to do this? The only method I saw in pair.cpp that does this was ev_tally_xyz_full but I don’t know if this is the ideal way to calculate its virial contribution.

I cannot answer this question, but there is the fix numdiff/virial command — LAMMPS documentation which can be used to check the results.

Hello Alex!

Thank you for the advice! I’ll look into this fix.

I guess if I were to ask my question in a more simple manner, is there any way to add a single net force per atom into the virial calculation? Or will I have to find a way to break down the net force into four smaller forces given four nearest atom neighbors and feed it into ev_tally_xyz ?

I cannot answer this, since this is out of my area of expertise. The way how the virial is computed in LAMMPS should be described in the following paper. Perhaps that has the information you seek.

Thompson, Plimpton, Mattson, J Chem Phys, 131, 154107 (2009).

Hello again!

Thank you again! I’ll jump in and take a look at this as well, thanks again!