Possible bug in reported force constants after vibrational analysis

[Leopold Talirz]

Dear developers,

during an MD tutorial using CP2K I was looking at the force constants k 
reported after vibrational analysis. Below you find an example for the 
water molecule

VIB|                         1                    2                    3
 VIB|Frequency (cm^-1)  1589.789609          3637.026185         
 3741.860069
 VIB|Intensities           0.000000             0.000000             
0.000000
 VIB|Red.Masses (a.u.)     1.083688             1.044444             
1.081600
 VIB|Frc consts (a.u.)     0.017723             0.467889             
0.542862

I thought that the force constants are defined such that k = mu * omega^2, 
where mu is the reduced mass and omega the vibrational frequency (?). But 
this is obviously not the case.

Looking into vibrational_analysis.F , the formula 

    konst(i) = 
SIGN(1.0_dp,H_eigval2(i))*2.0_dp*pi**2*(ABS(H_eigval2(i))/massunit)**2*rmass(i)

for the force constant on line 406 (r14024) involves H_eigval2^2, which is 
proportional to omega^4 instead of omega^2.
Is this a bug?

Best,
Leopold



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20140523/9365fd77/attachment.htm>