Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revised implementation of a new atomic stress definition for correct computation of heat flux with many-body interactions #1704

Merged
merged 35 commits into from Nov 19, 2019

Conversation

@donatas-surblys
Copy link
Collaborator

donatas-surblys commented Oct 8, 2019

Summary

This PR is a revised version of PR #1462 according to advice from @athomps. I created a new branch, as little code was shared with the previous PR.
Details and discussion can be found in PR #1462, but in brief, this implements a new "centroid atomic stress" definition, which enables compute heat/flux to produce correct values for many-body interactions, such as angles, dihedrals and impropers, while keeping the original atomic stress intact. This is based on the work by our group [1], and is also equivalent to a recent work and implementation by Boone et al. [2,3], although they did not refer to it as a new atomic stress definition.

Related Issues

Supersedes PR #1462

Author(s)

Donatas Surblys (Tohoku University) surblys.donatas at gmail dot com

Licensing

By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1).

Backward Compatibility

Full backward compatibility.

Implementation Notes

vatom array is left untouched, and a new per-atom 3x3 "centroid" virial array, cvatom, is created in angle, dihedral and improper classes to store values computed by the new definition. The pressatomflag in compute classes can now have a new value "2", to indicate that the compute calculates per-atom centroid virial. The vflag in integrate.cpp or min.cpp can also have value "8" (and other combinations setting the appropriate bit), to indicate that per-atom centroid virial needs to be computed. Therefore, cvatom is only allocated and computed when needed.

Currently, only angles, dihedrals and impropers are implemented. This should work well enough, unless pair styles with many body interactions, such as AIREBO are used. The USER-OMP versions are also implemented. Because I needed to add an argument to the ThrOMP::ev_setup_thr, I had to update all of the pair/bond/angle/dihedral/improper files in USER-OMP, so the number of changed files is quite large.

A new compute centroid/stress/atom, that has 9 components per atom, has been constructed and is mostly identical to compute stress/atom. Compute heat/flux has been modified to accept the new centroid/stress/atom compute.

Before continuing to implement centroid virial for tersoff/AIREBO and similar, I would like to get some feedback from the developers if this approach is acceptable.

Post Submission Checklist

Please check the fields below as they are completed after the pull request has been submitted. Delete lines that don't apply

  • The feature or features in this pull request is complete
  • Licensing information is complete
  • Corresponding author information is complete
  • The source code follows the LAMMPS formatting guidelines
  • Suitable new documentation files and/or updates to the existing docs are included
  • The added/updated documentation is integrated and tested with the documentation build system
  • The feature has been verified to work with the conventional build system
  • The feature has been verified to work with the CMake based build system
  • A package specific README file has been included or updated
  • One or more example input decks are included

Further Information, Files, and Links

[1] https://doi.org/10.1103/PhysRevE.99.051301
[2] https://doi.org/10.1021/acs.jctc.9b00252
[3] https://github.com/WilmerLab/lammps/tree/corrected_heatflux/src

@donatas-surblys donatas-surblys requested review from akohlmey and sjplimp as code owners Oct 8, 2019
@athomps

This comment has been minimized.

Copy link
Contributor

athomps commented Oct 8, 2019

I like this a lot better. However, I am not convinced that it is necessary to add cvatom as a separate variable. Why not reuse vatom, since it is also a 3x3 array? You would just need to have switches within angle.cpp etc., to control whether the centroid or regular per-atom stress is computed.

@akohlmey

This comment has been minimized.

Copy link
Member

akohlmey commented Oct 8, 2019

I like this a lot better. However, I am not convinced that it is necessary to add cvatom as a separate variable. Why not reuse vatom, since it is also a 3x3 array? You would just need to have switches within angle.cpp etc., to control whether the centroid or regular per-atom stress is computed.

correction. vatom is dimensioned nmax by 6.

i am in favor of having the storage separately. reusing existing storage locations for two different purposes is just asking for trouble, since we don't have proper encapsulation in LAMMPS and keeping track of what is what would become a maintenance nightmare. conserving every little byte of storage is not as much of an issue than it used to be.

@athomps

This comment has been minimized.

Copy link
Contributor

athomps commented Oct 8, 2019

I stand corrected. I misread the code. I also agree that reusing the vatom array is not necessary. I thought it might simplify the code, but it does not. Therefore, I approve this request.

@donatas-surblys

This comment has been minimized.

Copy link
Collaborator Author

donatas-surblys commented Oct 9, 2019

Great! I'll finish implementing this for pair styles and will update the documentation. Might need a bit of time to run some simulations to make sure everything is OK.

@akohlmey

This comment has been minimized.

Copy link
Member

akohlmey commented Oct 9, 2019

Great! I'll finish implementing this for pair styles and will update the documentation. Might need a bit of time to run some simulations to make sure everything is OK.

