LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Implementation of DFT-Calculated Force Constants as interatomic potentials in LAMMPS
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Implementation of DFT-Calculated Force Constants as interatomic potentials in LAMMPS

From: "Zahid Rashid" <Z.Rashid@...2403...>
Date: Tue, 17 Apr 2018 16:09:04 +0200

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.