Editing Intensities.sh

go back to [[Main Page]], [[Computational Resources]], [[Scripts]], [[Scripts for VASP]] or [[VASP]]

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

== DESCRIPTION ==

Although vibrational analysis is a kind of routine task, sometimes it is needed to get fast the intensities. 

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: <TT>IBRION=5</TT>, <TT>NFREE=2</TT>, <TT>POTIM="small_number"</TT> (whereas <TT>"small_number" < 0.05</TT>). Sometimes it's okay to proceed the "infrared" analysis even at lower ''k''-point grid in spite of speed.

In older VASP versions (4.6 or lower), 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'' ([http://cms.mpi.univie.ac.at/vasp/Berry_phase/node4.html 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 [http://dx.doi.org/10.1063/1.466753 1] and [http://dx.doi.org/10.1103/PhysRevB.57.223 2], theory [http://dx.doi.org/10.1103/RevModPhys.73.515 3], further reading: [http://othes.univie.ac.at/10117/ David's PhD thesis], chapter 2.10.3).

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

<p align="left" style="width:100%">
{| border="1"
 |- 
 | <TT>IBRION = 7</TT>
 | switches on the DFPT vibrational analysis (with no symmetry constraints) ([http://cms.mpi.univie.ac.at/vasp/vasp/Optical_properties_density_functional_perturbation_theory_PT.html man])
 |- 
 | <TT>LEPSILON = .TRUE.</TT>
 | enables to calculate and prints BEC tensor ([http://cms.mpi.univie.ac.at/vasp/vasp/LEPSILON_static_dielectric_matrix_ion_clamped_piezoelectric_tensor_Born_effective_charges.html man])
 |- 
 | <TT>NSW = 1</TT>
 | reduces # of ionic steps to 1; unlike for <TT>IBRION=5</TT>, here the setting is very important
 |}
</p>

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

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

The script <TT>intensities.sh</TT> itself:

== CODE ==

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

--[[User:Dkarhanek|Dkarhanek]] 16:44, 30 September 2010 (CEST)