@donatas-surblys i just enabled additional automated testing for this PR. it will run one set of tests, where the check is if the executable compiles with certain settings and runs over a set of inputs completes with different parallelization settings, and a second, more thorough set of tests, where also the results are compared to a reference. it roughly covers the code paths in LAMMPS for which there are example inputs and a few other basic cases. that will make sure that at least for common operations, everything will remain working.

@donatas-surblys donatas-surblys force-pushed the donatas-surblys:many-body-atomic-stress-rev branch from 7dfd93d to 2d75e6b Nov 12, 2019
@donatas-surblys donatas-surblys requested a review from rbberger as a code owner Nov 12, 2019
@donatas-surblys

This comment has been minimized.

Copy link
Collaborator Author

donatas-surblys commented Nov 12, 2019

I’ve rebased the PR to the newest master.

Missing support for angle/dihedral/improper hybrid was added.

Initial support for centroid atomic stress for pair styles has been added. A new flag “cntratmstressflag” has been added to the Pair class. Value of 1 indicates that the centroid atomic stress is identical to normal atomic stress, and the values are stored in vatom. Value of 2 indicates that centoid atomic stress has been implemented for the pair style, and that the values are stored in cvatom. Value of 4 indicates that centroid atomic stress has not been implemented for that pair style. The default has been set to 4. I initially wanted a value of 0 to indicate a lack of implementation, but it made it difficult to interact gracefully with pair hybrid. The “cntratmstressflag” might also have a value of 3, indicating to set both vflag_atom and cvflag_atom to true, as sometimes the tally functions are not called otherwise.

Initial support has also been added to USER-OMP pair styles, but unfortunately I did not find a way to access “cntratmstressflag” from inside ThrOMP::ev_setup_thr, therefore it is checked if cvatom is a null pointer instead. This has an unfortunate side effect, that USER-OMP pair styles can no-longer tell the difference between cntratmstressflag = 2 and cntratmstressflag = 3, so the current plan is to always at least allocate vatom array when centroid atomic stress is needed.

I’ve set cntratmstressflag = 1 to non-kspace pairwise-only pair styles in OPT, USER-OMP, USER-FEP and CLASS2, though I’m sure I’ve missed a few. Kspace is kind of tricky, because while it should be OK for total system heat flux, it becomes theoretically shaky when applied to a local control volume. I have some vary preliminary results that looked OK for local heat flux with kspace, but there’s no rigid theoretical or empirical work to back it up. I guess we might as well enable them, as this would be in line with compute stress/atom for pairwise interactions.

I’ve added a check to the compute centroid/stress/atom, to error out if a unsupported pair style is used when the pair virial is requested. The flag logic is set up so that if the error line was changed to a warning or deleted, then unimplemented styles would simply fallback to the old behaviour.

I’ve tried to update the documentation for compute stress/atom and compute heat/flux to include the information on the new compute centroid/stress/atom. I’ve converted all the equation images to equations and added inline math where possible. I’ve tried to keep as much of the original text and structure as possible, but had to rearrange some things. I am considering adding a warning that when computing local heat flux (or pressure), only control volumes in bulk regions are likely to give reliable results. Let me know if more clarity or details are needed.

I am also not sure about the citation policy. I’ve added our group’s paper [1] that gives details on the implementation to the documentation of compute stress/atom, and also added it to the documentation of compute heat/flux together with the paper by Boone et al. [2] inside a warning about many-body interactions, as it is closely related. I haven’t added the work by Fan et al. [3], as the pair styles they deal with are no yet implemented.

In fact, we should set things up so that both possible stress definitions can be accessed and thus test of re-doing old simulations remains possible.

@akohlmey The current PR implementation allows for this. You can use compute stress/atom for the original atomic stress and compute centroid/stress/atom for the new definition simultaneously. compute heat/flux will do the right thing when passed either one of them.

[1] https://doi.org/10.1103/PhysRevE.99.051301
[2] https://doi.org/10.1021/acs.jctc.9b00252
[3] https://doi.org/10.1103/PhysRevB.92.094301

@akohlmey

This comment has been minimized.

Copy link
Member

akohlmey commented Nov 12, 2019

How about changing the somewhat awkward name of the flag in pair.h from cntratmstressflag to the more palatablecentroidstressflag?

Copy link
Contributor

athomps left a comment

@donatas-surblys very nice work. I particularly like how you added the centroid stress flag to Pair class.

  1. I agree with @akohlmey about centroidstressflag.
  2. On doc page, change forth to fourth.
  3. In pair.h and maybe elsewhere, change pairwise to twobody, for clarity.
@sjplimp

This comment has been minimized.

Copy link
Contributor

sjplimp commented Nov 12, 2019

