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

Implement HMA compute in LAMMPS #1503

Merged
merged 41 commits into from Aug 21, 2019

Conversation

@ajschult
Copy link
Collaborator

commented Jun 10, 2019

Summary

Harmonically Mapped Averaging (HMA) allows crystal properties (energy, pressure, heat capacity) to be computed more efficiently (more precision in the same time or the same precision in less time). Other benefits include reduced finite-size effects, reduced truncation effects, reduce timestep inaccuracy, shorter decorrelation time and faster equilibration.

Related Issues

n/a

Author(s)

Andrew Schultz, University at Buffalo (ajs42@buffalo.edu)
Apoorva Purohit, University at Buffalo (apoorvap@buffalo.edu)
David Kofke, University at Buffalo (kofke@buffalo.edu)

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

n/a

Implementation Notes

HMA (for crystals) is implemented as a new compute (hma). The compute needs to be created when the atoms are at their lattice sites, but data need not be collected until desired. HMA pressure requires input of a harmonic pressure. Heat capacity requires the second derivative (with atom position); this capability is added for the Pair class and implemented in the PairLJSmoothLinear class. Actually obtaining the heat capacity requires combining two of the terms produced by the compute.

Results of HMA in LAMMPS were verified against "conventional" LAMMPS results (averaging the potential energy or pressure directly) and also against MD and MC with our own Etomica simulation package. Several of the conventional vs. HMA comparisons are part of the online presentation linked below.

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

https://link.aps.org/doi/10.1103/PhysRevE.92.043303
http://rheneas.eng.buffalo.edu/~andrew/lammps-hma/

@ajschult ajschult requested a review from sjplimp as a code owner Jun 10, 2019
@akohlmey akohlmey added the needs_work label Jun 10, 2019
@akohlmey

This comment has been minimized.

Copy link
Member

commented Jun 10, 2019

@ajschult thanks for your pull request submission. please note, that your changes seem to be based on outdated files and are thus undoing several corrections in the current master. please review and correct.

also, i am concerned about adding the Pair::single2() function and then the fact that it is added for just one rather infrequently pair style. that makes your contribution at this stage not very usable.
@sjplimp what do you think about this?

@ajschult

This comment has been minimized.

Copy link
Collaborator Author

commented Jun 10, 2019

@ajschult thanks for your pull request submission. please note, that your changes seem to be based on outdated files and are thus undoing several corrections in the current master. please review and correct.

OK, it looks like there was a bad merge at some point. I've unbroken the files with issues.

also, i am concerned about adding the Pair::single2() function and then the fact that it is added for just one rather infrequently pair style. that makes your contribution at this stage not very usable.

Right, our plan is to implement single2 for more pair styles. Also, only the heat capacity component of HMA depends on single2. The energy and pressure work fine with any pair style.

@ajschult ajschult requested review from junghans and rbberger as code owners Jun 11, 2019
@sjplimp

This comment has been minimized.

Copy link
Contributor

commented Jun 11, 2019

Right, our plan is to implement single2 for more pair styles.

Per Axel's Q, is there any way to do this (accurately) w/out a new method per pair style?
E.g. finiite difference for the Hessian term?

Copy link
Member

left a comment

Please find some nitpicks and questions associated with code to make things more consistent with the rest of LAMMPS for both, users and developers.

src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/USER-HMA/compute_hma.cpp Outdated Show resolved Hide resolved
src/pair.cpp Outdated Show resolved Hide resolved
src/pair_lj_smooth_linear.cpp Outdated Show resolved Hide resolved
@akohlmey

This comment has been minimized.

Copy link
Member

commented Jun 11, 2019

@ajschult to sort out the failing integration check for the docs, you need to add entries to doc/utils/sphinx-config/false_positives.txt. You should check the output of make -C doc spelling or look at: https://ci2.lammps.org/job/lammps/job/pull-requests/job/build-docs-pr/432/warnings3Result/

@ajschult

This comment has been minimized.

Copy link
Collaborator Author

commented Jun 12, 2019

Per Axel's Q, is there any way to do this (accurately) w/out a new method per pair style?
E.g. finiite difference for the Hessian term?

In theory, yes. In practice, it would be expensive. The analytic approach costs O(N * n) (N atoms, n neighbors, similar to an MD step), while a finite difference approach would cost O(N * N) (move each atom, measure force on each atom; this appears to be the approach taken in USER-PHONON/dynamical_matrix.cpp). Given that the point of HMA is to increase the efficiency (better precision from the same configurations), adopting an expensive approach is likely to wipe out HMA's advantage at many conditions (better precision with more expense; a user could get the same by simply running more steps). This would get worse with larger systems.

In terms of accuracy, we haven't attempted a finite difference approach, and so I'm not sure about if accuracy would be an issue. If accuracy was an issue, there are ways to improve the accuracy of finite differences (using multiple steps), but those would also increase the expense.

@Adrian-Diaz

This comment has been minimized.

Copy link
Collaborator

commented Jun 12, 2019

