# 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 **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 *a*_{0} = 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*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*/(9*a*_{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 built-in 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**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.