LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Implementation of DFT-Calculated Force Constants as interatomic potentials in LAMMPS

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

 From: Steve Plimpton Date: Tue, 17 Apr 2018 08:58:23 -0600

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
terms.

Steve

On Tue, Apr 17, 2018 at 8:09 AM, Zahid Rashid wrote:
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.

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most