Dear LAMMPS forum,
I would appreciate an explanations on inclusion of ZBL potential, using the flag zbl=1 for the MEAM potential.
As far, as I could see, the flag zbl(I,J) = 1 allows intrinsically to blend the MEAM I-J pair potential with the ZBL potential for small atom separations. Well, it looks simple but how does LAMMPS “know” what the “small atom separations” are? I use B.-J. Lee, J.-H. Shim, and M.I. Baskes (2003), Physical Review B, 68(14), 144112 MEAM version and did not find any indication of the inner and outer distances for MEAM to ZBL switching function in the MEAM data files “library.meam or/and Al.meam”.
The only explanation that I’ve found in one of the papers is : “The MEAM potentials are splined to a ZBL potential in the standardized method implemented in to LAMMPS”
Thank you very much
So to figure things like this out, you typically have to contact the authors for a detailed explanation, but that may be difficult with a publication that is 20 years old, since people move on (especially grad students and postdocs, who do the bulk of this kind of research work) and then it is difficult to remember or even to get hold of them.
The next best option is to look at the source code and reverse engineer it. Since you put this into the beginner category, I will give a short demonstration (otherwise, I would just have pointed out the file).
Mind you, whatever I am showing here, I just reverse engineered myself after reading your post. And while I am somewhat familiar with the LAMMPS code and what it looks like, I have only a very minimal understanding of the MEAM model. So take everything I write with a grain of salt.
To start reverse engineering, you need a “hook”, i.e. something to look for that is a known entity. In this case, the “hook” would be to find the section with the ZBL function is computed. If you’re lucky, it will be signaled by the function name, but ZBL can also be recognized by the numerical constants used in the potential function (see: pair_style zbl command — LAMMPS documentation). So if the former is the case, you can confirm it by the latter. In this case there is a function MEAM::zbl() that looks like a good candidate and indeed it matches. Looking at the function it is evident that this computes the “full” ZBL potential using the distance and core charges of the two atoms in question.
The next step is now to look for places, where this function is called, and - fortunately - there is just one location and that is somewhat down in the MEAM::compute_pair_meam() function. There we have the following block of code that is executed if the zbl flag is set for a pair of atom types ( if (zbl_meam[a][b] == 1) {).
// For Zbl potential:
// if astar <= -3
// potential is zbl potential
// else if -3 < astar < -1
// potential is linear combination with zbl potential
// endif
if (zbl_meam[a][b] == 1) {
astar = alpha_meam[a][b] * (r / re_meam[a][b] - 1.0);
if (astar <= -3.0)
phir[nv2][j] = zbl(r, ielt_meam[a], ielt_meam[b]);
else if (astar > -3.0 && astar < -1.0) {
frac = fcut(1 - (astar + 1.0) / (-3.0 + 1.0));
phizbl = zbl(r, ielt_meam[a], ielt_meam[b]);
phir[nv2][j] = frac * phir[nv2][j] + (1 - frac) * phizbl;
}
}
So this first computes the parameter A^*(r) which takes the \alpha parameter and multiplies it with a scaling factor that starts from -1 at r=0, is 0 at r=r_e (which is the closest equilibrium distance between two atoms computed from the lattice constant and the lattice geometry). Looking this up in the library.meam file for, say fcc-copper we get: \alpha = 5.10570729 and r_e = 3.62 \cdot \sqrt{2} = 5.119453. So for very close distances (A^* <= -3) the full ZBL function is applied and then a linear interpolation while ($A^* < -1). Beyond that, the regular MEAM pairwise interaction is applied. So for the copper example you have full ZBL until about 2.11 Angstrom and then linear switching until about 4.12 Angstrom. Due to the use of the two material specific parameters that are related to the size and electron density of the elements, the switching is material specific.
Dear Axel, Thank you very much for your exact and exhaustively complete answer to my question. Also, I am very sorry that I forced you to do the work that I should have done. Of course, you’ve found the answer much faster and more precisely than I would done. Thanks again and regards, Victor
Don’t worry. I didn’t need to do it. I did it, because it was a nice opportunity and a very good example to showcase the process. Since the forum is archived, this will hopefully be found by others in the future that are looking for ways to figure out something through the same procedure.
Unlike other people that have trouble following the conventions for such a forum, you did three very important steps correctly: 1) you did research on your own before posting and documented that, 2) you described your problem concise and sufficiently complete, 3) you flagged your post as a “beginner level” topic, thus my “beginner level” response.
I am always happy to spend time and responding to such topics, because it can be a positive example for how this forum can work (unlike arguing with people what kind of information the must provide or what it is that they want to do and them not explaining what they did so far), and it is by far less effort to sort out a well posed question that is easy to understand, or at least make a decent attempt.
Dear akohlmey, I have a further question about this. How will I adjust the relevant calculation parameters of zbl？
I set zbl (1,2)=1 for the following potential function, but I think the resulting force and potential are still smaller than the values of the zbl potential function itself. So I want to change the relevant parameters of the zbl in meam.
FeN.zbl.meam (773 Bytes)
library.meam (423 Bytes)
Fe-N_calculated_by_zbl.txt (26.4 KB)
results_meam.txt (24.9 KB)
What parameters are you talking about? And how would you want modify them?
ZBL has no parameters outside the number of electrons which is based on the choice of element.
Dear akohlmey, I want to change inner and outer in pair_style zbl command — LAMMPS documentation.
In previous attempts, changing these two values would result in different results for force and potential energy.
Sorry, I raised a very foolish question. I just realized that these two parameters only affect the magnitude of potential energy, not their impact. Thanks for your answer very much.
FWIW, pair style zbl has nothing to do with the ZBL functionality in pair style mean. Both are completely independent code.
Thanks very much.