add a complicated force in lammps

Dear all,

I am intending to add an additional force after the standard mechanical force(eam) calculation, then do energy minimization for the system. However, the calculation of this complicated force requires the input of the coordinates of all the atoms in the system and some other input parameters, then compute the interactions, instead of using/searching only the atoms in the neighbor list, No cutoff in this case.

I have been thinking to add the force directly in the pair_eam.cpp code, so the energy minimization can automatically account for this additional force. Another strategy I have been thinking is to add the force to fix_addforce.cpp.

I am not quite familiar with lammps, therefore I hope to stay within its current class, so things won’t get too complicated. So just wondering which should be more promising, or any other way to work this out, or do I need to write a separate .cpp (.h) to add into lammps class.

Appreciate help! Many thanks!

Xue

If you want the force to be applied during the
minimization (or subsequent dynamics), then
it needs to be a fix. You can use fix addforce
if you can write a atom-style variable formula
for the force. See the addforce and variable
doc pages for details. If not (and it sounds
like not from your email), then you would have
to extend fix addforce or write your own fix
to input more details from a file.

Steve

Hi Steve, lammps-users,

Thanks a lot for your help. So yes, I should definitely need to write my own fix codes, to add an additional force before each energy minimization step. However, I am not so experienced in this. Are there any tutorials or example codes that I can follow?

I have written a pseudo-program for what I want to do. This is not a periodical structure, just a cluster composed of thousands of atoms(FCC structure).

This fix should take in two input parameters from the input script.

I don’t know whether I should extend fix_addforce, or write new fix code for the force calculation part. The purpose is to add this force to the eam mechanical force, then do the energy minimization.

For the force calculation: it is as the following:

Define properties based on the two input parameters.

Gather coordinates of every atoms in the system.

Based on the coordinates and the properties defined above, compute the interaction matrix “Mat” for the system. (this is a very large matrix, on the order of 4N*4N)

Perform some linear algebra calculation on “Mat” to compute forces on each atoms in the system.

Add this forces to the conventional mechanical force.

Would it be possible that I can have a simple example code that can do gather atoms coordinates info, assign forces to atoms, add them, then do energy minimization? So I can just fill in my own force algorithm part into the code? Really appreciate your help.

Thanks,

Xue

Are there any tutorials or example codes that I can folw?

See Section_modify.html in the manual.

Based on your pseudo-code below, the hardest part of
what you want to do is solve the linear system. If N
is large, the matrix has to be distributed. Which means
a parallel linear algebra solve, presumably of a sparse
matrix, thus presumably by an iterative method. LAMMPS
has no code to perform such a solve. There are parallel
libraries (Trilinos, PetSC, etc) that do, but we don’t link
to them (they are larger than LAMMPS!). So you will
need to figure out how to interface your fix to some
such library.

Steve

You may also want to consider using Lammps as a library to build your initial solution to the problem.

There may be many potential hurdles for you to face, i.e., proper calculation of the matrix elements, gathering of the new forces for the matrix transformations, etc, thus for troubleshooting purposes you (I would) may find easier to deal with the new methodology outside of Lammps. Only once you are convinced that you have mastered the new method, incorporating it into Lammps will follow naturally.

Carlos

Hi Carlos, Steve, lammps-users,

Thanks for your replies.

The calculation of the new force I wanna add to the eam mechanical force is very time-consuming and memory-demanding. I have written a vanilla style C code for this force, and used GPU computing cuda and cula library, which is working. So I am thinking to incorporate this into lammps code. So lammps can automatically invoke my new force, take care of boundary conditions of the system and do energy minimization. This seems pretty hard, I understand. A lot of details I need to deal with.

I did thought about using lammps as a library. Use lammps to get eam mechanical force, and use my code for the addtional force, then add them together to be the total force on each atom. The process should be like this:

1.tell lammps lib the coordinates of atoms, calculate eam force, go back to my main code.
2.use my code to calculate my force, go back to my main code. (takes long long time)
3.add the two force together, tell lammps lib this new force value, use lammps lib to do energy minimization for one step.
4. go back to 1, repeat until the minimum energy value is reached.

Then my question is first, whether I can use lammps library to do energy minimization with this new force, second, whether this communication be very time-consuming.

I am not sure which approach(directly add force into lammps or use lammps as lib) is better or more efficient. Looking forward to your reply. Many many thanks,

Xue

I could not tell you which approach is more efficient. My expertise embedding code into the Lammps backbone is very limited. Steve, Axel or any other of the developers are the voices to listen to in this case. My advice was just to help you split the potential sources of errors. But if you already have a tested version of your technique that works then better wait to hear what they have to tell you.

