Add new fix to compute force and torque due to electric potential

[gsalkuin]

Summary

This PR adds a new fix that computes the force/torque on both charges and dipoles given an electric potential. Currently, the fix efield command does not compute the force on dipoles due to the gradient of the electric field. Instead of requiring the components of the electric field and its gradients, this fix only takes a Lepton expression of the form V(x,y,z).

Related Issue(s)

N/A

Author(s)

Gabriel Alkuino (Syracuse University) - gsalkuino@outlook.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

No known issues.

Implementation Notes

Post Submission Checklist

• 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
• Suitable tests have been added to the unittest tree.
• A package specific README file has been included or updated
• One or more example input decks are included

Further Information, Files, and Links

[akohlmey]

@gsalkuin Thank you for your submission. Also, many thanks for making the effort to adapt to the LAMMPS programming style, providing (and integrating!!) documentation. You've made my job a lot easier. I've made a few more cosmetic changes on top of that and added a couple of tests for the automated testing. So, from the formal point of view this is ready to go it. I'll leave it to other LAMMPS developers that are more competent in this area of simulations to check the implementation itself.

[gsalkuin]

Thank you for the positive feedback and for making the additional changes. This is my first PR and I'm still learning C++, but I hope to make more contributions in the future.

[ndtrung81]

@gsalkuin thanks for the contribution! It's important to see that fix efield with a non-uniform e-field misses a force component. I have only a comment in the source code to make the for loop(s) over nlocal more efficient, and easier to port to accelerator styles.

[gsalkuin]

@ndtrung81 Thanks for the feedback. You are right, I found that the loop is very inefficient when compared to fix efield for the uniform efield case. I think the main cause is lines 195-212 in fix_epot_lepton.cpp. Putting this outside the loop and adding a has_ref flag similar to pair_lepton_coul.cpp made it much faster when not all derivatives are nonzero.

[akohlmey]

@gsalkuin it looks like we have a portability issue with Windows and Microsoft Visual C++
Do you have access to a Windows machine with Microsoft Visual Studio so that you can have a look for yourself?
FWIW, the "Community Edition" of Visual Studio is available cost free to individuals, thus it is straightforward to install the compiler (and IDE) and LAMMPS is well set up to be compiled straight from checking out the git repository, since git and CMake are supported directly.

Otherwise, I can take a look (I have a Windows Virtual machine that I can use for that purpose).

[gsalkuin]

@akohlmey I think the issue is due to a missing #include <array>. I am able to compile it on my Linux and Windows machine now. Hopefully, it passes the checks.

[srtee]

Hi @gsalkuin , this is a great addition to the LAMMPS codebase -- even though I've marked my review as needing addressed, this is really well done (as well as including good documentation and tests!)

My main question is around the intended treatment of periodicity -- should the Lepton expression act on wrapped coordinates or unwrapped? There are problems either way -- acting on wrapped coordinates would impose artificial force-mismatches on molecules straddling periodic boundaries. But acting on unwrapped coordinates would result in a non-periodic simulation, where particles with similar box coordinates might experience completely different forces depending on their image flags. I am not familiar with periodic simulations involving inhomogeneous fields, so I am happy to defer to your expertise -- but periodicity should be more clearly documented.

My other main questions were around having lots of RAIIs being built and destroyed inside the post_force loop, and similar operations -- but if the performance hit is minimal, then you may consider my comments as premature optimisation and reply accordingly :)

Cheers,
Shern

[gsalkuin]

Hi @srtee,

Thanks for the detailed feedback and suggestions! I will try to answer what I can based on my limited expertise. I mostly use LAMMPS for mesh-free simulations in solid mechanics, so I rarely use PBC. The intended use case for this fix is to simulate electric/magnetic control, which usually has the field localized in a region. But that is an excellent point on the force becoming mismatched or unbounded in traditional MD use cases. I used unwrapped coordinates to avoid possible discontinuities in the potential, but perhaps this is not an issue for LEPTON, and wrapped coordinates will be better. I will defer this to your and other MD experts' judgment. I am also not that familiar with the many box-changing commands, but as long as the expression can be written as a function of "x", "y", "z" I think it should be pretty flexible.

As for the post_force() variable-reference checking, I think this is needed due to LeptonUtils::substitute(expr, lmp). If a LAMMPS variable evaluates to zero at a particular time step, LEPTON may remove a reference during optimize() which will cause an error. This is indeed a very costly step. A quick fix I have is to define another function to check if there are no or only constant LAMMPS variables. The performance improvement is pretty significant, but it made the original code more bloated. Maybe there is a way this could be added in lepton_utils.cpp instead.

