Meam_force questions

Questions concerning the meam_force.f file:
1. What information does fmap hold? It just depends on the atom type, correct?
2. Line 109: What information do scrfcn and fcpair store?
3. What are eflag_either, eflag_global, eflag_atom and how are they set?
4. Line 121: How is rdrar set? What is it?

Greg Wager (CCd) would have to answer that level of detail. He is the author of the code.


Hi Anne,

I can answer these questions. In general, it’s helpful to look at mean_force.F in the context of where it gets called in pair_meam.cpp (under the MEAM directory in the main LAMMPS source). A lot of the quantities you’re asking about get set either by LAMMPS itself, or by other calls to the meam library (like meam_dens_init and meam_dens_final) from the same file. Also, you’ll probably want to refer to the Gullett et al. report referenced on the meam documentation pages on the LAMMPS web site.

To your specific questions:

  1. You probably know that LAMMPS allows atoms to be grouped as “types”. In the input file line command that defines the meam potential allows you to assign an element (e.g. Si or C) to each type. See These elements in turn get enumerated in the meam parameter file you provide. So, fmap contains a map from the atom type to the particular element, so that fmap(itype) = ielement. So, if types 1, 2, and 3 are silicon atoms and type 4 are carbon atoms, then fmap = [1,1,1,2] where Si=1 and C=2.

  2. The array scrfcn contains the screen function values defined as \bar{S}_ij in equation 4.11b-4.11e in the Gullett report. Note that scrfcn is a flattened array that contains the values for all i-j pairs in the system, and it gets passed from pair_meam.cpp to include an offset that indexes into the part of the array containing the screen function values for all of the neighbors (j) of atom i. Similarly, fcpair stores the cutoff functions for the i-j pairs, f_c in equation 4.11a. The product you see computed on line 109 of meam_force.F corresponds to equation 4.11a. Both scrfcn and fcpair get computed for the i-j pairs in the meam_dens_init() subroutine called from pair_meam.cpp.

  3. These are flags set by LAMMPS itself: eflag_global = 1 means compute the total potential energy of the system, eflag_atom = 1 means compute the energy of each individual atom., and eflag_either = 1 if either one of these are set. This is set in the main LAMMPS code itself (see for example ev_setup() in integrate.cpp), but at a glance I don’t see exactly where this ultimately comes from.

  4. The pair potential phi® is tabulated as a function of r at the beginning of the run for a discrete set of values of r, separated by increment dr. During the computation of the energy and force, we use a cubic spline interpolation of these tabulated values, rather than computing the value of phi® from scratch. This also allows us to easily compute the derivative phi’®. This spline interpolation is what’s going on on lines 120-130 in meam_force.F. We hard-code the number of discrete values in this interpolation to 1000, and compute dr so that we span the maximum cutoff radius. Finally, rdrar is set in meam_setup_done.F to be the reciprocal of dr. Line 130 of meam_force.F is simply using rdrar to compute appropriate index (pp) into the tabulated values.

Hope that helps!


Greg Wagner
Manager, Thermal/Fluid Science and Engineering Department
Sandia National Laboratories, Livermore, CA
Tel: (925) 294-2180 Fax: (925) 294-3410
Email: [email protected]…3…