9610
Comment:

9618
converted to 1.6 markup

Deletions are marked like this.  Additions are marked like this. 
Line 39:  Line 39: 
source /gbar/hald/home3/s05/s053338/.10302.sh  source /xbar/nas2/home3/rnl/s053338/.10302.sh 
Line 44:  Line 44: 
source /gbar/hald/home3/s05/s053338/.10302.csh  source /xbar/nas2/home3/rnl/s053338/.10302.csh 
Line 84:  Line 84: 
.. _`Alfccsingle.py`: attachment:Alfccsingle.py  .. _`Alfccsingle.py`: [[attachment:Alfccsingle.py]] 
Line 205:  Line 205: 
.. _Alfcc.py: attachment:Alfcc.py  .. _Alfcc.py: [[attachment:Alfcc.py]] 
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 a_{0} = 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 a_{0} = 4.05 Å. For the first example, we choose 300 eV as plane wave cutoff (400 eV density cutoff) and 4 x 4 x 4 kpoints. Copy the script:
Alfccsingle.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 Alfccsingle.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:
Alfccsingle.nc Alfccsingle.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 Alfccsingle.txt either by opening it in emacs or by using:
more Alfccsingle.txt
You can conveniently monitor some variables by using the grep utility. By typing:
grep conv Alfccsingle.txt
you see the changes in density, occupation numbers, energy etc. during the SCF cycle. By typing
grep DFT Alfccsingle.txt
you see the total energy, as it evolves during the SCF cycles. The energy for the XCfunctional in use is marked as selfcons and the others as nonselfcons. As we have not indicated any XCfunctional in the script, Dacapo uses its default XCfunctional. Which functional is that?
By typing:
grep EIG Alfccsingle.txt
we get the eigenvalues for the bands at the different kpoints. 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('Alfccsingle.nc') >>> e = bulk.GetPotentialEnergy() >>> calc = bulk.GetCalculator() >>> density = calc.GetDensityArray() >>> from ASE.IO.Cube import WriteCube >>> WriteCube(bulk, density, 'Al.cube') >>> [hit CTRLd] $ vmd Al.cube
Try to make VMD show an isosurface of the electron density.
3.2 Convergence in kpoints and planewave energy cutoff
Now, we will investigate the necessary kpoint sampling and kinetic energy plane wave cutoff for bulk fcc Aluminum at the experimental lattice constant a_{0} = 4.05 Å; this is a standard first step in all plane wave calculations.
Copy the script Alfcc.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 Alfcc.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 kpoint sampling.
Do you expect that the kpoint / cutoff energy test is universal for all other Aluminum structures than fcc? What about other chemical elements ?
Are both kpoint / 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 kpoint 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 a_{0} = 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 a_{0} and the energy curve minimum E_{0} = E(a_{0}). From your fit, calculate the bulk modulus
B = M/(9a_{0})d^{2}E/da^{2}
for a = a_{0}, 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 builtin fitting tools
 data > Transformations > Regression
 data > Transformations > Non linear curve fitting
 data > Transformations > Splines
Compare your results to the experimental values a_{0} = 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:
b_{1} = a_{bcc} (0.5, 0.5, 0.5)
b_{2} = a_{bcc} (0.5, 0.5, 0.5)
b_{3} = a_{bcc} (0.5, 0.5, 0.5)
Make a qualified starting guess on a_{bcc} 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 a_{bcc} 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 a_{bcc}, using a few points. Is it a good idea to use the same kpoint 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 kpoints, make sure that your kpoint sampling is dense enough.