Getting started with Dacapo
Contents
1 Atomic Simulation Environment
It is a good idea to quickly skim through the first 5 sections of the ASE manual
2 Setting up your UNIX environment
First, find out what type of shell you are using. Type:
echo $SHELL
If you are running a bash shell, include this line in your .bashrc file:
source /xbar/nas2/home3/rnl/s053338/.10302.sh
If you are running a tcsh or csh shell include this line in your .tcshrc or .cshrc file:
source /xbar/nas2/home3/rnl/s053338/.10302.csh
The next time you log in, everything will be ready to go! For this session, you will have to do:
source .bashrc
or source .cshrc or source .tcshrc.
3 Bulk Al
We will initially look at bulk Aluminum as a model system. The experimental equilibrium structure for Aluminum is fcc with cubic lattice constant a0 = 4.05 Å.
3.1 Running a calculation
Now we are ready to run the first dacapo calculation. We will look at bulk fcc aluminum and make a single energy calculaion for the experimental lattice constant a0 = 4.05 Å. For the first example, we choose 300 eV as plane wave cutoff (400 eV density cutoff) and 4 x 4 x 4 k-points. Copy the script:
Al-fcc-single.py
to a place in your file area. Read the script and try to get an idea of what it will do. Run the script by typing:
python -i Al-fcc-single.py
The option -i stands for interactive and means that the program has to be finished with CTRL+D manually. At the end of the program, a vmd window with the bulk structure pops up. Verify that the structure indeed is fcc.
Notice that the program has generated two output files:
Al-fcc-single.nc Al-fcc-single.txt
In general, when you execute a Dacapo electronic structure calculation, you get two files:
A binary data file (conventional suffix .nc) containing input/output data. You may convert this file to readable format by the command:
ncdump filename.nc
or:
ncdump filename.nc | more
if it is a big file.
An ASCII formatted log file (conventional suffix .txt) that monitors the progress of the calculation.
Try to take a look at the file Al-fcc-single.txt either by opening it in emacs or by using:
more Al-fcc-single.txt
You can conveniently monitor some variables by using the grep utility. By typing:
grep conv Al-fcc-single.txt
you see the changes in density, occupation numbers, energy etc. during the SCF cycle. By typing
grep DFT Al-fcc-single.txt
you see the total energy, as it evolves during the SCF cycles. The energy for the XC-functional in use is marked as selfcons and the others as non-selfcons. As we have not indicated any XC-functional in the script, Dacapo uses its default XC-functional. Which functional is that?
By typing:
grep EIG Al-fcc-single.txt
we get the eigenvalues for the bands at the different k-points. How many bands are occupied? What would you expect for Aluminum?
One can also monitor other variables, and these will be introduced in the coming exercises, as we need them.
The binary file contains all information about the calculation. Try typing the following from the Python interpreter:
>>> from Dacapo import Dacapo
>>> bulk = Dacapo.ReadAtoms('Al-fcc-single.nc')
>>> e = bulk.GetPotentialEnergy()
>>> calc = bulk.GetCalculator()
>>> density = calc.GetDensityArray()
>>> from ASE.IO.Cube import WriteCube
>>> WriteCube(bulk, density, 'Al.cube')
>>> [hit CTRL-d]
$ vmd Al.cube
Try to make VMD show an isosurface of the electron density.
3.2 Convergence in k-points and planewave energy cutoff
Now, we will investigate the necessary k-point sampling and kinetic energy plane wave cutoff for bulk fcc Aluminum at the experimental lattice constant a0 = 4.05 Å; this is a standard first step in all plane wave calculations.
Copy the script Al-fcc.py to a place in your file area. Read the script and get an idea of what it will do. Then run the script by typing:
python -i Al-fcc.py
A view of a piece of the crystal will pop up. Verify that the structure is indeed fcc by rotating the structure with the mouse.
Estimate the necessary values of the kinetic energy plane wave cutoff and k-point sampling.
Do you expect that the k-point / cutoff energy test is universal for all other Aluminum structures than fcc? What about other chemical elements ?
Are both k-point / cutoff energy convergence covered by the variational principle, i.e. are all your calculated energies upper bounds to true total energy?
Why do you think Al was chosen for this exercise?
4 Equilibrium lattice properties
Having determined the necessary values of the plane wave cutoff and k-point sampling, we now proceed to calculate some equilibrium lattice properties of bulk Aluminum.
First map out the cohesive curve E(a) for Al(fcc), i.e. the total energy as function of lattice constant a, around the experimental equilibrium value of a0 = 4.05 Å. Notice that the vacuum energy level E(∞) is not zero. Get four or more energy points, so that you can make a fit.
Fit the data you have obtained to get a0 and the energy curve minimum E0 = E(a0). From your fit, calculate the bulk modulus
B = M/(9a0)d2E/da2
for a = a0, where M is the number of atoms per cubic unit cell. Make the fit using your favorite math package (Mathematica/MatLab/Maple/Python/...) or xmgrace. xmgrace may take a data file with (x, y) in columns. Start xmgrace like:
xmgrace data_file_name
In xmgrace you may use one of the built-in fitting tools
- data -> Transformations -> Regression
- data -> Transformations -> Non linear curve fitting
- data -> Transformations -> Splines
Compare your results to the experimental values a0 = 4.05 Å and B = 76 GPa. Mind the units, when you calculate the bulk modulus. What are the possible error sources, and what quantity is more sensitive, the lattice constant or the bulk modulus?
5 Equilibrium lattice properties for bcc
Setup a similar calculation for bcc, in the minimal unit cell; a symmetric set of Bravais lattice vectors is:
b1 = abcc (-0.5, 0.5, 0.5)
b2 = abcc (0.5, -0.5, 0.5)
b3 = abcc (0.5, 0.5, -0.5)
Make a qualified starting guess on abcc from the lattice constant for fcc, that you have determined above. One can either assume that the primitive unit cell volumes of the fcc and bcc structure are the same or that the nearest neighbor distances are the same. Find a guess for abcc for both assumptions. Later, you can comment on which assumption gives the guess closer to the right lattice constant.
- Check that your structure is right by repeating the unit cell. You
can do this by plotting the structure in VMD like in the Al_fcc.py script.
Map out the cohesive curve E(a) for Al(bcc) and determine abcc, using a few points. Is it a good idea to use the same k-point setup parameters as for the fcc calculations? Calculate the bulk modulus, as it was done for fcc, and compare the result to the fcc bulk modulus. What would you expect?
Using the lattice constants determined above for fcc and bcc, calculate the fcc/bcc total energies at different plane wave cutoffs: 200 eV and 450 eV, i.e. four calculations. Compare the structure energy differences for the two cutoffs. Generally, in plane wave calculations, energy differences converge much faster with plane wave cutoff than total energies themselves. Further notice that the energy zero in pseudopotential calculations does not have physical significance. This exercise is sensitive to the number of k-points, make sure that your k-point sampling is dense enough.