We (meaning my student Jacob Breese) noticed that in some simulations using the angle_style cosine/delta, many of the bonds would tend to straighten out towards 180 degrees, even though that was not the angle setting and the force constant was high. This is in the current 29 Aug 2024 version but was also present in a prior versions.
Unfortunately, this issue was brought up in 2018 here Angle_Style cosine/delta but we had not noticed that post until I went to post this note. That thread came to the same conclusion as below, but without details.
To get around the issue, one can use the alternative funciton angle_style cosine/shift from the EXTRA-MOLECULE package, which seems to work correctly, except note that that style does not include the 1/2 factor in the energy constant, and the energy is shifted such that its maximum is at 0 rather than its minimum.
The issue, it seems, is that the inverse of a sin term was accidentally taken twice in angle_cos_delta.cpp, leading to a part of the force calculation to include sin(theta) instead of 1/sin(theta)
Line 92-94, where s is sin(theta):
s = 1.0/s;
cot = c/s;
Deleting line 92 (and line 131 in angle_cos_delta_omp.cpp) or making line 94 cot=c*s; should fix the problem and bring the calculation there to equal that in angle_cosine_shift.cpp. I verified that angle_cosine_shift.cpp appears correct with my pencil and paper calculations, though not by printing any intermediate results from the code itself.
Based on this, we can see that, if the angle setting (theta0) is 90 degrees and the force constant is large such that the bond stays near 90, then the function will perform well anyway. The further you get from 90, the further the force will get from the expected value (the difference is not directly proportional because there are two terms and the error only impacts one of them). At first, the force curves near the equilibrium point are similar to the expected curves, but shifted to create an equilibrium angle somewhat further from 90 than the setting (if the setting is > 90 it is larger, or if <90 it is smaller). At angle setting > about 116 degrees (or similarly low angles, but I am not sure anyone would want that), the force curves are qualitatively different than expected. This is based on simply comparing the functions calculated as shown in the file, but we also noticed in testing that a setting of 116 degrees does not lead to angles near 116 degrees.
As context, we typically simulate coarse-grained polymers where there is a history of using a cos(theta) angle potential (with equilibrium bond angle of 180 degrees but a low force constant) to stiffen the polymers in a generic way without impacting other properties too much. That’s why when we wanted to set a specific angle instead, we used the similar potential cos(theta-theta0) to compare more easily to prior work, whereas others may have used a harmonic potential for this purpose. I don’t know how many people may be using these angle styles for other reasons.
@halllm Thanks for your detailed and thorough bug report. I’ve applied the corresponding changes to pull requests that will update both, the stable and the development version of LAMMPS.
1 Like
Thanks very much, appreciate everything you do for this community!