Added Axilrod-Teller manybody potential

[sergeylishchuk]

Purpose

Adds Axilrod-Teller-Muto manybody potential (pair_style atm).

Author(s)

Sergey Lishchuk.

Backward Compatibility

Does not break backward compatibility for inputs.

Implementation Notes

The code is adapted from in-house software.

Post Submission Checklist

Please check the fields below as they are completed

• The feature or features in this pull request is complete
• Suitable new documentation files and/or updates to the existing docs are included
• One or more example input decks are included
• The source code follows the LAMMPS formatting guidelines

Further Information, Files, and Links

Put any additional information here, attach relevant text or image files, and URLs to external sites (e.g. DOIs or webpages)

The code was ported to LAMMPS from the in-house code used in the following publication:

S. V. Lishchuk, "Role of three-body interactions in formation of bulk viscosity in liquid argon", J. Chem. Phys. 136 (2012) 164501.

DOI: 10.1063/1.4704930

[sjplimp]

Hi Sergey - I made several changes/comments in the src and doc (txt) file
where I was unclear on what you were explaining or trying to encode.

[sjplimp]

Hi Sergey - I looked at the latest changes for the Axilrod-T-M
I think I understand
the rules on nu and the 6-way symmetry better now.
The diagram helps.

So I pushed changes to 2 files that I think are correct,
but you should check them carefully:

pair_atm.cpp
pair_atm.txt

The doc file is to explain how to set K with
the wildcards, like I,J allow it, as explained
in pair_coeff.txt. It also explains the
cutoff criterion more clearly. Which is
different that what the code was doing
(I changed the code as well). I don't think
it is possible to use the criterion that
r12r23r31 < cut^6. A particular I,J pair
will not be in the neighbor list if it
if rIJ > cut, regardless of how small
rIJ or rKK are.

I also changed some things in the code
to insure the nu = 0.0 initialization
is done correctly and the restart info
is stored correctly. I think you had
it happening for every pair_coeff command
which would wipe out previous settings.

The answers in the log.9Jul18 log files
changed with the new version. I think
that is probably due to the changed
cutoff criterion. But again you
should check everything carefully.

Also, for computational efficiency
I wonder if it is possible to use
a half neighbor list, and not re-compute
the AMT term 3 times for every triplet
(once per I,J,K). Is it possible
to just compute each I,J,K triplet
once (maybe with a value of nu that
is 3x larger) and get the same force
on every atom? Or are the forces not
identical for all 3 permutations of
I,J,K due to the different angles?
I.e. only the sum of the 3 terms is
"symmetric" ?

Another thing you should check is
if you get the same answers when
the newton flag on on versus off.
And you should check that restarting
works as it should.

Thanks,
Steve

[lammps-jenkins]

Can't set status; build succeeded.

[akohlmey]

@sjplimp there seems to be a typo in the pair_atm.txt file (looks like it has ":per" instead of ":pre"):

Converting pair_atm.txt ...
ERROR: Unrecognized command per
Makefile:110: recipe for target 'pdf' failed
make: *** [pdf] Error 1

[sergeylishchuk]

Hi Steve,

I made some changes to source and documentation:-

1. The right way to calculate cutoff in ATM interactions is to limit the product r12 r23 r31 which enters the denominator of the potential (see e.g. E.Rittger, Comp. Phys. Comm. 67 (1992) 412). Rigorously, this is not possible in finite-size simulations, so the additional cutoff on maximum distance between any pair of atoms is introduced, typically half box length if there are no neighbour lists or the cutoff of the accompanying two-body potential, which is typically larger than three-body cutoff. Either way, the additional cutoff value does not affect the results much because it is very rare (Boltzmann factor!) event that two particles come close enough so that the cutoff length for the third particle becomes very large, larger than eg. 2-body cutoff or distance between neighbours in the neighbour list. On the other hand, use of proper three-body cutoff significantly reduces computation due to skipping computations for which each distance (r12 r23 r31) is close to 2-body cutoff, so that their product is large and the potential (~r^{-9}) is very small. Your code if (r6 > cutoffsq) continue; makes no sense to me because the quantities of different dimensions are compared (L^6 and L^2). So I returned cutoff criterion to include only particles within the neighbour list for which r12 times r23 times r31 < cut^6 is satisfied, the doc file now says: "The potential is calculated if atoms J and K are in the neighbor list of atom I and the distances between atoms satisfy rIJ rJK rKI > cutoff^3.".

