Dear all,
I think I found a bug in the calculation of forces for the angle_style cosine/periodic when the multiplicity is one. In such a case (as a result of the fact that the loop at lines 133-137 is not entered), computed forces are only half of what they should be.
I suggest to modify line 120
un = 1.0;
into
if(m==1){un=2.0;}else{un=1.0;}
All code lines refer to the latest stable version (2Aug2023).
Best regards,
Paolo
When I input the posted potential formula into angle style lepton and use the angle_write command to compare the potential energy and force for the case n=1 and n=2 to that of angle style cosine/periodic, the potential energy and forces agree. Here is the input I am using:
atom_style molecular
units real
boundary p p p
region box block -10 10 -10 10 -10 10
create_box 1 box angle/types 1
mass 1 1.0
shell rm table1.txt table2.txt table3.txt table4.txt
angle_style lepton
angle_coeff 1 180 "2.0*C*(1-cos(theta)); C=75.0"
angle_write 1 100 table1.txt Lepton1
angle_coeff 1 180 "0.5*C*(1-cos(2*theta)); C=75.0"
angle_write 1 100 table2.txt Lepton2
angle_style cosine/periodic
angle_coeff 1 75 1 1
angle_write 1 100 table3.txt Periodic1
angle_coeff 1 75 1 2
angle_write 1 100 table4.txt Periodic2
And the two resulting plots
Correction: the angle_write command is not a sufficient test because it only uses the potential energy from the Angle::single() function and then generates the force through numerical differentiation.
Here is a proper test that, again, compares cosine/periodic against lepton and this one does show the factor 2 difference for the multiplicity = 1 case, so your suggested bugfix is required.
atom_style molecular
units real
boundary p p p
region box block -10 10 -10 10 -10 10
create_box 1 box angle/types 1 extra/angle/per/atom 2 extra/special/per/atom 4
mass 1 1.0
create_atoms 1 single 1.0 0.0 0.0
create_atoms 1 single 0.0 0.0 0.0
create_atoms 1 single -1.0 0.0 0.0
create_bonds single/angle 1 1 2 3
pair_style zero 1.0
pair_coeff 1 1
timestep 1.0
group move id 1
fix 1 move move rotate 0.0 0.0 0.0 0.0 0.0 1.0 100.0
thermo_style custom step pe eangle fmax
thermo 10
angle_style lepton
angle_coeff 1 0 "2.0*C*(1-cos(theta)); C=75.0"
dump 1 all custom 10 lepton-n=1.lammpstrj id x y z fx fy fz
run 100 post no
undump 1
reset_timestep 0
angle_coeff 1 0 "0.5*C*(1-cos(2*theta)); C=75.0"
dump 1 all custom 10 lepton-n=2.lammpstrj id x y z fx fy fz
run 100 post no
undump 1
reset_timestep 0
angle_style cosine/periodic
angle_coeff 1 75 1 1
dump 1 all custom 10 cosine-n=1.lammpstrj id x y z fx fy fz
run 100 post no
undump 1
reset_timestep 0
angle_coeff 1 75 1 2
dump 1 all custom 10 cosine-n=2.lammpstrj id x y z fx fy fz
run 100 post no
undump 1
reset_timestep 0
For m=1 I get with angle style lepton:
Step PotEng E_angle Fmax
0 0 0 0
10 28.647451 28.647451 159.49703
20 103.64745 103.64745 186.74237
30 196.35255 196.35255 142.65848
40 271.35255 271.35255 88.167788
50 300 300 5.9081468e-27
60 271.35255 271.35255 88.167788
70 196.35255 196.35255 142.65848
80 103.64745 103.64745 186.74237
90 28.647451 28.647451 159.49703
100 0 0 0
and with angle style cosine/periodic:
Step PotEng E_angle Fmax
0 0 0 9.6487359e-14
10 28.647451 28.647451 79.748513
20 103.64745 103.64745 93.371186
30 196.35255 196.35255 71.329239
40 271.35255 271.35255 44.083894
50 300 300 2.412184e-14
60 271.35255 271.35255 44.083894
70 196.35255 196.35255 71.329239
80 103.64745 103.64745 93.371186
90 28.647451 28.647451 79.748513
100 0 0 9.6487359e-14
For multiplicity m=2 the two variants agree. Here is angle style lepton data:
Step PotEng E_angle Fmax
0 0 0 0
10 25.911863 25.911863 129.03581
20 67.838137 67.838137 57.706566
30 67.838137 67.838137 44.083894
40 25.911863 25.911863 71.329239
50 0 0 5.9081468e-27
60 25.911863 25.911863 71.329239
70 67.838137 67.838137 44.083894
80 67.838137 67.838137 57.706566
90 25.911863 25.911863 129.03581
100 0 0 0
and here angle style cosine/periodic:
Step PotEng E_angle Fmax
0 0 0 1.9297472e-13
10 25.911863 25.911863 129.03581
20 67.838137 67.838137 57.706566
30 67.838137 67.838137 44.083894
40 25.911863 25.911863 71.329239
50 0 0 4.8243679e-14
60 25.911863 25.911863 71.329239
70 67.838137 67.838137 44.083894
80 67.838137 67.838137 57.706566
90 25.911863 25.911863 129.03581
100 0 0 1.9297472e-13
Dear Axel,
thank you for your prompt and thorough reply. I hope this bugfix can be useful.
Best regards,
Paolo