Improper_class2 angle mismatch between source code and doc page

LAMMPS LAMMPS Mailing List Mirror
1

Dear All,

For Improper_class2 force filed, we found a mismatch for the improper angle definition between the latest version source code and the doc page.

Please see the figure attached for the detail.

Since we are interested in calculating the force on atoms associated with Improper_class2 term, and we want to make comparison with Lammps results, we want to know which definition Lammps will choose in the calculation.

Currently, we are following the definition in doc page, so our result is different from Lammps.

Thank you for your help,

Lili Zhang

Improper_class2-angle mismatch.jpg

2

Hi Lili

Well, the documentation is wrong; the code, right. The “Wilson” out-of-plane definition is that defined in Wilson, E. B., Decius, J. C., and Cross, P. C. (1980) Molecular Vibrations (Dover, New York) p 59. and is clearly the angle between one bond and the plane defined by the other two bonds. This term is asymmetric with respect to the outer atoms, but as was clearly spelled out in the documentation of the original implementation in Discover (attached) or in the original paper on the Class 2 forcefields (Maple, Hagler & Dinur PNAS 85, 5350 (1988) ) the three possible out-of-plane terms are summed, leading to a symmetric result. It is always wise to refer to the original papers.

Paul

psaxe@…1795…

Phone: +1.760.495.4924 x 102
Cell: +1.505.603.8182

This email is covered by the Electronic Communications Privacy Act, 18 U.S.C. Section 2510-2521. This message and any attachments hereto may contain confidential information intended only for the use of the individual or entity named above. If you are not the intended recipient(s), or the employee or agent responsible for delivery of this message to the intended recipient(s), you are hereby notified that any dissemination, distribution or copying of this email message is strictly prohibited. If you have received this message in error, please immediately notify the sender and delete this email from your computer. The sender does not waive any privilege in the event this message was inadvertently disseminated.

MedeA® and Materials Design® are registered trademarks of Materials Design, Inc.

class2_oop_definition.pdf (163 KB)

3

Incidentally, on Feb 28th, John Jasa posted a notice saying he saw
some evidence of numerical problems using "improper class2". I don't
know if it's true, and I'm not trying to defame "class2". The example
he posted does not conserve energy, however this could happen if he
did not set the timestep carefully (In fact, it does not seem like he
set the timestep. I did not have time to look into this.) Anyway,
I'm curious to know if there is in fact a problem with "improper
class2".

I have included his example (see attachments), and forwarded his post
below. (For some reason I can't find it on the mail list archives.)

Cheers
Andrew
-------------John Jasa 's post from 2013-3-4 -------------------
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)

4

Hi Andrew

We are using the class 2 forcefield extensively and have seen no problems. Some time ago we also compared the energy and forces against Discover, which is the reference implementation of the class 2 forcefields, and found no problems.

We have seen a small problem with soft angles in notably O-Si-O that lead to 180º angles occurring occasionally during dynamics. This leads to a numerical issue causing the calculation to blow up. While there is code in LAMMPS to protect from this, it does not appear to work properly. We are looking into it. I will note that this seems to be a general feature of all the torsion terms as they all use roughly the same algorithm and same way to attempt to fix the problem when an angle goes to 180. So it does not appear to be specific to class 2 forcefields.

Paul

psaxe@…1795…

Phone: +1.760.495.4924 x 102
Cell: +1.505.603.8182

This email is covered by the Electronic Communications Privacy Act, 18 U.S.C. Section 2510-2521. This message and any attachments hereto may contain confidential information intended only for the use of the individual or entity named above. If you are not the intended recipient(s), or the employee or agent responsible for delivery of this message to the intended recipient(s), you are hereby notified that any dissemination, distribution or copying of this email message is strictly prohibited. If you have received this message in error, please immediately notify the sender and delete this email from your computer. The sender does not waive any privilege in the event this message was inadvertently disseminated.

MedeA® and Materials Design® are registered trademarks of Materials Design, Inc.

5

Dear Paul and Andrew,

Thank you for your help and explanation. John Jasa and I are working together. John is in charge of MD simulation, and I am working on multi-modeling.

Now, we are calculating force per atom from another arbitrary atom for class2 force field. We will report any issues during our investigation.

Thanks,

Lili