Hi Brian,

I think I see your point, and yes, it does look the ZBL piece should be accounted for in a few more places for consistency.

Just to make sure we have the same understanding of this, I'll try to summarize what's going on (even though you've already done this nicely). The compute_pair_meam subroutine loops over all possible a-b pairs (where a and b are element types), then over all discrete r values we're using for the spline fit. We call phi_meam(r,a,b), which computes what we can think of as the "uncorrected" pair potential. Then, if we're using second nearest-neighbors, we have to do some more work to "correct" this. Finally, we mix in the ZBL potential.

There are two general ways in which, for a given a-b pair, we may need to compute the a-a or b-b potential as an intermediate step. First, within phi_meam() itself, if the reference structure involves nearest neighbors of like element types, those a-a pairs come directly into the potential. This is true for the C11 and L12 reference structures. For some other alloy structures, like the B1 and B2 structures, all nearest neighbors are of opposite type, so this issue doesn't apply.

Second, the a-a pair is needed for alloys when we're using second nearest-neighbors (NN2). Currently, we don't include the ZBL potential piece here either; instead, we compute the entire potential as if there were no ZBL, and then add in the ZBL potential at the end. This probably matters less, since it only makes a difference when the second neighbors are close enough for the ZBL piece to kick in (so that first neighbors are extremely compacted).

For the first case, here's my original argument for doing things the way they're currently done: the phi_meam(r,a,b) subroutine computes the pair potential that obeys the Rose energy function for the a-b reference structure. The ZBL potential is different form the Rose energy function, and therefore it makes sense to compute everything within phi_meam() (including those a-a terms) as if there is no ZBL potential. Now, however, I think that's wrong (and that you're right). The simplest way to see that this is true is to think about the case where we want to use ZBL for the a-a pair, but not the a-b alloy (so that a-b should obey the Rose energy). If we don't account for the ZBL piece of the a-a potential, the resulting a-b energy will actually not obey the Rose energy because of the perturbation in a-a. Assuming that makes sense, do you agree?

A similar argument holds for the NN2 case, I think.

So, it appears that everything could be made consistent by moving the ZBL correction to the end of the phi_meam() subroutine, instead of the end of compute_pair_meam. Again, do you agree?

If this makes sense, I'll try to make this change after some testing. It should only change results for the C11 and L12 alloy systems, and for NN2 systems under extreme compaction. Thanks for taking such a close look at the code.

Regards,

Greg