2. Each triplet is now calculated only once. This was my original intention, but I overlooked that other permutations were included.

3. I enforced the requirement for newton flag to be on, as in other 3-body interactions.

4. In the doc I added the comment about the need to repeat pair_coeff command after restart when using ATM with hybrid/overlay. Once it is done, it works properly.

Thanks,
Sergey.

[sjplimp]

Hi Sergey @sergeylishchuk - thanks for the changes, we're close I think. However I
still think there is a conceptual problem with this:

The potential is calculated if atoms J and K are in the neighbor
list of atom I and the distances between atoms satisfy rIJ rJK rKI >
cutoff^3.

First, I think that you want to insure the same set of IJK triplets
are computed (on a given timestep), independent of the neighbor list
skin size. If the only criterion is rIJ rJK rKI < cutoff^3, then with
a large skin value, you could compute triplets that would not show up
in a neigh list with a small skin.

So I think you need to compare rIJ and rIK to the force cutoff used to
create the neighbor list (cutoff_global?) and not consider the triplet
if either exceeds it.

Second, I assume you want to produce answers that are rotationally
invariant. I.e. if you rotated the simulation box so that x->y, y->z,
z->x, you want the same set of IJK triplets to be computed (on a given
timestep) as without the rotation. But that rotation will change
which atom is "I" in the outermost loop of Pair::compute(). I.e. the
one to the "left" of the other 2, due to your geometric test. Which
is fine, except I think that means the cutoff criteria actually needs
to also be symmetric, i.e. each of the 3 distances < cutoff. Which
means include rJK < cutoff with the first 2 tests above. If you only
check 2 of the 3 distances, then I think an IJK triplet could be
disallowed for some rotations and not for others?

Note that if you end up having to check all 3 distances < cutoff, than
testing rIJ rJK rKI < cutoff^3 is no longer needed.

Steve

[akohlmey]

@sergeylishchuk @sjplimp i've updated the branch to the latest master and thus resolved the indicated merge conflicts. in addition, i've integrated to documentation properly into the (refactored) manual and made a few formal corrections. there are still some SJP comments in the code, so i assume, those still need to be worked on. if not, please remove. thanks.

[sjplimp]

will update with latest version of Sergey and my files,
should address all the remaining Qs

[sjplimp]

@sergeylishchuk I just pushed what I think is our final changes, using the latest files you sent me. You should check to see if you agree. @akohlmey you can check if you think it is ready to merge

[sergeylishchuk]

@sjplimp @akohlmey I agree with #44bda24, no more Qs.

[akohlmey]

i have one more question: there are significant differences in energy and pressure when running the example input with 1 MPI rank vs. using 4 MPI ranks or when i add a small translation via the offset keyword to the lattice command. is this an intrinsic property/problem of the methodology employed or is there perhaps something missing in the example input?

[sjplimp]

good catch @akohlmey - the offset case is likely different b/c the velocities are initialized based on atom position. However the 1 vs 4 proc difference, even on step 0, for the as-is in.atm script is troubling.
@sergeylishchuk any ideas? You might want to check that the same # of 3-body interactions are computed on step 0, independent of proc count. I.e. put a counter in the compute() loop in pair atm, and then sum the count across procs. I'm thinking that should be identical for all proc counts, every step, until round-off issues set in.

[sergeylishchuk]

@sjplimp File "pair_atm.cpp", line 121:
if (rij2 > cutoff_squared) continue;
The code reaches this line the same number of times at different proc counts, but passes the "if" test different number of times.

