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] += Fx
f[i] += Fy
f[i] += 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.
Shenzhen University, China.