439. Calculate frequencies from a hessian file from NWChem: example in Octave (matlab)

I wanted to calculate normal modes (frequencies) for specific atoms in a calculation, and so I had to write my own code.

This Octave code calculates frequencies for the first N atoms, where N is given in the input.mass file.

---

Background
The format that NWChem uses for the Hessian is that of a flat, triangular matrix i.e. a triangular matrix such as

1  
2 3 
4 5 6

is represented as

1
2
3
4
5
6

The Hessian is symmetric around the diagonal, so the full Hessian matrix is

1 2 4
2 3 5
4 5 6

The Hessian is independent of the masses of the atom pairs, while the frequencies are heavily dependent on the masses (isotope effects are quite visible for light elements).

To get the mass-weighted matrix we divide by the square root of the product of the masses (H * /(sqrt(m1*m2))). Note that the matrix reported in the nwchem output ("MASS-WEIGHTED NUCLEAR HESSIAN (Hartree/Bohr/Bohr/Kamu)") is multiplied by 1,000.

Once you have the mass-weighted hessian you need to calculate the eigenvalues, sort them and convert them to cm-1 using a scaling factor.

That's it.

---

The code:
See below for example input.mass and input.hess

%% prepare
clear;
format long

%%Calculate conversion factor from H/B/B/amu to cm-1

%% csi=299792458; %speed of light, m/s 
%% t2au=2.418884326505E-17; % seconds per a.u.
%% Better to do it by hand to avoid rounding errors:
cau=(2.99792458 * 2.418884326505)*1E-9; %c in metres per t(a.u.)

%% 1 electron (au)=9.10938291E-31 kg
%% 1 amu = 1.66053892E-27 kg
%% Better to do by hand to avoid rounding errors:
amu2au=(1.66053892/9.10938291)*1E4;% 1 amu in a.u. (via kgs)
%% For clarity
cmtom=1/100; %m per cm
%% And finally we get our scaling factor:
scaling=cmtom*(1/(2*pi*cau*sqrt(amu2au))); %( m/cm * 1/((m/au) * au) = m/cm * 1/m = 1/cm)


%%read masses
% The mass file contains the masses of the atoms
% The first line is the number of atoms in the file
% The remaining lines are the atom masses in the same order
% as the atoms are given in the nwchem input
protomasses=fopen("input.mass");
natoms=str2num(fgetl(protomasses));
for i = 1:natoms
 mass(end+1)=str2num(fgetl(protomasses));
end
fclose(protomasses);

%% Read and construct hessian from flat hessian in .hess file
%% The .hess file provided by nwchem is flat (i.e. one
%% dimensional) and is the triangular form (i.e half) of 
%% the full hessian. We use fgetl/str2num so that we can deal 
%% with instances of scientific notation in the hessian file.
%% While we"re at it we construct the mass-weighted force matrix too.
protohessian=fopen("input.hess"); 
hessian=zeros(3*natoms);
massweighted=zeros(3*natoms);

for i = 1:3*natoms
 for j=1:i
  hessian(i,j)=str2num(fgetl(protohessian));
  massweighted(i,j)=hessian(i,j)/sqrt( mass(ceil(i/3))*mass(ceil(j/3)));
 end
end

for i=1:3*natoms
 for j=1:i
  hessian(j,i)=hessian(i,j);
  massweighted(j,i)=massweighted(i,j);
 end
end

%% Diagonalize and compute frequencies in cm^{-1}
eigen=sort(eig(massweighted));
freqs=sqrt(eigen).*scaling;

%% Make imaginary frequencies negative and store them 
%% in a new array
for n=1:size(freqs,1)
 if imag(freqs(n))==0
  frequencies(end+1,1)=real(freqs(n));
 else
  frequencies(end+1,1)=-imag(freqs(n));
 end
end

%% Echo frequencies to stdout
printf("%10.4f \n",frequencies)
%% Save frequencies as well to modes.outs
outfile=fopen("normal.out","w");
fprintf(outfile,"%i \n",natoms);
fprintf(outfile,"%10.10f \n",frequencies);
fclose(outfile);
%save 'modes.out' -ascii  frequencies

---

input.mass (for water):

3
1.5994910D+01
1.0078250D+00
1.0078250D+00

---

input.hess (this one has imaginary frequencies as well):

     6.6177469151D-01
    -5.8658669668D-12
    -1.0013075598D-05
     1.0754299967D-09
     4.5060920407D-10
     3.6644723357D-01
    -3.3088202114D-01
     2.1099357839D-10
     1.6617441386D-01
     3.6163164885D-01
     2.5270659061D-12
     4.0920019206D-06
     3.2209366184D-11
     1.6382988861D-11
     8.3427731090D-07
     2.3904755566D-01
    -2.2311539742D-10
    -1.8322029567D-01
    -2.0261099118D-01
     1.1292349908D-10
     1.7796238990D-01
    -3.3088202212D-01
    -2.4469194991D-10
    -1.6617441477D-01
    -3.0749615389D-02
     1.2368245322D-10
    -3.6436594678D-02
     3.6163164980D-01
     2.5272503844D-12
     4.0920029550D-06
     3.2022391582D-11
     1.6289326095D-11
    -4.9229359909D-06
    -1.5407535297D-11
    -1.8816632580D-11
     8.3427660670D-07
    -2.3904755666D-01
    -2.2750774181D-10
    -1.8322029575D-01
     3.6436523006D-02
     1.1512005611D-10
     5.2580053385D-03
     2.0261099171D-01
     1.1238793371D-10
     1.7796238961D-01

---

Output:

 
  -11.0036 
   -1.6327 
    3.1676 
    3.9298 
    7.5811 
   12.2862 
 1619.0207 
 3616.0904 
 3781.1341

c.f.

 ----------------------------------------------------------------------------
 Normal Eigenvalue ||                 Infra Red Intensities
  Mode   [cm**-1]  || [atomic units] [(debye/angs)**2] [(KM/mol)] [arbitrary]
 ------ ---------- || -------------- ----------------- ---------- -----------
    1      -11.004 ||    0.426523           9.840       415.796      59.477
    2       -1.633 ||    0.000029           0.001         0.028       0.004
    3        3.168 ||    0.000003           0.000         0.003       0.000
    4        3.930 ||    0.000700           0.016         0.682       0.098
    5        7.581 ||    0.134394           3.101       131.014      18.741
    6       12.286 ||    0.000000           0.000         0.000       0.000
    7     1619.021 ||    0.070174           1.619        68.409       9.786
    8     3616.091 ||    0.004517           0.104         4.404       0.630
    9     3781.135 ||    0.009065           0.209         8.837       1.264
 ----------------------------------------------------------------------------