Carlos

there are two different issues at hand here:

1) integrating your method to compute forces into an MD engine like LAMMPS
2) making it work efficiently

i think you need to solve point 2) first and not point 1). from the
limited description, it sounds like you have a badly scaling
algorithm. you keep mentioning that you need the information for all
atoms. that sounds like it is not a good approach for a distributed
memory parallel program like LAMMPS. do you have some benchmark data
about how much time it takes for your calculation by itself compared
to the computation of LAMMPS (in serial)? from your description it
sounds like the time spent on the LAMMPS calculation would be small,
if not negligible and should get irrelevant when you scale your
problem size up.
at that point, it won't make much sense to run LAMMPS in parallel and
*then* solving point 1) is easy: compile lammps as a serial library
and/or use the python script interface.

if you *do* want to run in parallel, you have a lot more work ahead of
you, because you first need to look at how to run in parallel. but
personally, i would first try to explore options, to make your
algorithm scale better or modify it, so that it doesn't need all atoms
and turn this into some complex linear algebra issue. your description
sounds like it would scale worse than some simple DFT based quantum
chemistry (which can be approximated and done even faster in different
ways, too), so why would one even *want* to use such a badly scaling
empirical method?

axel.

Thanks for Axel’s reply!

The time spent on my algorithm is quite long(need matrix inversion and multiplications), if I use CPU computing technique. So I modified the code into 1 CPU+1GPU usage, and it works much faster. And the lammps on mechanical force cal time is negligible.

I don’t need lammps to run parallel, 1 CPU is fine. For my added force, 1 CPU is also good, because I will tell the GPU to do all the cal work. I would not wanna deal with MPI stuff… Then in this case, is it easy to integrate my force into LAMMPS to get the energy minimized configuration of the system? I am thinking on modifying fix_addforce.cpp.

Many thanks,

Xue

Thanks for Axel's reply!

The time spent on my algorithm is quite long(need matrix inversion and
multiplications), if I use CPU computing technique. So I modified the code
into 1 CPU+1GPU usage, and it works much faster. And the lammps on
mechanical force cal time is negligible.

I don't need lammps to run parallel, 1 CPU is fine. For my added force, 1
CPU is also good, because I will tell the GPU to do all the cal work. I
would not wanna deal with MPI stuff... Then in this case, is it easy to
integrate my force into LAMMPS to get the energy minimized configuration of
the system? I am thinking on modifying fix_addforce.cpp.

that is not such a good idea in my opinion.
it would probably be much better to use fix_external.cpp via the
library interface:

http://lammps.sandia.gov/doc/fix_external.html

that was designed for such applications. you can then enforce in the
library interface, that only one MPI task is assigned to LAMMPS.

axel.

Hi Axel, and all lammps-users,

Very help advice! I have read through the fix_external, and also the example provided in the COUPLE folder. I saw the code input the DFT force into lammps for MD energy minimization. What I need is not only the external force but also lammps’ own pair style force added together for energy minimization. Would this be possible by using fix_external, instead of fix_addforce ?

Thanks a lot!

Xue

Hi Axel, and all lammps-users,

Very help advice! I have read through the fix_external, and also the example
provided in the COUPLE folder. I saw the code input the DFT force into
lammps for MD energy minimization. What I need is not only the external
force but also lammps' own pair style force added together for energy
minimization. Would this be possible by using fix_external, instead of
fix_addforce ?

read the code and read it right. there is a point where asking too
many questions can get you into as much trouble as not asking any
questions. if you cannot see immediately the answer to your own
question, you should not be doing what you are trying to do.

axel.

Fix external does not get invoked by a minimization, though
you could edit the code to make it do that. However, you;d
also have to figure out how to make the energy of the
system consistent with the external forces (in a gradient sense).
Else minimization cannot work.

Steve

Fix external does not get invoked by a minimization, though
you could edit the code to make it do that. However, you;d

apparently "somebody" already did it about half a year ago:

8d815363 (sjplimp 2010-09-16 19:34:48 +0000 73) int
FixExternal::setmask()
8d815363 (sjplimp 2010-09-16 19:34:48 +0000 74) {
8d815363 (sjplimp 2010-09-16 19:34:48 +0000 75) int mask = 0;
fa8c7b68 (sjplimp 2013-01-04 17:15:32 +0000 76) if (mode ==
PF_CALLBACK || mode == PF_ARRAY) {
fa8c7b68 (sjplimp 2013-01-04 17:15:32 +0000 77) mask |= POST_FORCE;
fa8c7b68 (sjplimp 2013-01-04 17:15:32 +0000 78) mask |=
MIN_POST_FORCE;
fa8c7b68 (sjplimp 2013-01-04 17:15:32 +0000 79) }
8d815363 (sjplimp 2010-09-16 19:34:48 +0000 80) return mask;
8d815363 (sjplimp 2010-09-16 19:34:48 +0000 81) }

