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