Implementation of DFT-Calculated Force Constants as interatomic potentials in LAMMPS

Dear LAMMPS developers,

I am trying to extend LAMMPS pair_style to allow the use of DFT-Calculated harmonic and anharmonic force constants as interatomic potentials. The harmonic force-constants (3x3 matrices) are calculated using VASP+PHONOPY while the anharmonic force-constants (3x3x3 tensors) are calculated using VASP+ShengBTE packages, at the equilibrium geometry. Base on these force-constants, I have the following expression for force on an atom i;

F_i = - sum_j Phi_ij u_j - 1/2 sum_j sum_k Psi_ijk u_j u_k ----------- EQ(1)

Phi_ij = harmonic force-constants (3x3 matrices) between atom i and j (j=i allowed) and
Psi_ijk = anharmonic force-constants (3x3x3 tensors) among i, j and k (j=k=i allowed)
u_j = r_j(t) - r_j(eq) represents the displacement of atom j from its equilibrium position.

In the implementation, I am facing some problems. So I would like to ask two questions.

1- EQ(1) gives 3 components of force F on atom i due to displacement of atom i itself and its neighbouring atoms j and k,
which I add to f[i] in the compute() method;
f[i][0] += Fx
f[i][1] += Fy
f[i][2] += Fz

Do I need to add the same amount of F (but with opposite signs) to other atoms j (for harmonic part) and j and k (for anharmonic part) due to Newton’s third law?
If so, how to splite the contribution for j and k in the last term of EQ(1)?

My goal is to use MD based Green-Kubo method to calculated the thermal conductivity.
For that I have the following expression for heat current (with out convection)

J_i = sum_j ( r_i(t) - r_j(t) ) [ Phi_ij u_i + 1/2 sum_k Psi_ijk u_i u_k ] . v_j (t)

2- Given (r_i - r_j) and Phi_ij u_i + 1/2 sum_k Psi_ijk u_i u_k which method, i.e., ev_tally_xxxx(), should I call to calculated the per atom varial?

Any help will be appreciated.

Kind regards.

Zahid Rashid
Shenzhen University, China.

Not understanding all the details of this new pair style, but

for Q1, I think how you accumulate forces in atoms i,j,k will

depend on what kind of neighbor list your pair style requested.

If it is a half neighbor list and the newton flag is on (default),

then the IJ pair will appear only once in the collective neigh

list of all procs. So all forces that pair generates need to

be accumulated on all affected atoms. If you request a full

neighbor list, then IJ and JI will appear in the list, and you

can do something different.

For Q2, you can think of one interaction term in the potential

as producing force on some small cluster of atoms. The

correct virial for that cluster of N atoms is simply Sum(Fi x Ri) where

“x” is an outer product to produce 9 (or 6 symmetric) terms.

You can just accumulate 1/N of that sum to each of the N atoms.

The various ev_tally() functions do that for various N-body