# How to define a piecewise-continuous force?

In my system, I want to add a force to the z component for each atom, changing as a function of radial distance r. In a simple example:

1/r r < r_c
0 r >= r_c

How do I define such a piecewise-continuous force as a function of the radial distance so that it is compatible with the addforce command? Thanks.

In my system, I want to add a force to the z component for each atom, changing as a function of radial distance r. In a simple example:

1/r r < r_c
0 r >= r_c

How do I define such a piecewise-continuous force as a function of the radial distance so that it is compatible with the addforce command? Thanks.

you can use logical expressions in atom style and equal style expressions. here you need to use, of course, an atom style variable. those logicals would evaluate to either 1.0 (if true) or 0.0 (if false) and thus if you multiply those with the force contribution for a given section and sum each segment, you can define a piecewise force.

axel.

Hi Samuel.
I have nothing useful to add except that force you described has a non-zero curl. LAMMPS does not mind, but your potential energy numbers will be meaningless. (But this is true whenever you use fix addforce without the “energy” keyword.)
Don’t know if that matters for your application.
Cheers

-andrew

I ended up adding a force in z (see below) that should be non-zero only below a certain cutoff. However, when I run the simulation — inserting a test particle that interacts with a substrate — I see that the energy as given by evdwl for the test particle doesn’t go to zero past this cutoff. Why not? Also, the energy is extremely large despite that I thought it should be just the energy for the one test particle. How to correct this?

Thanks again.

variable Fz atom “(z<9.033713001)(-0.13310(-(6/5)(3.011237667^9/z^10)+3(3.011237667^3/z^4))) + (z>=9.033713001)*0.0”

variable Ez atom “(z<9.033713001)(0.13310((2/15)*(3.011237667/z)^9-(3.011237667/z)^3)) + (z>=9.033713001)*0.0”

fix 4 particle addforce 0.0 0.0 v_Fz energy v_Ez

LAMMPS will compute and do what you ask it to do. there are no know bugs for any of the features that you are using. so if you don’t get out what you expect, you need to debug/check that would you put in is what you expect. since this is for a single particle, this should be easily computed manually and you can output the numbers that LAMMPS is computing via a custom dump.

axel.

How do I look at just the energy of the test particle with the added force rather than the energy of the entire system?

1. variable Fz atom “(z<9.033713001)(-0.13310(-(6/5)(3.011237667^9/z^10)+3(3.011237667^3/z^4))) + (z>=9.033713001)*0.0”

2. variable Ez atom “(z<9.033713001)(0.13310((2/15)*(3.011237667/z)^9-(3.011237667/z)^3)) + (z>=9.033713001)*0.0”

3. fix 4 testparticle addforce 0.0 0.0 v_Fz energy v_Ez

4. fix_modify 4 energy yes

5. thermo_style custom step temp epair emol ecoul evdwl

you can output atom style variables in custom dump files.

axel.

Hi Samuel.
I have nothing useful to add except that force you described has a non-zero curl. LAMMPS does not mind, but your potential energy numbers will be meaningless. (But this is true whenever you use fix addforce without the “energy” keyword.)

But this is not why your energy numbers will be meaningless. Any external force that only acts in the z direction but depends also on the x and y position will not conserve energy. You have created a perpetual motion machine. Perhaps that’s okay. I don’t know if this detail matters for your application… (Why do you need to compute the energy?)

Anyway, once you have a conservative force with a meaningful expression for energy, then you can add it to the total energy of the system using “fix_modify energy…”, as you are doing.

To make a lazy attempt to answer your question: In that case, I suspect you could calculate the sum of the energies on your particles using something like “compute Eext all reduce sum v_Ez”, and print them to a file using “dump” with the “custom” keyword, and “c_Eext”). (Or use the “print” command.) There might be other ways to calculate it as well.

https://lammps.sandia.gov/doc/compute_reduce.html

(Forgive me for being terse. I haven’t tried these commands, and I might have got the syntax wrong. I’m using my phone to reply to your post.)

I hope this helps.
-andrew