I will try to address your other concerns when time permits. Thanks!

[akohlmey]

so I rarely use PBC. The intended use case for this fix is to simulate electric/magnetic control, which usually has the field localized in a region. But that is an excellent point on the force becoming mismatched or unbounded in traditional MD use cases. I used unwrapped coordinates to avoid possible discontinuities in the potential, but perhaps this is not an issue for LEPTON, and wrapped coordinates will be better.

If you do not have PBC, wrapped and unwrapped coordinates are the same. Recent versions of LAMMPS will drop image flags != 0 for non-periodic dimensions when reading coordinates from data files.

The code should print a warning if there is a periodic dimension, and there should be an explanation in the documentation that using it for periodic boundaries is likely going to cause problems and unphysical behavior for atoms crossing the boundary.

[stanmoore1]

When the comments/feedback from other reviewers are addressed then I approve this PR. I don't have any additional comments or concerns.

[sjplimp]

Thanks @gsalkuin for this addition. I have 3 cosmetic suggestions.
(1) The first paragraph of the doc page description for both this command and fix efield should mention the other command as an alternative. This is to guide users. E.g. fix efield could say you should use this other command if your efield has a gradient b/c that term is missing. This command could say use fix efield if you have no gradient, or if the gradient effect is small and you want something faster.
(2) Adding an example input script or two in the examples dir, illustrating use of this command, would also be helpful to users.
(3) Naming of the command. Since this command is converting a potential formula to an efield which is then applied to particles (same final step as fix efield), my preference would be to name this command fix efield/lepton. So that when users scan the list of fixes, they appear together alphabetically. If you don't like that, then fix epotential/lepton would be my 2nd choice, since epot is not as clear.

[gsalkuin]

Thanks, @stanmoore1 and @sjplimp!

@sjplimp, I'll work on making the changes you suggested. For (3), I didn't name it efield/lepton to avoid confusion since it doesn't follow the same syntax as efield and efield/tip4p. But if this isn't an issue, I can rename it to efield/lepton.

[sjplimp]

@gsalkuin Having the same or similar syntax is more an issue when the new class is derived from a parent class. In this case it is more a conceptual motivation. Have all the fixes which apply an efield to charged particles, start with the same name. Kind of like all the compute temp variants, many of which have different syntax, but all conceptually calculate a temperature.

Let's see if Stan or Axel have feedback on the (3) name change.

[akohlmey]

Let's see if Stan or Axel have feedback on the (3) name change.

I agree with calling it fix efield/lepton. It helps people look in the right places for what they need. All /lepton styles are by their very nature different from others, so there is not a problem with having a different syntax than fix efield.

[akohlmey]

@gsalkuin I have submitted a pull request in your fork with some updates for this branch.

@srtee can you please confirm that all your concerns have been addressed.
Otherwise, I will merged this tomorrow or as soon as my changes are included.

[srtee]

@akohlmey nothing major from me, just suggestions to documentation and code style (perhaps minor performance improvements). I hope to finish those within 24 hours from this comment but if they're not in before then that's on me, I'm happy to merge.

@gsalkuin very good job on moving the heavy Lepton machinery out of the inner loops!!

Notes to self:

• clarify role of itinerant polarisation in efield (affects energy, shouldn't affect force or else unphysical)
• make those "x, y, z" static constexpr and give better names?

[srtee]

Hi @gsalkuin,

I've put my proposed changes in a PR to your branch, with more detailed comments there. In brief:

• I suggested some code changes that will hopefully make the code (slightly) more maintainable in future
• I moved those if (atom->q_flag) and if (atom->mu_flag) checks back into the nlocal loop -- there's a lot of code duplication when those are moved outside the loop. This was definitely a mistake on my end suggesting moving those checks out, and your original design choice to have those checks inside the loop makes for much more readable code.
• I did suggest a different way of getting and storing the Lepton references -- I feel less confident about this one improving the code maintainability, compared to the other changes I suggested.

I did think about adding things to the documentation but -- it's on users to know the physics and realism of their simulations.

I also checked that the regression test you added does, in fact, fail when I randomly flip one of the force signs in the source code, good job!

To be clear -- if the changes I suggested make the code less sensible to you, rather than more, just let me know in this comment thread and I'm happy to approve the PR as-is. I think we should have this PR merged in pretty soon :)

[gsalkuin]

Hi @srtee,

The changes you made are excellent—they make the code much more concise. I’ve merged your PR. Thanks!