Simulated infrared spectra

Simulated infrared spectra#

At http://homepage.univie.ac.at/david.karhanek/downloads.html#Entry02 there is a recipe for computing the Infrared vibrational spectroscopy intensities in VASP. We are going to do that for water here. First, we will relax a water molecule.

Note

that link seems to be dead.

from ase import Atoms, Atom
from vasp import Vasp

atoms = Atoms([Atom('H', [0.5960812,  -0.7677068,   0.0000000]),
               Atom('O', [0.0000000,   0.0000000,   0.0000000]),
               Atom('H', [0.5960812,   0.7677068,   0.0000000])],
              cell=(8, 8, 8))

calc = Vasp(label='molecules/h2o_relax',
          xc='PBE',
            encut=400,
            ismear=0,  # Gaussian smearing
            ibrion=2,
            ediff=1e-8,
            nsw=10,
            atoms=atoms)
print('Forces')
print('===========================')
print(atoms.get_forces())

Forces
===========================
[[-3.813e-05  5.328e-05  0.000e+00]
 [ 7.625e-05  0.000e+00  0.000e+00]
 [-3.813e-05 -5.328e-05  0.000e+00]]

Next, we instruct VASP to compute the vibrational modes using density functional perturbation theory with IBRION=7. Note, this is different than in Vibrational frequencies where finite differences were used.

from vasp import Vasp
from ase.io import read
import os

# read in relaxed geometry
contcar = 'molecules/h2o_relax/CONTCAR'
if os.path.exists(contcar):
    atoms = read(contcar, format='vasp')
    print(f'Loaded H2O with {len(atoms)} atoms')
    print(f'Positions:\n{atoms.positions}')
else:
    # Use molecule builder as fallback
    from ase.build import molecule
    atoms = molecule('H2O')
    atoms.center(vacuum=6)
    print('Using ASE molecule builder (no relaxed geometry found)')

! rm -fr molecules/h2o_vib_dfpt
# Define vibration calculation (requires VASP to actually run)
calc = Vasp(label='molecules/h2o_vib_dfpt',
            xc='PBE',
            encut=400,
            ibrion=7,  # DFPT
            nfree=2,
            potim=0.015,
            nsw=1,
            lepsilon=True,
            atoms=atoms)

calc.calculate()

Loaded H2O with 3 atoms
Positions:
[[ 5.96211047e-01  7.23151102e+00  0.00000000e+00]
 [ 5.96211047e-01  7.68488978e-01  0.00000000e+00]
 [-2.59694976e-04  0.00000000e+00  0.00000000e+00]]

To analyze the results, this shell script was provided to extract the results.

#!/bin/bash
# A utility for calculating the vibrational intensities from VASP output (OUTCAR)
# (C) David Karhanek, 2011-03-25, ICIQ Tarragona, Spain (www.iciq.es)

# extract Born effective charges tensors
printf "..reading OUTCAR"
BORN_NROWS=`grep NIONS OUTCAR | awk '{print $12*4+1}'`
if [ `grep 'BORN' OUTCAR | wc -l` = 0 ] ; then \
   printf " .. FAILED! Born effective charges missing! Bye! \n\n" ; exit 1 ; fi
grep "in e, cummulative" -A $BORN_NROWS OUTCAR > born.txt

# extract Eigenvectors and eigenvalues
if [ `grep 'SQRT(mass)' OUTCAR | wc -l` != 1 ] ; then \
   printf " .. FAILED! Restart VASP with NWRITE=3! Bye! \n\n" ; exit 1 ; fi
EIG_NVIBS=`grep -A 2000 'SQRT(mass)' OUTCAR | grep 'cm-1' | wc -l`
EIG_NIONS=`grep NIONS OUTCAR | awk '{print $12}'`
EIG_NROWS=`echo "($EIG_NIONS+3)*$EIG_NVIBS+3" | bc`
grep -A $(($EIG_NROWS+2)) 'SQRT(mass)' OUTCAR | tail -n $(($EIG_NROWS+1)) | sed 's/f\/i/fi /g' > eigenvectors.txt
printf " ..done\n"

# set up a new directory, split files - prepare for parsing
printf "..splitting files"
mkdir intensities ; mv born.txt eigenvectors.txt intensities/
cd intensities/
let NBORN_NROWS=BORN_NROWS-1
let NEIG_NROWS=EIG_NROWS-3
let NBORN_STEP=4
let NEIG_STEP=EIG_NIONS+3
tail -n $NBORN_NROWS born.txt > temp.born.txt
tail -n $NEIG_NROWS eigenvectors.txt > temp.eige.txt
mkdir inputs ; mv born.txt eigenvectors.txt inputs/
split -a 3 -d -l $NEIG_STEP temp.eige.txt temp.ei.
split -a 3 -d -l $NBORN_STEP temp.born.txt temp.bo.
mkdir temps01 ; mv temp.born.txt temp.eige.txt temps01/
for nu in `seq 1 $EIG_NVIBS` ; do
 let nud=nu-1 ; ei=`printf "%03u" $nu` ; eid=`printf "%03u" $nud` ; mv temp.ei.$eid eigens.vib.$ei
done
for s in `seq 1 $EIG_NIONS` ; do
 let sd=s-1 ; bo=`printf "%03u" $s` ; bod=`printf "%03u" $sd` ; mv temp.bo.$bod borncs.$bo
done
printf " ..done\n"

