Class 2 Improper Force Calculations

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

My team and I are working on obtaining the per-atom forces using a post-processing program after LAMMPS is run. To do this, we took the derivatives of the energy formulas posted on each part of the doc pages to obtain equations for forces, then used the forcefield coefficients and outputted positions to calculate forces. Specifically, we are using class2 for all types.

Using small nonperiodic systems of a few atoms we have verified that the pair, bond, angle, and dihedral force contributions are working as intended by comparing them to the LAMMPS output for forces. However, when an improper is introduced into the system, the forces are not computed correctly via the program. We have error-checked it many times over now, but the problem still exists. I figure there’s probably something small wrong with it, but in the meantime I’m checking all possible sources for mistakes.

My question is concerning the improper force calculations done in LAMMPS. Is there any chance there’s more going on than the improper and angleangle terms that we perhaps missed? Following the doc pages should encompass everything, but I fear there is a missed term or some contribution not accounted for.

We are looking for a way to examine the per-atom forces on a particular subset of atoms while cutting off relations to other outside the set, and using this post-processing tool seemed to be the best mode of action since the suggested commands in LAMMPS could not achieve what we wanted.

Thanks for your time,
John Jasa
Mechanical and Materials Engineering
University of Nebraska-Lincoln

2

Hi all,

My team and I are working on obtaining the per-atom forces using a
post-processing program after LAMMPS is run. To do this, we took the
derivatives of the energy formulas posted on each part of the doc pages to
obtain equations for forces, then used the forcefield coefficients and
outputted positions to calculate forces. Specifically, we are using class2
for all types.

Using small nonperiodic systems of a few atoms we have verified that the
pair, bond, angle, and dihedral force contributions are working as intended
by comparing them to the LAMMPS output for forces. However, when an
improper is introduced into the system, the forces are not computed
correctly via the program. We have error-checked it many times over now,
but the problem still exists. I figure there's probably something small
wrong with it, but in the meantime I'm checking all possible sources for
mistakes.

it is difficult to discuss this without an example demonstrating the issue.

My question is concerning the improper force calculations done in LAMMPS. Is
there any chance there's more going on than the improper and angleangle
terms that we perhaps missed? Following the doc pages should encompass
everything, but I fear there is a missed term or some contribution not
accounted for.

the only definite answer lies in the source code.

We are looking for a way to examine the per-atom forces on a particular
subset of atoms while cutting off relations to other outside the set, and
using this post-processing tool seemed to be the best mode of action since
the suggested commands in LAMMPS could not achieve what we wanted.

huh? what suggested commands? you have to go a bit more into detail.
there are a *lot* of things that can be done through creative use of
the existing commands. not all of that may be obvious or explicitly
documented.

axel.

3

Hi there,

Copied below is the data file for a simple test system used:

LAMMPS 2005 data file for TinyForce

4 atoms
0 bonds
0 angles
0 dihedrals
1 impropers

1 atom types
1 improper types

-10 10 xlo xhi
-10 10 ylo yhi
-10 10 zlo zhi

Masses

1 12.011150

Pair Coeffs

1 0.0540000000 4.0100000000

Improper Coeffs

1 1000 180

AngleAngle Coeffs

1 550 342 645 131 151 141

Atoms

1 1 1 0.000000 0 0 0
2 2 1 0.000000 0 2 0
3 3 1 0.000000 2 0 0
4 4 1 0.000000 0 0 2

Impropers

1 1 1 2 3 4

It’s obviously quite unrealistic, but I am using it to diagnose the problem with force computations. Here, the forces per atom as computed from LAMMPS should match the post-processing program used, but they do not.

I’ll scour through the source code to make sure that our calculations coincide with what LAMMPS is doing.

We are trying to compute only the forces on atoms inside of a region due to atoms within it. As such, it was suggested we group the atoms outside the region and the ones inside and turn off pairwise interactions, but that would not account for the bond, angle, and dihedral interactions that straddle the boundary. By computing the contributions from each interaction individually, we would be able to select which are included for the region, our end goal. If you have any suggestions about how to go about doing this, they’d be much appreciated!

Thank you for your help,
John Jasa
Mechanical and Materials Engineering
University of Nebraska-Lincoln

4

Just a quick email:
First of all, thanks for reporting this.
If there is a bug in the improper_class2.cpp file, these tests might
help locate it:

1) Is energy conserved?
To make it easier for other people to observe the same error you are
running into, try creating a simple DATA file, and run LAMMPS using
"fix nve". Is the total energy conserved? I am guessing the answer
is no. That would be a very useful test. (I attached a small DATA
file you can modify. See below.)

2) Does this problem occur when you use improper-energy functions
which are even?
(Try setting the "X0" parameter to 0. Does the problem still occur?)

3) Do your atoms lie well-within within the periodic boundary box?
(Try increasing the periodic boundary box size in all directions.)

