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