Possible bug in reported force constants after vibrational analysis

[Matthias Krack]

Hi Leopold,

I think the code is correct, since the eigenvalues epsilon of the Hessian 
returned in H_eigval2 after its diagonalisation are equivalent to 
SQRT(konst/rmass). The frequencies, however, are proportional to 
SQRT(epsilon), i.e. SQRT(H_eigval2).

Matthias

Am Freitag, 23. Mai 2014 16:11:57 UTC+2 schrieb 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/20140525/646be5c4/attachment.htm>