Compute local/pair for EAM potential

Dear LAMMPS developers,

I would like to compute the local pairwise distances and forces for EAM potential, as I want to calculate the local pressure using Irving-Kirkood’s definition.

I’ve searched the LAMMPS mailing list, and what I found is that Steve mentioned that the force is only computed for the pairwise-potential part of the EAM (http://lammps.sandia.gov/threads/msg13125.html). But this is a very old post, and I do not see any restriction of using the pair/local in the Doc page. I wonder if the compute local/pair will work for EAM potential now. I have simply tested one case by substituting the pairwise distance into the pairwise potential part only to calculate the force. I didn’t get the same force computed by LAMMPS. However, I can get the same results when only using the pairwise potential function. So I guess “compute pair/local” also considers the density part of EAM. Could you please tell me if the LAMMPS can compute the local pair forces for EAM?

Thank you for your time.

Regards,
Pengyu

Nothing has changed in LAMMPS in this respect.
Pair local computes the energy/force between a pair
of atoms. It works as you expect for pairwise pair styles,
but for a many-body potential like EAM, there is no
simple expression for the total force between a pair
of atoms. The embedding term is a non-linear
functional of a quantity that is a sum over all an
atom’s neighbor’s contributions to the density.

Steve

Dear LAMMPS developers,

I would like to compute the local pairwise distances and forces for EAM
potential, as I want to calculate the local pressure using Irving-Kirkood's
definition.

I've searched the LAMMPS mailing list, and what I found is that Steve
mentioned that the force is only computed for the pairwise-potential part
of the EAM (http://lammps.sandia.gov/threads/msg13125.html). But this is
a very old post, and I do not see any restriction of using the pair/local
in the Doc page. I wonder if the compute local/pair will work for EAM
potential now. I have simply tested one case by substituting the pairwise
distance into the pairwise potential part only to calculate the force. I
didn't get the same force computed by LAMMPS. However, I can get the same
results when only using the pairwise potential function. So I guess
"compute pair/local" also considers the density part of EAM. Could you
please tell me if the LAMMPS can compute the local pair forces for EAM?

​you may be able to bypass this entire issue, by writing your pressure​
computation as a USER-TALLY compute.
please have a look at code in src/USER-TALLY/compute_stress_tally.cpp and
src/USER-TALLY/compute_stress_tally.h which computes the per atom virial
stress from pairwise forces, similar to compute stress/atom, but limited to
the interactions between two groups. this is done by inserting a callback
function into the force+energy tally that is called for each interacting
pair *during* the regular pair force computation. this will work for EAM
even though it is a manybody potential, since it uses the same tally
function as pairwise additive potentials.

dumping pairwise information to a file and post-processing it, can become a
troublesome issue, since those files will become very large very rapidly
and thus accumulating the information directly from inside the code may
save you a lot of trouble and effort.

​axel.​

Dear Steve and Axel,

Thank you for your reply and suggestions.

I would like to provide more details about my simulation and updates of what I found from using “compute pair/local”. I actually use the pair_style eam/alloy to implement manybody DPD interactions, with a tabulated file. In manybody DPD, the force and energy have an additive term and a density-dependent term. The force between atom i and j can be expressed in terms of distance between them, the local density contribution to atom i and the local density contribution to atom j. Hence we can write the total energy on atom i as a function like EAM potentials. From the results of “compute pair/local”, the calculated pair energy is just equal to the additive term, while the calculated pair force is actually equal to the value which I calculated using the analytical force function (as now I know how many atoms around atom i and j respectively). I think that the calculation of force from pair/local is correct in this case, probably because force between two atoms has a “simple” expression. Would you please common on this result?

I agree with Axel that dumping pairwise information will generate lots of data. I would try to do this computation inside the LAMMPS if possible. I had a look at the codes you mentioned and callback function and I think they will be useful. However, I never code in C++ before so I guess it would be hard for me.

Thank you again for your help.

Regards,
Pengyu

You can look at the single() function in pair_eam.cpp and
see what it returns. It is this:

The force between atom i and j can be expressed in terms of distance between them, the local density contribution to atom i and the local density contribution to atom j.

I.e. the pair term and the density dependent term. But the
latter is not just dependent on i,j but on all the atoms that
contributed to density_i and density_j.

Steve

Dear Steve and Axel,

Thank you for your reply and suggestions.

I would like to provide more details about my simulation and updates of
what I found from using "compute pair/local". I actually use the pair_style
eam/alloy to implement manybody DPD interactions, with a tabulated file. In
manybody DPD, the force and energy have an additive term and a
density-dependent term. The force between atom i and j can be expressed in
terms of distance between them, the local density contribution to atom i
and the local density contribution to atom j. Hence we can write the total
energy on atom i as a function like EAM potentials. From the results of
"compute pair/local", the calculated pair energy is just equal to the
additive term, while the calculated pair force is actually equal to the
value which I calculated using the analytical force function (as now I know
how many atoms around atom i and j respectively). I think that the
calculation of force from pair/local is correct in this case, probably
because force between two atoms has a "simple" expression. Would you please
common on this result?

from looking at the source, ​the single function in the EAM potential is a
very special case. ​in general, many-body potentials cannot be represented
by pairwise additive terms. as steve mentioned, for EAM, there is the
derivative of the embedding energy which includes contributions from all
neighbors. however, this contribution is collected into a per-atom term and
then used to compute a pairwise force and this per-atom term is cached in
the internal storage for the EAM potential. hence the single function will
reproduce the exact pairwise force and energy contribution (indeed the
expressions for the computation of energy and force are the same modulus
the scaling factor for use with fix adapt), but *only* if the atom
positions are unchanged.

I agree with Axel that dumping pairwise information will generate lots of
data. I would try to do this computation inside the LAMMPS if possible. I
had a look at the codes you mentioned and callback function and I think
they will be useful. However, I never code in C++ before so I guess it
would be hard for me.

​i think that is a very shortsighted view. if you are a good programmer,
then it should not take you much time to teach yourself sufficient C++.
keep in mind, that you already have a template for most important steps and
only need to fill in the physics. i grant you that this is a bit more
complex due to the way the stress term is defined, that you want to
implement, but it is not insurmountable. and you may get help from the list
here. for a standalone code, you have to write the entire processing from
scratch, you have to write a parser for the dumped data, it will not be
parallel and there is a massive overhead and slowdown due to the reading,
storing and writing of files. in my book, the time and effort to enable you
to implement a cleaner and better solution is small compared to the
long-term disadvantages and overhead of doing this externally.

out of curiosity, what programming language were you planning to use for
your post-processing software?

​axel.​

Dear Steve and Axel,

Thanks again for all useful information.

After reading the pair_eam.cpp file, the pair force calculated from single() function is the desired value. I’ve also looked at the ev_tally() function, it seems like the fpair is the same pair force that I want. Therefore, I will try to make use of it and see if I can do the computation in LAMMPS.

To answer your question, Axel. I have actually written the simple codes in Matlab. Since the data is large, I have a function to load a portion of the dumped data at a time (and save the data into binary files for later use), and another function to process this portion of data and then accumulate the results. As you expect, it is slow and most of the time was spent on loading the data. So I will try to learn how to code in c++ and implement it in LAMMPS. After all, it is a good practice and I can learn more about LAMMPS.

Thanks and regards,

Pengyu