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
----------------------------------------------------------------------------
Do you have any information with regards to computational cost? As compared to the standard calculation?
ReplyDeleteThere's no noticeable computational cost using the script above. All the computational work was done generating the Hessian. I presume that by standard calculation you mean letting NWChem compute the frequencies instead?
DeleteCurious if you were able to match the eigenvectors of these frequencies with what NWChem outputs? This was extremely helpful btw!!!
ReplyDelete