LAMMPS minimization with external force field that depends on particle positions

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am working on a project involving energy minimization for 3D Lennard-Jones binary glass with an external force field. Here below is the description:
(1) The force field is a bending force field, the function of which depends on positions of particles and depends on three variables, a, b and c. At any given configuration during an energy minimization step, I can solve linear equations to update a,b and c, in order to make total force and total torque equal to zero. To summarize, given a configuration, I can use an analytical expression from positions of particles to calculate a,b and c, and then update my force field.
(2) I have tried to use LAMMPS with commands: fix addforce, variable (atom), and compute reduce. I’m struggling with errors and I wonder if I could write a LAMMPS input file using these functions for it or I need to edit the source C++ file.
(3) Do you have any similar examples that I could study from?

Thank you!
Best,
Meng

2

Hi,
I am working on a project involving energy minimization for 3D Lennard-Jones binary glass with an external force field. Here below is the description:
(1) The force field is a bending force field, the function of which depends on positions of particles and depends on three variables, a, b and c. At any given configuration during an energy minimization step, I can solve linear equations to update a,b and c, in order to make total force and total torque equal to zero. To summarize, given a configuration, I can use an analytical expression from positions of particles to calculate a,b and c, and then update my force field.

Have you taken a look at "fix momentum"?
https://lammps.sandia.gov/doc/fix_momentum.htmlfix_momentum.cpp

Even if it doesn't do what you want, it might be a good starting point
for writing a new fix. (Terminology "fix" refers to any code which
controls the way atoms move during a simulation.) The C++ code for
"fix momentum" is located in these two files "fix_momentum.cpp" and
"fix_momentum.h". If you want to write your own custom fix, copy
these files and give them new names

Also check out "fix recenter"
https://lammps.sandia.gov/doc/fix_recenter.html
...and "fix spring"
https://lammps.sandia.gov/doc/fix_spring.html
https://lammps.sandia.gov/doc/fix_spring_self.html

--- how to create your own custom fix (if you need to) ---

replace "FixMomentum" with "FixNAMEofNEWfix" everwhere it appears in
"fix_momentum.h" and "fix_momentum.cpp".
Also change "momentum" to "name_of_new_fix" everywhere it appears.
If your fix requires numerical parameters to be supplied by the user,
then edit the constructor
(IE, FixMomentum::FixMomentum() near lines 33-67 of "fix_momentum.cpp")

The actual code that modifies the movement of the atoms is located in
FixMomentum::end_of_step()
...however if you are simply adding forces to the existing forces
acting on these atoms, then perhaps instead you should put this code
in
post_force() instead of end_of_step(). That's what "fix addforce" does.
The code for fix addforce is "fix_addforce.cpp"

For details, see:

https://lammps.sandia.gov/doc/Developer.pdf

here is another link:
https://lammps.sandia.gov/threads/msg77779.html
(probably irrelevant to what you are trying to do)

There may be other changes you need to make, but hopefully this gives
you some idea what to do.
Basically, find a fix which is similar to what you want to do, copy
the corresponding .h and .cpp files, rename the fix. You don't have
to understand what every line of code does. Just use these files as a
template.

Hope this helps get you started.

Andrew
P.S. Also, if you are unfamiliar with C++, perhaps these web sites
will help. (I did not choose these web sites very carefully)
http://www.cplusplus.com/doc/tutorial/classes/
http://www.learncpp.com/cpp-tutorial/89-class-code-and-header-files/

3

I think you need to write a new fix, either in C++ or Python, to solve
your linear eqs for a.b.c for each atom and apply the forces.

Steve

4

Hi,

Background:
I’m implementing energy minimization for 3D Lennard-Jones binary glass with an external force field. The force field depends on positions of each particle. I have a MD cpp library of our research group and the code can successfully implement Conjugate Gradient on my system, but it’s not fast enough, which makes me want to try LAMMPS. In LAMMPS, I’ve successfully implemented the force field by writing my own fix add force file. Then I use minimize to minimize the total energy of Lennard-Jones plus the force field induced potential energy.

Problem:
I have a problem when using minimize (CG), with error line search alpha is zero, as following:

Minimization stats:

Stopping criterion = linesearch alpha is zero

Energy initial, next-to-last, final =

-5.92534765781 -5.92534768347 -5.92534768347

Force two-norm initial, final = 0.0128405 0.24935

Force max component initial, final = 0.000619363 0.0462415

Final line search alpha, max atom move = 0.000226556 1.04763e-05

Iterations, force evaluations = 18 57

By the way:
1, Also, for different force field amplitude, the minimization sometimes work, reaching ftol = 1e-8;
2, Changing CG to FIRE style, it always reached maxiter even with maxiter = 100000, i.e., the minimum is hard to find in general;

Thoughts and questions:
I’ve read other posts and knew that it was due to that my potential is too complicated. The minimizer cannot further decrease energy in the current Conjugate direction.
(1) However, is there a way to let the algorithm try another new Conjugate direction to explore rather than directly jump out the iterations?
(2) Is it possible to implement my own minimize file? Are there any examples to start from?

Best,
Meng

5

Question:
Do you include the energy contributions from the fix in your total energy?
Axel

6

And is the potential continuous in force and energy at
the cutoff? Even a simple LJ potential will not converge
will if it is not shifted to 0 at the cutoff.

Steve

7

Hi,

The LJ potential is a standard shifted-force Lennard Jones potential, which is 1st-order and 2nd-order continuous at cutoff. However, the potential of the force field may not. Here below are the details of the force field:

The force field is a bending force field, and the force on any particle at position (x,y,z) only has z-direction component:
F_z = A[cos(pi*x/L)-2/pi] + ax + by +c,
Where A is a constant amplitude, L is a constant box length, and a, b , and c are close-to-zero numbers that are fixed at current configuration but will vary whenever the configuration changes. At any given configuration during an energy minimization step, I solve linear equations to update a,b and c, in order to make total force and total torque equal to zero. To summarize, given a configuration, I can use an analytical expression from positions of particles to calculate a,b and c, and then update my force field.

I integrate the force function in z direction to calculate its corresponding potential energy. I suppose this force field is “discontinuous” since a, b and c are varying, i.e., my Energy Landscape is varying itself. Thus, I need a Conjugate Gradient algorithm that can constantly update new conjugate directions when the line search along the old direction cannot find a minimum.

Best,
Meng

8
9

Hi, Steve:

Are there any thoughts or suggestions on this?

Thanks!
Best,
Meng

10

You didn’t answer Axel’s Q as to whether you have used fix_modify energy yes to add
the energy for minimization. I assume by your original email you are defining
a atom-style variable for the z force which encodes your equation as an atom-style
variable (as a func of atom coord) that LAMMPS computes on the fly. You must also then use the fix addforce
energy option and encode an equation for the corresponding energy as a 2nd atom-style
variable (as a func of atom coord). If you have done that correctly (and used fix_modify energy yes),
then the minimization should work. If it’s not possible to express your force/energy in that
way (you said something about a fix addforce file - don’t see how you could read in these
forces/energies from a file when they depend on the current atom coords), then

you probably need to write your own fix, similar to fix addforce.

Steve