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)

where

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.