<div dir="ltr">Hi Leopold,<br><br>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).<br><br>Matthias<br><br>Am Freitag, 23. Mai 2014 16:11:57 UTC+2 schrieb Leopold Talirz:<blockquote class="gmail_quote" style="margin: 0;margin-left: 0.8ex;border-left: 1px #ccc solid;padding-left: 1ex;"><div dir="ltr">Dear developers,<div><br></div><div>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<br></div><div><br></div><div><div>VIB|                         1                    2                    3</div><div> VIB|Frequency (cm^-1)  1589.789609          3637.026185          3741.860069</div><div> VIB|Intensities           0.000000             0.000000             0.000000</div><div> VIB|Red.Masses (a.u.)     1.083688             1.044444             1.081600</div><div> VIB|Frc consts (a.u.)     0.017723             0.467889             0.542862</div></div><div><br></div><div>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.</div><div><br></div><div>Looking into vibrational_analysis.F , the formula </div><div><br></div><div>    konst(i) = SIGN(1.0_dp,H_eigval2(i))*2.0_<wbr>dp*pi**2*(ABS(H_eigval2(i))/<wbr>massunit)**2*rmass(i)</div><div><br></div><div>for the force constant on line 406 (r14024) involves H_eigval2^2, which is proportional to omega^4 instead of omega^2.</div><div>Is this a bug?</div><div><br></div><div>Best,</div><div>Leopold</div><div><br><div><br></div><div><br></div></div></div></blockquote></div>