Intensities.sh

go back to Main Page, Computational Resources, Scripts, Scripts for VASP

An article about VASP vibrational intensities using DFPT (Density-Functional Perturbation Theory), a.k.a. LRT (Linear Response Theory).

Although vibrational analysis is a kind of routine task, sometimes it needs a bit of playing around in order to be done really properly, or to get some extra information (yes, the intensities, vide infra).

Mostly we just want to "check" if our optimization finished to a reasonable minimum or that the transition state we found is a pretty saddle point we wanted to have... well, in these cases, we can be satisfied with changing the INCAR settings to be: IBRION=5, NFREE=2, POTIM="small_number" (whereas "small_number" < 0.05). Sometimes it's okay to proceed the "infrared" analysis even at lower k-point grid in spite of speed.

In older VASP versions (< 5.*), the intensities of vibrational modes could only be calculated (with obstructions) based on the change of dipole moment for each vibration direction separately. This is the same implementation as it is being routinely used while implemented in software suites like GAUSSIAN - where the corresponding source code comes from 1970's.

As of VASP 5.* version, the DFPT linear response calculations are available. Among others, the user can obtain the matrix of Born effective charges (BEC), which refers to change of atoms' polarizabilities w.r.t. an external electric field. The BEC tensor is a key to calculate the vibrational intensities using the most modern method available (as far as I know!), using the formula by Gianozzi&Baroni (aplications 1 and 2, theory 3, further reading: David's PhD thesis, chapter 2.10.3).

For calculations with VASP DPFT formalism, include the following in the INCAR file:

IBRION = 7	switches on the DFPT vibrational analysis (with no symmetry constraints) (man)
LEPSILON = .TRUE.	enables to calculate and prints BEC tensor (man)
NSW = 1	reduces # of ionic steps to 1; unlike for IBRION=5, here the setting is very important

Note In the INCAR file, please take care that you also obtain the electronic density properly converged with NELMIN=10 in combination with NELMIN=120 and always SWITCH OFF THE SYMMETRY with ISYM=0 tag.

In case that you want to extract intensities from OUTCAR, submit the below script without any parameters (it reads the OUTCAR in the current directory). The output is written into "intensities" subdirectory, the most imporant output is to be found in the file "intensities/results/results.txt", whereas "intensities/results/exact.res.txt contains intensities before normalization".

The script intensities.sh itself:

---

#!/bin/bash
# A utility for calculating the vibrational intensities from VASP output (OUTCAR)
# (C) David Karhanek, 2009-08-05

# extract Born effective charges tensors
printf "..reading OUTCAR"
BORN_NROWS=`grep NIONS OUTCAR | awk '{print $12*4+1}'`
grep "in e, cummulative" -A $BORN_NROWS OUTCAR > born.txt

# extract Eigenvectors and eigenvalues
EIG_NVIBS=`grep cm-1 OUTCAR | wc -l`
EIG_NIONS=`grep NIONS OUTCAR | awk '{print $12}'`
EIG_NROWS=`echo "($EIG_NIONS+3)*$EIG_NVIBS+3" | bc`
grep "Eigenvectors and eigenvalues" -A $EIG_NROWS OUTCAR | 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)" | bc -l`
 done
 freq=`awk '(NR==1){print $8}' temps02/eigens.vib.$nuu`
 printf "%3u %11.6f %16.14f\n" $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 "%3u %11.6f %16.14f\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*

# list results
cat results/results.txt

---

--Dkarhanek 13:28, 29 September 2010 (CEST)