The code works fine for very small pair cutoffs, e.g.
pair_style hybrid/overlay lj/cut 1.5 atm 1.5 2.5
but the problem appears if any of them is large, e.g. either
pair_style hybrid/overlay lj/cut 4.5 atm 1.5 2.5
or
pair_style hybrid/overlay lj/cut 1.5 atm 4.5 2.5

Why this might be?

[sjplimp]

@sergeylishchuk Take a look at the pair_atm.cpp I just pushed. The issue with your 2 counters is that the 1st inner loop was to jmax-1, not jmax. This logic was correct (b/c the 2nd inner loop did the right thing), but when running in parallel, the order of what atoms appear in each atom's list of neighs will be permuted, depending on how many procs you run on. Thus you could get a different count for my count1 vs count2 (your previous message). Changing the first loop to jmax gives identical counts
for count1, count2 on all proc counts.

However the real issue is with count3 = # of triplets computed,
which was still different. The change to the 2nd loop, to loop over all neighs, skipping jj==kk, seems to fix the problem. Now energy/pressure are identical for any input cutoffs on any proc count, for up to 50 steps. So the logic of trying to do the inner 2 loops with half the computation is at fault I think, Probably has to be 2x more (brute force double loop). Probably b/c the order of neighs for a single atom is not independent of proc count. I think the order you see the 2 atoms (J,K) matters when doing your tests for which has the smaller x value. If you see them as (K,J) (and not J,K b/c you only loop over half of them), then the code doesn't decide to compute the same triplets?

[sergeylishchuk]

@sjplimp Thank you for changes and explanation, now it looks working. Originally I used the loop structure from "pair_sw".

[sjplimp]

ok @sergeylishchuk - here I think is a clearer explanation of the problem.
Assume atoms IJK have x coords xI < xJ < xK. You want that triplet to be
computed exactly once. The code is setup to compute it on the proc
that owns I. But the short list of neighs of I can be in any order.
Doing the double loop over neighs the "half" way we were (i2 = 1 to
iNeigh), (i3 = i2+1 to iNeigh), will not find the triplet if K comes
before J in the list b/c then xK > xJ. Doing the double list the full way
it is now (i2 = 1 to N), (i3 = 1 to N, skip i2 = i3), will find it
exactly once.

However it should also be possible to do the half double loop and find
it exactly once if we check that xI < xJ and xI < xK, but do no check
for xJ < xK. We would still have to treat ties for xI = xJ or xI = xK
carefully. Do you see any problems with that logic?

[sergeylishchuk]

@sjplimp How about leaving the previous loop structure but replacing rjk test by rik test (attached)? I hope this implements what you suggest. Energy and pressure at different nprocs look fine for me.

pair_atm.cpp.gz

[sjplimp]

yes, that is a match to the end of my last comment - I cleaned up this version (no counting)
and pushed it, along with new log files - now we should be ready to go @akohlmey

[stanmoore1]

@akohlmey this style never got added to the .gitignore file, can we do this before the stable release?

[akohlmey]

On Fri, Dec 7, 2018 at 6:28 PM Stan Moore ***@***.***> wrote: @akohlmey <https://github.com/akohlmey> this style never got added to the .gitignore file, can we do this before the stable release?

yup. just been discussing something similar with @sjplimp. as for changes to the code, there would have to be a serious bug to warrant making changes there. ...and then the clock would be reset and we need to wait at least another week. but it looks like we already have some other "cosmetic" changes pending, i.e. no changes to the code itself. so we may have one more set of only cosmetic merges next week, step the version again, and declare that version immediately as stable. i think this makes sense after all the other cleaning up and making things consistent. that also gives me an opportunity to update a few more log files (not that that is crucial, but due to some corrections to the nose-hoover code, there are some small differences, especially for fix npt). axel.

…

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#987 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARp46H_t_OoRg9ILI0wuc-J7RzWUmfGks5u2vmggaJpZM4VCU72> .

-- Dr. Axel Kohlmeyer akohlmey@gmail.com http://goo.gl/1wk0 College of Science & Technology, Temple University, Philadelphia PA, USA International Centre for Theoretical Physics, Trieste. Italy.