Dear all,

We have found and fixed a bug in fix_shake when it is used to fix the bond lengths and angle in a three atom cluster. Below is a minimum pair of input files that reproduces the problem and a one-line change to fix_shake.cpp that fixed the bug.

If you have a three-atom cluster, X-Y-Z, where the bond lengths X-Y and Y-Z are fixed, and the angle X-Y-Z is fixed, this is implemented in LAMMPS as three distance constraints for X-Y, Y-Z and X-Z. Currently, the square of the X-Z distance is calculated in line 380 of fix_shake.cpp as follows:

rsq = 2.0*bond_distance[bond1_type]*bond_distance[bond2_type] *

(1.0-cos(angle));

This is correct only when the X-Y and Y-Z bond lengths are equal. The correct formula, irrespective of the bond lengths, is:

rsq = bond_distance[bond1_type]*bond_distance[bond1_type]

- bond_distance[bond2_type]*bond_distance[bond2_type]

- 2.0*bond_distance[bond1_type]*bond_distance[bond2_type]*cos(angle);

In the attached input files, we set up a single methanol-like molecule (with unified atoms), where the bond angle should be 108.5 degrees. If you run this input file, fix_shake sets that angle to 103.313 degrees. After implementing the above fix, the angle is correctly set to 108.5 degrees.

This bug was discovered by Lunna Li, who is a first-year graduate student in our group, and fixed jointly between us.

All the best,

Patrick Varilly

Post-doc in Daan Frenkel’s group at University of Cambridge

Two Files:

config-EC2.dat