function [freqs]=nwhess(hessfile,massfile)

%% 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);

%% Original incorrect form: scaling=(1/(2*pi*csi*t2au*sqrt(1837)))

%%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

%% 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

end