# parse deviation vectors (eig)
printf "..parsing eigenvectors"
let sad=$EIG_NIONS+1
for nu in `seq 1 $EIG_NVIBS` ; do
 nuu=`printf "%03u" $nu`
 tail -n $sad eigens.vib.$nuu | head -n $EIG_NIONS | awk '{print $4,$5,$6}' > e.vib.$nuu.allions
 split -a 3 -d -l 1 e.vib.$nuu.allions temp.e.vib.$nuu.ion.
 for s in `seq 1 $EIG_NIONS` ; do
  let sd=s-1; bo=`printf "%03u" $s`; bod=`printf "%03u" $sd`; mv temp.e.vib.$nuu.ion.$bod e.vib.$nuu.ion.$bo
 done
done
printf " ..done\n"

# parse born effective charge matrices (born)
printf "..parsing eff.charges"
for s in `seq 1 $EIG_NIONS` ; do
 ss=`printf "%03u" $s`
 awk '{print $2,$3,$4}' borncs.$ss | tail -3 > bornch.$ss
done
mkdir temps02 ; mv eigens.* borncs.* temps02/
printf " ..done\n"

# parse matrices, multiply them and collect squares (giving intensities)
printf "..multiplying matrices, summing "
for nu in `seq 1 $EIG_NVIBS` ; do
 nuu=`printf "%03u" $nu`
 int=0.0
 for alpha in 1 2 3 ;  do            # summing over alpha coordinates
  sumpol=0.0
  for s in `seq 1 $EIG_NIONS` ; do   # summing over atoms
   ss=`printf "%03u" $s`
   awk -v a="$alpha" '(NR==a){print}' bornch.$ss > z.ion.$ss.alpha.$alpha
   # summing over beta coordinates and multiplying Z(s,alpha)*e(s) done by the following awk script
   paste z.ion.$ss.alpha.$alpha  e.vib.$nuu.ion.$ss | \
   awk '{pol=$1*$4+$2*$5+$3*$6; print $0,"  ",pol}' > matr-vib-${nuu}-alpha-${alpha}-ion-${ss}
  done
  sumpol=`cat matr-vib-${nuu}-alpha-${alpha}-ion-* | awk '{sum+=$7} END {print sum}'`
  int=`echo "$int+($sumpol)^2" | sed 's/[eE]/*10^/g' |  bc -l`
 done
 freq=`awk '(NR==1){print $8}' temps02/eigens.vib.$nuu`
 echo "$nuu $freq $int">> exact.res.txt
 printf "."
done
printf " ..done\n"

# format results, normalize intensities
printf "..normalizing intensities"
max=`awk '(NR==1){max=$3} $3>=max {max=$3} END {print max}' exact.res.txt`
awk -v max="$max" '{printf "%03u %6.1f %5.3f\n",$1,$2,$3/max}' exact.res.txt > results.txt
printf " ..done\n"

# clean up, display results
printf "..finalizing:\n"
mkdir temps03; mv bornch.* e.vib.*.allions temps03/
mkdir temps04; mv z.ion* e.vib.*.ion.* temps04/
mkdir temps05; mv matr-* temps05/
mkdir results; mv *res*txt results/
let NMATRIX=$EIG_NVIBS**2
printf "%5u atoms found\n%5u vibrations found\n%5u matrices evaluated" \
       $EIG_NIONS $EIG_NVIBS $NMATRIX > results/statistics.txt
  # fast switch to clean up all temporary files
  rm -r temps*
cat results/results.txt

Note that the results above include the rotational and translational modes (modes 4-9). The following shell script removes those, and recalculates the intensities. Note that it appears to just remove the last 6 modes and req compute the intensities. It is not obvious that will always be the right way to do it as the order of the eigenvectors is not guaranteed.

#!/bin/bash
# reformat intensities, just normal modes: 3N -> (3N-6)
printf "..reformatting and normalizing intensities"
cd intensities/results/
nlns=`wc -l exact.res.txt | awk '{print $1}' `; let bodylns=nlns-6
head -n $bodylns exact.res.txt > temp.reform.res.txt
max=`awk '(NR==1){max=$3} $3>=max {max=$3} END {print max}' temp.reform.res.txt`
awk -v max="$max" '{print $1,$2,$3/max}' temp.reform.res.txt > exact.reform.res.txt
awk -v max="$max" '{printf "%03u %6.1f %5.3f\n",$1,$2,$3/max}' temp.reform.res.txt > reform.res.txt
printf " ..done\n..normal modes:\n"
rm temp.reform.res.txt
cat reform.res.txt
cd ../..

..reformatting and normalizing intensities ..done
..normal modes:

..reformatting and normalizing intensities ..done
..normal modes:

The interpretation of these results is that the mode at 3713 cm^-1 would be nearly invisible in the IR spectrum. Earlier we interpreted that as the symmetric stretch. In this mode, there is only a small change in the molecule dipole moment, so there is a small IR intensity.

See also cite:giannozzi:8537. For HREELS simulations see cite:0953-8984-22-26-265006.

The shell script above has been translated to a convenient python function in mod:vasp.

from vasp import Vasp
from ase.io import read
import os

contcar = 'molecules/h2o_vib_dfpt/CONTCAR'
outcar = 'molecules/h2o_vib_dfpt/OUTCAR'

if os.path.exists(outcar):
    atoms = read(contcar, format='vasp')
    calc = Vasp(label='molecules/h2o_vib_dfpt', atoms=atoms)
    try:
        print('mode  Relative intensity')
        for i, intensity in enumerate(calc.get_infrared_intensities()):
            print('{0:02d}     {1:1.3f}'.format(i, intensity))
    except Exception as e:
        print(f'IR calculation not available: {e}')
else:
    print('DFPT vibration calculation not found')
    print('Run VASP with IBRION=7 and LEPSILON=.TRUE. for IR intensities')

mode  Relative intensity
   0.234
   0.045
   0.166
   0.616
   0.000
   0.198
   0.088
   0.374
   0.169