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:
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.
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.
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.
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.
apos.dat is definitely the most annoying file to create. Here's a script to do the annoying parts for you:
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.
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").