----- small test case -----
I doubt this is particularly helpful, but I have a attached some tiny
input files on a simple 4-atom system. (There was a bug in the
dihedral_class2.cpp, and I used these files to find that bug.) Feel
free to modify the example to use "improper_style class2" instead of
"dihedral_style class2". Also, we want to use NVE conditions, so
please edit the "in.test" file and comment out the line that says "fix
1 all langevin 300 1.0 100.0 48279". (That line was used to anneal
the conformation to the lowest energy. But here we want to conserve
energy.)

I think these are the only modifications that you have to make to this
file. If you like, you can move the 4 atoms around and see if the
problem goes away for certain angles (positive or negative).

Thank you for reporting the error. There have been errors like this
in the past. If you have time to create a simple test case for us,
that would be helpful.
Cheers!

Andrew

P.S.
---OPTIONAL CHANGES:---
The molecule used in that example is a linear chain of 4 beads.
(1-2-3-4) The "Bonds" section (in "butane.data") is:

Bonds

    1 1 1 2
    2 1 2 3
    3 1 3 4

Although LAMMPS does not care about the bonds if you like, you can
change the bond topology to resemble a branching hub (like "BF3", or
boron hydride)

  2
  >
3-1-4

To do that, replace the "Bonds" section with something like this:

Bonds

    1 1 1 2
    2 1 1 3
    3 1 1 4

Cheers!

test_system_class2.tar.gz (1.47 KB)

5

oops, you posted an example before I finished my email.

6

Can you run "fix nve" on this for us and tell us if energy is conserved?
(That will help us rule out the possibility there might be an error in
your post-processing script. -Just a possibility.)
Thanks

Andrew

7

the only definite answer lies in the source code.

yes - why not look in src/improper_class2.cpp. Everything
is there in one place.

Steve

8

Hi there,

Thanks much for the help! Going through your steps, I got this output for NVE for the simple system with only pair and a single improper:

Step Temp E_pair E_mol TotEng Press Volume
0 75.318871 -0.0053814401 -1157.3088 -1156.6407 15.409808 9856.1106
1000 1203.1993 -0.050567446 -1141.1285 -1130.4195 66.081203 11749.538
2000 3280.8694 0 -1146.7563 -1117.4174 22.308446 67552.829
3000 3381.1489 0 -1143.8259 -1113.5902 5.6520083 274604.02
4000 3872.6671 0 -1152.4135 -1117.7824 2.5031675 676718.7
5000 4235.652 0 -1156.9435 -1119.0665 1.1762173 1562098.8
6000 3785.9493 0 -1152.8176 -1118.962 0.54918059 2999107.4
7000 3757.9019 0 -1152.2252 -1118.6204 0.3276495 5105683.9
8000 4060.3406 0 -1155.1911 -1118.8817 0.21769128 8038926.2
9000 4089.5768 0 -1155.6786 -1119.1078 0.14488222 11884200
10000 4402.8871 0 -1158.612 -1119.2395 0.10896606 16853927

So the total energy appears to start at -1156.6407 then come to rest around a value of -1119.2395. This is with arbitrarily chosen coefficients for each of the relations.

Changing the value of chi0 to 0 yielded:

Step Temp E_pair E_mol TotEng Press Volume
0 253.32252 -0.022262517 -1593.7227 -1591.4796 5.8869811 237953.25
1000 3431.7074 0 -1597.3416 -1566.6538 4.6246316 309276.67
2000 3164.4039 -0.0045368486 -1595.3152 -1567.0223 3.222129 405380.67
3000 3151.1313 0 -1596.2679 -1568.0891 1.8733732 693272.88
4000 3305.3063 0 -1590.8601 -1561.3027 1.1017327 1226781.4
5000 3728.5643 0 -1598.1024 -1564.7599 0.71864144 2126120.3
6000 3462.4435 0 -1592.3131 -1561.3504 0.40161397 3574892.5
7000 3500.4681 0 -1592.7041 -1561.4014 0.23239199 6239732.8
8000 3913.3955 0 -1597.0418 -1562.0465 0.17765942 9044915.8
9000 3696.002 0 -1594.767 -1561.7158 0.11508685 13520692
10000 3931.9689 0 -1596.6941 -1561.5327 0.085038176 19220146

Again, the total energy starts at a certain amount (-1591.4796) before leveling off around -1561.

Yes, the atoms are well within the boundaries of the system, however for this simple test case I’m using shrink-wrapped non-periodic boundaries to validate the forces felt by only the atoms in one image. Should I only use periodic boundary conditions in this case?

I sent my code and the source code to a colleague so hopefully another pair of eyes can help shed light on possible discrepancies or errors. Attached are the files used for the force investigation.

Thanks again,
John Jasa

in.TinyForce (402 Bytes)

TinyForce.data (741 Bytes)