I've been working with the spring constraints, and the numbers seemed to be coming out a bit wrong.
Example output:
adding spring # 1
spring parameters (i,j,k,r0): 1 12
6.420000000000000E-006 380.000000000000
spring length and energy : 9.80868047825847 0.879807155776173
spring forces : -4.753256542659162E-003 0.000000000000000E+000
0.000000000000000E+000
spring deriv : 4.753256542659162E-003
Doing the math to verify the energy is calculated correctly:
delta_r = r0 - r = 370.19
E = 0.5 * k * (delta_r ^ 2) = 0.4399
Above it's stated as 0.879, or double what the energy should be. So I figured I would look at the code.
In the file src/cons/cons_springs.F,
the subroutines are cons_spring_energy and cons_spring_force
Line 491 (and at line 461):
Line 491 and 461:
energy = k * (r - r0)**2
Line 464:
do i=1,3
f(i)=2*k*(r-r0)*
& (r1(i)-r2(i))/r
end do
So, it seems (because it was put in twice), that either it was copy-pasted blindly or that someone decided it was correct to double the energy of the spring. The only reason I can think of doubling the energy of the spring for the calculation is to account for the energy of the stretched molecule as well as the spring. Because the deformed-geometry energy will be acounted for already when the SCF is calculated (task energy), this seems redundant and incorrect.
Oddly enough, on line 464 (also shown above), the spring constant is doubled for calculating the force. This at least keeps the force consistent with the energy, which is nice and makes it usable. For now I'll just divide the spring constant I want to use by 2 before inputting it.
Let me know if I'm interpreting this wrong, or if there is some other reason for this behavior.
|