documentation and/or error handling for eam/alloy
demkowicz opened this issue · 12 comments
Summary
The documentation on pair_style eam/alloy does not specify how LAMMPS handles situations where the electron density at an atom site exceeds the maximum value specified by the potential file. The code doesn't crash or give an error. I suspect it applies the value of embedding energy at the maximum value of electron density specified in the potential file for all higher values of electron density. I suggest a) changing the code to give an error and exit or b) updating the documentation to make your convention explicit.
LAMMPS Version and Platform
23June2022
Details
I encountered this issue while computing binding energy curves (cold curves) for FCC Cu using the EAM potential of Mishin et al. (PRB 63, 224106 [2001]). The potential file that comes with LAMMPS (Cu_mishin1.eam.alloy) has a max electron density of 1.6402, which is exceeded at lattice parameters of about 3.21A and below. However, LAMMPS gives values of the cold curve for this potential down to lattice parameters as low as 2.5A (I didn't try any lower). When I substituted the embedding energy data in Cu_mishin1.eam.alloy with the analytical expression in PRB 63, 224106 [2001], I found that the cold curve differs from the one obtained using the potential file for atomic electron densities above 1.6402 (i.e., lattice parameters of about 3.21A and below).
I suspect it applies the value of embedding energy at the maximum value of electron density specified in the potential file for all higher values of electron density. I suggest a) changing the code to give an error and exit or b) updating the documentation to make your convention explicit.
Please note the following comment in the EAM pair style source code:
// if rho > rhomax (e.g. due to close approach of two atoms),
// will exceed table, so add linear term to conserve energy
[...]
if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);I don't think we should error out, but it may be worth documenting this and possibly print a (single?) warning if the table was exceeded.
Is fp(i) the embedding energy here? At rhomax? MJD
It is the derivative of the embedding energy (phi) at the location of atom i. Please have a look at:
https://github.com/lammps/lammps/blob/develop/src/MANYBODY/pair_eam.cpp#L226C3-L227C41
My question concerns cases where the local value of rho exceeds the range used to define the embedding energy. If the embedding energy is not defined, how then is the derivative of the embedding energy computed? At rhomax?
Alex, I read your comment the first time. Did you read mine?
Alex, I read your comment the first time. Did you read mine?
Yes.
So fp is obtained by extrapolating a spline fitted at rhomax? Trying to interpret the code in lines 234-239…
Is fp(i) the embedding energy here? At rhomax? MJD
rho[i] is the calculated embedding energy at atom I (due to surrounding atoms)
this if test is testing for the Q you are asking (what happens if actual rho exceeds limit of tabulated rho)
as the comment says, it is not an "error", but the code accounts for the possibility
by using the derivative (fp) at the end point to add extra energy to phi
If you think of the rho functional as a curve ending at the point rhomax, then for rho values beyond
that point, the curve is extended with a straight line with slope equal to the slope at the endpoint
So effectively, rho can have an arbitrarily large value. If that isn't what the user wants, they
need to use an EAM potential file which tabulates rho to larger values. Note that the pairwise
term in EAM typically supplies a large repulsive force as 2 atoms get close together. I.e. that
is not the job of the embedding term.
This could be explained in the doc page. I don't think any change to the code is needed.
Makes sense, Steve. I think the key is to let the user know that the max electron density has been exceeded, so perhaps a message from the code would do? The user can decide what do do after that.