Compute stress / atom - Group form expression (PBC)

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I am trying to program the calculation of per atom virial stress tensor for the pair potential as mentioned in the compute stress/atom documentation. For the periodic boundary condition case, it has been mentioned in (Thompson) Thompson, Plimpton, Mattson, J Chem Phys, 131, 154107 (2009), that the group form expression is the suitable one.

  1. But I am unclear about how the per atom virial contribution can be extracted from the group summation?

  2. Whether the summation of rij dot Fij is done only with the atoms in a particular group or with the atoms within its own cut off sphere for a reference atom in that group?

Thanks in advance !

Best Regards,
Vimal Muthusamy

2

Dear LAMMPS users,

I am trying to program the calculation of per atom virial stress tensor

program where and how?

for the pair potential

for what pair potential?

as mentioned in the compute stress/atom documentation. For the periodic boundary condition case, it has been mentioned in (Thompson) Thompson, Plimpton, Mattson, J Chem Phys, 131, 154107 (2009), that the group form expression is the suitable one.

  1. But I am unclear about how the per atom virial contribution can be extracted from the group summation?

virial contributions are always computed by looking at contributions from the closest periodic copies of each pair or group of atoms (for interactions containing 3 or more atoms).

  1. Whether the summation of rij dot Fij is done only with the atoms in a particular group or with the atoms within its own cut off sphere for a reference atom in that group?

contributions are either tallied individually as part of the pair styles (e.g. in manybody potentials), or using F dot r before doing the reverse communication of the force contributions back to the local atoms (in case of newton pair being enabled). for details, you have to look at the source code of the individual pair styles. you can see flags in the constructor whether the virial may be computed using F dot r (e.g. not possible for TIP4P styles due to geometry issues) and calls to various “tally” functions that increment per atom stress, if requested.

axel.

3

Dear Axel Kohlmeyer,

Thank you so much for your quick response.

I am doing it for the EAM potential. For the non-periodic case, using F dot r for the atoms within the cut-off sphere works fine. But for the periodic case, even after including the contributions from periodic copies, the stress is non-zero but uniform for the case of a single homogeneous crystal at 0 K.

Best regards,
Vimal

4

you have to provide more details of what it is exactly what you are doing (and programming and programming how) and what you are comparing to.

why should the stress be zero in the case you describe? it doesn’t have to be in general? it would only be in case you have a crystal as exactly the lattice constant that is consistent with the potential in use.

axel.

5

Dear Axel Kohlmeyer,

I apologize for that. My routine calculates the interatomic forces for the EAM potential by using atomic coordinates from the LAMMPS dump file. The outputs (resultant force per atom & pair-wise interactions for all atoms) from the force routine has been verified with the results from LAMMPS for the same set of inputs (Potential and boundary conditions).

Now, using the results from the force routine, stress per atom is calculated in a separate routine. The stress calculated deviates from the stress results (which is zero) in LAMMPS for the above stated conditions.

In LAMMPS, the sample has been completely relaxed using cg style and convergence was achieved for the force tolerance of 1e-8.

Best Regards,
Vimal Muthusamy

6

Dear Axel Kohlmeyer,

I apologize for that. My routine calculates the interatomic forces for the EAM potential by using atomic coordinates from the LAMMPS dump file. The outputs (resultant force per atom & pair-wise interactions for all atoms) from the force routine has been verified with the results from LAMMPS for the same set of inputs (Potential and boundary conditions).

Now, using the results from the force routine, stress per atom is calculated in a separate routine. The stress calculated deviates from the stress results (which is zero) in LAMMPS for the above stated conditions.

In LAMMPS, the sample has been completely relaxed using cg style and convergence was achieved for the force tolerance of 1e-8.

Best Regards,
Vimal Muthusamy

7

so what you are saying is, that LAMMPS does compute things correctly, but your code doesn’t. that is too bad, but ultimately your problem. i have already explained in broad strokes of what LAMMPS is doing. it is your job to figure out where you are deviating. my guess would be, that you are not using the correct value of “r” for F dot r when using newton’s 3rd law to compute forces.

axel.

8

Dear Axel Kohlmeyer,

Thank you very much for the clarification.

Best regards,
Vimal