Using TAMkin with CP2K - how is CP2K diagonalizing the Hessian matrix?

[Torstein Fjermestad]

Dear all, 

I would like to use TAMkin (http://molmod.github.io/tamkin/) to compute thermodynamic properties from the vibrational frequencies computed with CP2K.

As I have understood it, the way TAMkin works is to read and diagonalize the Hessian matrix printed in the CP2K output.

When I compare the frequencies obtained with TAMkin with those obtained with cp2k, the values are not the same (Compare column 2 and 3 in the table below). 

The first task is then to understand how CP2K computes the frequencies. To investigate this, I copied the Hessian matrix from the CP2K output (printed below the string "Hessian in cartesian coordinates") and pasted it into Excel and saved it as a csv file. This file I read into python and diagonalized the Hessian with numpy. The frequencies I get are different from those printed in the cp2k output file (compare column 3 and 4 in the  table below). 

I have assumed that the units of the Hessian matrix in the cp2k output file are Ry*bohr-2*kamu-1 (kamu = 1000*amu). This assumption is based on trial and error; it was what got me closest to the cp2k frequencies. 

Table 1: The eight lowest frequencies. There are 48 in total.

frequencies	TAMkin / cm-1	cp2k / cm-1	Diag., cp2k Hessian / cm-1
1	31.1	27.3	32.6
2	38.3	40.8	40.1
3	47.7	51.0	50.2
4	74.7	74.0	78.2
5	82.0	80.0	85.8
6	111.0	109.9	116.3
7	223.7	224.1	234.3
8	245.6	246.3	257.3

What is wrong with my calculation?

What steps should be done to go from the Hessian to the frequencies?

The Hessian in csv format, the python script and the cp2k output file are attached. 

Thanks a lot for your help.

Regards,

Torstein Fjermestad

 

[Torstein Fjermestad]

Dear all, 

I would like to follow up on this question that I asked on June 10 this year.

I realize that I might have expressed myself in a very complicated way, and I therefore try to be clearer this time.

The problem is that when I diagonalize the Hessian printed in the cp2k output file with the following python code (The complete python code is attached)

from numpy import genfromtxt
import numpy as np
import sys

hessianFile=sys.argv[1]
hessian = genfromtxt(hessianFile,delimiter=',')
freq = np.linalg.eigvals(hessian)

 

, I get slightly different frequencies than what is printed in the output file (and the .mol file). The frequencies are compared in the table below.

Column 2: Frequencies from CP2K output file

Column 3: Frequencies obtained when diagonalizing the Hessian with numpy

Table 1: The eight lowest frequencies. There are 48 in total.

What am I doing wrong in the diagonalization?

[Krack Matthias (PSI)]

Dear Torstein

 

The Hessian printed after "Hessian in cartesian coordinates" is the Hessian within the Cartesian coordinate system of the input cell. Only the first (maximum) 9 “VIB| Cartesian Low frequencies” are printed for that Hessian. The eigenvalues printed thereafter as well as in the molden file are obtained from the Hessian after rotations and translations have been removed (after transformation with the D matrix). I think that the implementation follows more or less the approach described here.

 

HTH

 

Matthias

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e7e9984b-2035-4629-a64e-a896589aeba2n%40googlegroups.com.

[Eugene]

This is defenetly late responce but in order to keep information there is a reason for responce. I don't know units in CP2K output, but CP2K save Hessian in unformatted bynary file, which can be easily read and then diagonalized. The sqrt of eigenvalues is freqeuncy of vibrational modes wich can be converted to cm-1 by factor multiplying, where the factor is:

sqrt(Hartree/AMU/Bohr2Angstrom)*1/2/pi/Angstrom/c*sqrt(massunit/Bohr2Angstrom)/1)10**5

Hartree=4.359748536313086e-18

AMU=1.6605402e-27

Bohr2Angstrom=0.529177207423948

massunit=a_mass/e_mass

a_mass=1.660538782e-27

e_mass=9.10938215e-31

I've wrote convertor from cp2k hessian file format to FORCE_CONSTANTS file format (see attachment). Using PHONOPY one can write dynamical matrix, plot eigenvectors and lot of other things. 

Best wishes,
Eugene

среда, 10 июня 2020 г. в 16:56:53 UTC+3, tfjer...@gmail.com: