Calculating phonons with GoBaby
Hi there. How would you like to calculate phonons? Alright, then; here we goooooooooo!
The GoBaby script was written by Vidvuds Ozolins (who is currently at UCLA). I think that, in addition to calculating phonons, it is capable of lending some automation to running VASP jobs, but I'm not sure. It calculates phonons via diagonalization of the dynamical matrix. There are four scripts that make up the whole shebang: phon_init, phon_$systemName (where "$systemName" is the name of your system), phon_fin, and atomSx.
To begin, you need the following files: POTCAR, KPOINTS_phon, INCAR_phon, and apos.dat. POTCAR, KPOINTS, and INCAR are just the typical VASP input files with nothing special except a slight name change or two. apos.dat, on the other hand, is a special file. As of yet, it must be made by hand, but it sure seems like an easy thing to script. Here's a sample of what it might look like:
4.0034900755352907 0 0 #The three lattice vectors
0 4.0034900755352907 0
0 0 4.0034900755352907
1 1 1 # Number of unit cells lengths in each direction
2 1 # Number of species AND num. of q-points (=1, don't change)
6.941 1.008 # Atomic masses of species in a.m.u.
4 4 # Num of atoms of each species
0 0 0 # The cartesian version of the atomic positions
0 2.0017450377676453 2.0017450377676453
2.0017450377676453 0 2.0017450377676453
2.0017450377676453 2.0017450377676453 0
2.0017450377676453 2.0017450377676453 2.0017450377676453
2.0017450377676453 0 0
0 2.0017450377676453 0
0 0 2.0017450377676453
0.0 0.0 0.0 1.0 # Don't change these numbers
Here's a tip: You typically want the final supercell to be as cubic as possible, so finagle the unit cell lengths to make that happen. I haven't actually used GoBaby yet, so I don't know how seriously this consideration should be taken. Another tip: You can use DirectToCart to convert a POSCAR from direct to cartesian. DirectToCart was written by Kyle Michel, currently a graduate student in Professor Ozolins' lab at UCLA.
#Snippet from phon_init (referenced below)
...
23 # How many displacement points for a frozen phonon?
24 NSTEPS=5
25
26 # What is the step (in Angstroms) for a frozen phonon?
27 STEP=0.03
...
You may want to edit the variables STEP and NSTEPS in phon_init (shown above). STEP is the length of the perturbation. NSTEPS is the total number of positions that each perturbed atom will occupy in each symmetry-inequivalent direction. For instance, STEP=.03 and NSTEPS=5 will have an atom occupy its equilibrium position and positions .06, .03, -.03, and -.06 angstroms from equilibrium in each symmetry-inequivalent direction.
Running phon_init (via "./phon_init") will create all necessary files for phonon calculation, including POSCARs for each perturbation.
#Snippet from phon_$systemName (referenced below)
...
45 set qqversion =
46 set qqapp = "openmpi parallel"
47 set qqptasks = 4
48 set qqidir = /home/kylemichel/Danny/Phonons/sub$i
49 set qqjob = vasp
50 set qqodir = /home/kylemichel/Danny/Phonons/sub$i
51 cd /home/kylemichel/Danny/Phonons/sub$i
...
Now you should edit the directories in phon_$systemName ("qqidir", "qqodir", and the line under the "qqodir" line, shown above) to match the directory that holds your job. Then, run phon_$systemName, using "qsub phon_$systemName" because this is the most computationally intensive step. Each POSCAR will now have its own subdirectory. Run phon_fin now. This will do all the dynamical matrix stuff.
At this point, you have all the phonon information at your disposal. Check phonons.out. Look at all those numbers! That's the stuff that's gonna take your work from "theoretical" to "predictive". But we aren't done yet! As you can maybe see, a little bit more processing needs to be done.
#Snippet from a sample phonons.out, referenced below
1 8 1
2 0.00000000 0.00000000 0.00000000 1.0000000000
3 -1.28 -1.28 -1.28 282.57 282.57 282.57 #Look at the first three numbers on this line
4 282.57 282.57 282.57 381.62 381.62 381.62
5 400.37 400.37 400.37 701.46 701.46 701.46
...
So, rename phonons.out as "ev.out". We do this because the script we are about to use looks for ev.out. Open ev.out. Look at the first three numbers in the first row of six numbers (demonstrated above). Those are the standard deviations in the phonon frequencies, and they should be close to 0. To decide what "close" is, we can turn to math; the numbers should all be less than the reciprocal of the geometric mean of the masses. That is to say, where M_c is the mass of the cth species and the cell contains n species.
If you don't meet this criterion, something went wrong. If you do, though, you can continue. Change those three numbers to 0 before moving on. Run atomSx and redirect the output to a file (like "./atomSx > EnergyData").
This final file, EnergyData, is the goal. The numbers are arranged in four columns. The first is the temperature in Kelvin. The second is the entropy in kBoltzmanConstants/atom. The third is the heat capacity in kBoltzmanConstants/atom. And the fourth is the enthalpy in meV/atom.
Automation
apos.dat is definitely the most annoying file to create. Here's a script to do the annoying parts for you:
#!/bin/bash
gunzip OUTCAR.static.gz
cp CONTCAR apos.dat
cp CONTCAR.static apos.dat
sed -i '7,600d' apos.dat
sed -i '2d' apos.dat
sed -i '1d' apos.dat
grep PAW_GGA POTCAR | grep -v TITEL | awk '{print $2}' >> apos.dat
echo "#Number of supercell units, x,y,z" | cat >> apos.dat
echo "x 1 # of species in compound, keep 1 as it is. " | cat >> apos.dat
echo "#atomic weight" | cat >> apos.dat
echo "# like poscar, number of atoms in each" | cat >> apos.dat
cp OUTCAR OUTCAR.static
number=$(awk 'NR==1 { print $1}' DOSCAR)
#grep " position of ions in cartesian coordinates " OUTCAR.static -A $number
awk '/position of ions in cartesian coordinates/{c="'"$number"'";next}c{c--;print}' OUTCAR.static | cat >> apos.dat
echo "" | cat >> apos.dat
echo " 0.0 0.0 0.0 1.0 " | cat >> apos.dat
cp ../INCAR_phon .
cp ../kpgen.dat .
cp POTCAR POTCAR_phon
cp KPOINTS KPOINTS_phon
NAME=$(find *.q)
cp $NAME phon.$NAME
I call this script mapos. Just run mapos in the directory, and it should create a good apos.dat file.
Now, go into apos.dat, and replace the middle lines with your appropriate values. Make sure to format it so the syntax and format is correct!
Also, don't worry about the error messages, I've just programmed in redundancy so the file can be generalized.
Bonus!
You can use another of Kyle's utilities to quickly get the numbers you want at 300K. Just run MakeEnergies (which looks for a file called "EnergyData").