It could be very circuitous to create a strategy that does not simply define a virtual base function given that you seem to be implementing these directly into pair styles; which seems straightforward I suppose and may be the best way. However the function single2 can absolutely use a more formal name given its purpose.

@ajschult

This comment has been minimized.

Copy link
Collaborator Author

commented Jun 25, 2019

Is it necessary for the single hessian function to compute energy and force as well? why is this not relegated to the single function in your code to keep the functionality as specific as possible?

It could be; the new single hessian function needs the force but could call the existing single function to get it. But this would be less efficient because the distance would be computed in each function (and taking sqrt is more expensive than most of other operations in either function).

@Adrian-Diaz

This comment has been minimized.

Copy link
Collaborator

commented Jun 25, 2019

@ajschult why not pass the requisite computations for the hessian as an input to the function?

@ajschult

This comment has been minimized.

Copy link
Collaborator Author

commented Jun 25, 2019

@ajschult why not pass the requisite computations for the hessian as an input to the function?

@Adrian-Diaz I don't see how that helps. The things it needs (that single computes) are fforce, r (and rinv, r2inv, r6inv, which might be efficiently computed from r); single provides fforce, but not r. Calculation of r would have to be duplicated unless single were modified to return r.

@sjplimp

This comment has been minimized.

Copy link
Contributor

commented Jul 1, 2019

I think this is nearly ready-to-go. Please do change the function name to single_hessian().
Re: Can you provide a list of potentials that would be of most interest to LAMMPS users,
I suggest you put this on the command doc page, like fix adapt does. I.e. make
a small table (one-line at this point) of which pair styles have been done. And say
that it is easy to add new ones and request users contact you when they have one
they want. Also, you didn't do lj/cut (the most common pair style). Is that b/c it has
to be a pair style that goes smoothly to 0.0 for energy/force at the cutoff? Or would
lj/cut with pair_modify shift/yes work? In any event, explain that in the doc page as well.

ajschult added 3 commits Jul 2, 2019
Add a table of pair styles that implement single_hessian
@ajschult

This comment has been minimized.

Copy link
Collaborator Author

commented Jul 3, 2019

I think this is nearly ready-to-go. Please do change the function name to single_hessian().

done

I suggest you put this on the command doc page, like fix adapt does. I.e. make
a small table (one-line at this point) of which pair styles have been done. And say
that it is easy to add new ones and request users contact you when they have one
they want. Also, you didn't do lj/cut (the most common pair style). Is that b/c it has
to be a pair style that goes smoothly to 0.0 for energy/force at the cutoff? Or would
lj/cut with pair_modify shift/yes work? In any event, explain that in the doc page as well.

Right, lj/cut is not a good candidate due to truncation effects, even with pair_modify shift/yes (the second derivative would still be infinity at rc and so the resulting heat capacity would be inaccurate and probably inefficient. I've added a bit in the doc about the effect of truncation on efficiency as well as the list of pair styles that have single_hessian and a request to contact us if there are additional pair styles that people have particular interest in.

@athomps

This comment has been minimized.

Copy link
Contributor

commented Jul 9, 2019

@ajschult This is a great addition to LAMMPS. I have a couple of suggestions:

  1. The name Pair::pairTensor() is not very descriptive, and the Pair/pair bit is confusing. For any twobody energy function u(R) that only depends on radial distance R, it combines du/dR, d2u/dR2, and the vector (Rx,Ry,Rz) to form a kind of Hessian matrix d2u, hence I suggest hessian_twobody(). Later, you may need to add hession_threebody(), and so on.

  2. Regarding potentials, I suggest doing it for lj/cut and all the other LJ variants first, since it is essentially just cutting and pasting from PairLJSmoothLinear::single_hessian(). After that, I suggest some of the most popular potentials for real materials: EAM, Stillinger-Weber, Tersoff. Each of these will probably require you to add some extra functions to the Pair class, akin to Pair::hessian_twobody().

  3. The fact that the Hessian is discontinuous at Rcut for some potentials does not seem like a fundamental barrier to adding them. It might be good to indicate in the documentation which ones are well-behaved everywhere and which ones have discontinuities.

  4. This compute does analytically what the dynamic_matrix commands does using finite differences. Perhaps it makes sense to add that connection to both doc pages.

@stanmoore1

This comment has been minimized.

Copy link
Contributor

commented Aug 19, 2019

@akohlmey what is the status here? It looks like you requested changes, have they been implemented?

@akohlmey

This comment has been minimized.

Copy link
Member

commented Aug 19, 2019

@sjplimp

This comment has been minimized.

Copy link
Contributor

commented Aug 20, 2019

@ajschult All looks good - thanks for responding to the various suggestions.

@sjplimp sjplimp assigned akohlmey and unassigned sjplimp Aug 20, 2019
@akohlmey akohlmey added enhancement and removed needs_work labels Aug 20, 2019
akohlmey added 4 commits Aug 20, 2019
@akohlmey akohlmey merged commit cfa9179 into lammps:master Aug 21, 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
Projects
None yet
7 participants
You can’t perform that action at this time.