axel.

hmm … must have been a smart guy. But it
still won’t work unless you make energy and forces
consistent somehow.

Steve

Hi all,

Thanks for your suggestions! I have coupled my force codes with lammps EAM forces via fix_external. The only thing left is energy minimization part. Even though fix_external can invoke minimization(min_post_force) and total forces on the atom is updated correctly, lammps still know nothing about the extra energy due to my force. I tried setting eval as 0.0 in order to find local minimum only through force evaluation, but the pilot results seem wired, I will also keep on testing.

In the meanwhile I am wondering since I can compute the energy in a stand-alone code, whether it is possible to return the energy in the current version of lammps. Would this be possible only by modifying the driver program? Or I should extend the source code, such as fix_external.cpp and fix_modify.cpp ? Probably need to invoke mask |= MIN_ENERY" in the setmask ?

Best,

Xue

Hi all,

Thanks for your suggestions! I have coupled my force codes with lammps EAM
forces via fix_external. The only thing left is energy minimization part.
Even though fix_external can invoke minimization(min_post_force) and total
forces on the atom is updated correctly, lammps still know nothing about the
extra energy due to my force. I tried setting eval as 0.0 in order to find
local minimum only through force evaluation, but the pilot results seem
wired, I will also keep on testing.

well as a reference to compare to, you could run damped dynamics, i.e.
initialize velocities to some uniform distribution, and do an MD with
fix nve and fix viscous. that will not be the most efficient, but for
a well chosen viscosity parameter, it should get you pretty close to
the minimum that any of the minimizers will find. for complex systems,
it may even be competitive. some of the minimizers in LAMMPS are
similar. check out the docs.

In the meanwhile I am wondering since I can compute the energy in a
stand-alone code, whether it is possible to return the energy in the current
version of lammps. Would this be possible only by modifying the driver
program? Or I should extend the source code, such as fix_external.cpp and
fix_modify.cpp ? Probably need to invoke mask |= MIN_ENERY" in the setmask ?

i don't think so. this is only used by fix box relax. fix_modify
energy yes would be what i would do.
for that you need to implement compute_scalar() into fix_external, and
that should call another callback function from your code that returns
the energy instead of the forces. perhaps steve has a better
suggestion, but this is how i would go about it.

axel.

Hi Axel, and lammps-users,

Thanks for recommendation. I have been coupling GPU based calculation on my force with lammps lib, now it works well. Great!

Time to deal the energy minimization issue. Actually in my case, I do not need dynamics and thermo-effects. Pretty much a static case. My purpose is just to use min CG to find the stable configuration of the metallic cluster. I tried some runs, even though I set eval to be 0, still the minimization stops for energy evaluation(the EAM mechanical potential does not change much), and the forces are not below a specific standard. I may try to invoke only force evaluation for min, by commenting out the energy evaluation in the source code(min_cg) to see how things turn out.

I am wondering whether it would be sensible to use only force eval in this static case , would the results be reliable?

Still the best way should be returning the energy to lammps. I am working on this also. Appreciated if any one has more suggestions dealing with fix_external to return energy for minimize…
Thanks very much.

Xue

Hi Axel, and lammps-users,

Thanks for recommendation. I have been coupling GPU based calculation on my
force with lammps lib, now it works well. Great!

Time to deal the energy minimization issue. Actually in my case, I do not
need dynamics and thermo-effects. Pretty much a static case. My purpose is
just to use min CG to find the stable configuration of the metallic cluster.
I tried some runs, even though I set eval to be 0, still the minimization
stops for energy evaluation(the EAM mechanical potential does not change
much), and the forces are not below a specific standard. I may try to invoke
only force evaluation for min, by commenting out the energy evaluation in
the source code(min_cg) to see how things turn out.

I am wondering whether it would be sensible to use only force eval in this
static case , would the results be reliable?

if you do have the forces, why don't you run a damped dynamics as a
minimizer for starters?
just initialize the velocities to a reasonably high initial
temperature and dampen it down with this:
http://lammps.sandia.gov/doc/fix_viscous.html
you need to experiment a bit with the damping factor and you probably
need to play around with masses (it is usually often better to have
all masses be the same).
for some cases it can be quite competitive to other minimization
strategies and (potential) energy is only a diagnostic.

axel.

Or you could try min_style fire or quickmin which do a damped
dynamics. I think if you set the energy tolerance to 0.0,
they may work with only forces.

Steve