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 6is 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
----------------------------------------------------------------------------