@donatas-surblys Thanks for working thru the details on this, with all the feedback from the developers. It looks like a very nice addition to LAMMPS. 2 Qs/comments on the doc pages.
It looks like all the doc page changes are in 2 places, compute heat/flux and compute stress/atom.
(a) On the compute stress/atom page I think referring to the new styles as "many-body" terms may be confusing to LAMMPS users. Since in your usage it refers to bond/angle/dihedral, but not many-body pair styles (yet). The LAMMPS terminology is pair styles (which can encode pairwise or manybody interactions) versus bonded interactions (bond,angle, dihedral, etc with fixed topology). So if you use those terms in the doc page text, then you can explain for bonded interactions that the many-body nature of the stress computation is what is being altered/corrected in the new compute. Per your Restrictions note, you could also explain your formulation could also be applied to many-body pair styles, but that hasn't been done yet for any non-pairwise pair style (is that correct? does that also include EAM, which is technically many-body but pairwise-like?).
(b) really a question for @rbberger - does the new compute centroid/stress/atom still need to be added to other doc pages manually at this point, so it will appear in the table listing of all computes

@athomps

This comment has been minimized.

Copy link
Contributor

athomps commented Nov 12, 2019

What Steve said | sed s/pairwise/two-body/g :-)

@donatas-surblys

This comment has been minimized.

Copy link
Collaborator Author

donatas-surblys commented Nov 13, 2019

Thank you for the feedback.

I’ve changed the flag name.

Changed pairwise to two-body where appropriate in source code comments.

I’ve slightly modified both the compute stress/atom and compute heat/flux documentation. Basically got rid of using “many-body” to describe both bonded and non-bonded interactions, and mostly switched to simply using angle/dihedral/improper, as that is what is currently implemented. I’ve also update the restrictions section according to the advice of @sjplimp. Let me know if this is what you had in mind.

is that correct? does that also include EAM, which is technically many-body but pairwise-like?

@sjplimp The EAM potential is an interesting example. It all comes down on how you distribute per-atom potential. If you assume that only atoms i and j are attributed potential energy, even though it also depends on density, then you should able to treat it as a pairwise-like potential and it will probably work out. That said, I’d like to do some testing and a more rigid theoretical checks before enabling such cases.

It seems that the “lammps/pull-requests/cmake/cmake-kokkos-cuda-pr” check failed, but looking at the logs it seems to be a timeout error? Am I reading this right?

@akohlmey

This comment has been minimized.

Copy link
Member

akohlmey commented Nov 19, 2019

is that correct? does that also include EAM, which is technically many-body but pairwise-like?

@sjplimp The EAM potential is an interesting example. It all comes down on how you distribute per-atom potential. If you assume that only atoms i and j are attributed potential energy, even though it also depends on density, then you should able to treat it as a pairwise-like potential and it will probably work out. That said, I’d like to do some testing and a more rigid theoretical checks before enabling such cases.

Please note, that we just merged some changes to the EAM styles, that allow partitioning of the energy between pairs. The tricky part is the embedding energy term, which we now evenly distribute over all neighbors within the force cutoff as shown in PairEAM::single() (look for the use of numforce).
You can see the total changeset with:
git diff fd9da6f934b59ec9c0bd06371472caed1066128e..48894884124e319a89003741d967af422484c625 src/MANYBODY/pair_eam.{cpp,h}

It seems that the “lammps/pull-requests/cmake/cmake-kokkos-cuda-pr” check failed, but looking at the logs it seems to be a timeout error? Am I reading this right?

I will have closer look at the changes in the PR now and update it to the current master to make certain, there are no conflicts with the recent changes to the documentation.

akohlmey added 2 commits Nov 19, 2019
@akohlmey

This comment has been minimized.

Copy link
Member

akohlmey commented Nov 19, 2019

@sjplimp what do you think? should we include this now (and leave some possible improvements and cleanups for the next patch?) or rather wait some more? I would be in favor of merging it now since it will likely be another month or so until the next patch.

@athomps

This comment has been minimized.

Copy link
Contributor

athomps commented Nov 19, 2019

@sjplimp @akohlmey Yes, I think we should release this. I think once people in the heat transport community get a chance to try it out, @donatas-surblys will get useful feedback that can help guide further updates.

@sjplimp

This comment has been minimized.

Copy link
Contributor

sjplimp commented Nov 19, 2019

@akohlmey @athomps Agree with Aidan, ok to release.

@akohlmey akohlmey merged commit 8f36800 into lammps:master Nov 19, 2019
6 checks passed
6 checks passed
lammps/pull-requests/build-docs-pr head run ended
Details
lammps/pull-requests/cmake/cmake-serial-pr head run ended
Details
lammps/pull-requests/kokkos-omp-pr head run ended
Details
lammps/pull-requests/openmpi-pr head run ended
Details
lammps/pull-requests/serial-pr head run ended
Details
lammps/pull-requests/shlib-pr head run ended
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.