... let DFT be Daily Fun Training

Downloads / Utilities / Scripts - David Karhánek

Some useful content appears in this section from time to time.

Outline (quick jump):

  Entry #01:   Revert XDATCAR from 5.2.11+ VASP version to 5.2.2-
  Entry #02:   Calculation of vibrational (IR) intensities in VASP 5.*

Entry #01:   Revert XDATCAR from 5.2.11+ VASP version to 5.2.2-   (2011-03-01)

I've noticed a significant change in header of XDATCAR file within VASP 5.2.* subversions and written a script for reverting this back:

sed '2,5d;8i

Note that the newline character is important. Thanks to the new header, XDATCARs have become independent of POSCARs/CONTCARs, though it serves problems to some older scripts or programs (e.g. visualization in STRender).

The script with an if condition (for automatization) might look like this (save it):

# Reverting XDATCAR
len=$(awk '(NR==8){print length}' XDATCAR)
if [ $len -lt 5 ]
 sed '2,5d;8i\
 echo "  Changed XDATCAR written to YDATCAR."
 echo "  No changes on XDATCAR needed."

Entry #02:   Calculation of vibrational (IR) intensities in VASP 5.*   (2011-04-06)

This contribution was inspired by the CCL mailing list discussion on "IR intensities from VASP" initiated and subsequently summarized by Dr. Ralf Tonner. Many thanks to Ralf for having tackled this question concerning an approach that should be a free added value for the researchers community using VASP.

The following text handles with 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, one can be satisfied with changing the INCAR settings to be: IBRION=5, NFREE=2, POTIM="small_number" (whereas i.e. < 0.05; reasonable values are e.g. 0.02 for a gas-phase molecule and 0.05 for adsorbed molecule). 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 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 - testing in VASP [4] and [5]).

    [1] P. Giannozzi and S. Baroni, J. Chem. Phys. 100, 8537 (1994).   [DOI: 10.1063/1.466753]
    [2] K. Esfarjani, Y. Hashi, J. Onoe, K. Takeuchi, Y. Kawazoe, Phys. Rev. B 57, 223 (1998).   [DOI: 10.1103/PhysRevB.57.223]
    [3] S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).   [DOI: 10.1103/RevModPhys.73.515]
    [4] D. Karhánek, PhD Thesis, (2010). See chapters 2.10.3 (theory) and 5.3.3, 5.5.3 and 5.6.3 (applications).   Available online.
    [5] D. Karhánek, T. Bučko, J. Hafner, J. Phys.: Condens. Matter 22 265006 (2010).   [DOI: 10.1088/0953-8984/22/26/265006]

For calculation of IR intensities within DPFT formalism in VASP 5.*, include the following tags in the INCAR file:

 IBRION = 7  # switches on the DFPT vibrational analysis (with no symmetry constraints) (man)
 LEPSILON = .TRUE.  # enables to calculate and to print the BEC tensors (man)
 NSW = 1  # reduces nr. of ionic steps to 1; unlike for IBRION=5, here the setting is very important
 NWRITE = 3  # affects OUTCAR verbosity: explicitly forces SQRT(mass)-divided eigenvectors to be printed

Unlike classical finite differences calculations, the DFPT approach seems to be absolutely insensitive to the setup of NFREE and POTIM tags in the INCAR file.

(1) In the INCAR file, please take care that you also obtain the electronic density properly converged with NELMIN=10 in combination with NELM=120 and, if possible, always SWITCH OFF THE SYMMETRY with ISYM=0 tag.
(2) Accuracy of the IR intensities goes pretty wrong when using real-space projection operators. → Ensure that you are using LREAL = .FALSE. (reciprocal space projection) setting. (Thanks to Dr. Ralf Tonner to point this out.)
(3) Linear response calculation is for some reason incompatible with the Γ (gamma) point only VASP binaries. Please avoid using them.

The BASH script for evaluation of IR intensities within DFPT by VASP (save it):

# 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/
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 
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 
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
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
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`
 for alpha in 1 2 3 ;  do            # summing over alpha coordinates
  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}
  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`
 freq=`awk '(NR==1){print $8}' temps02/eigens.vib.$nuu`
 echo "$nuu $freq $int">> exact.res.txt
 printf "."
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/
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

In order to extract intensities, submit the above script without any parameters (it reads the "OUTCAR" file in the current directory):
The output is written into a freshly created "intensities" subdirectory, the most imporant output is to be found in the file "intensities/results/results.txt", whereas the file "intensities/results/exact.res.txt" contains intensities before normalization. All vibrational modes are listed.

Additionaly, if rotational-vibrational modes should NOT be taken into account, use the following script, which extracts the normal modes (3N-6) only and and renormalizes their relative intensities again (save it):

# 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 ../..

The normal modes are printed to screen and saved into the file "intensities/results/reform.res.txt".

Finally, enjoy your results!!!

For questions, remarks, hints, suggestions, bug reports, ... → contact me via e-mail david.karhanek(at)gmail.com.

I'd like to thank Dr. Doris Vogtenhuber, Dr. Martijn Marsman and Dr. Tomás Bučko from VASP developers group for their valuable help on the topic.

Visitor Counter

Veni, vidi, vici.