--------------------------- Getting started with Dacapo --------------------------- .. contents:: .. section-numbering:: .. |angst| unicode:: U+0212B .. ANGSTROM SIGN .. |infin| unicode:: U+0221E .. INFINITY .. |simeq| unicode:: U+02243 .. ASYMPTOTICALLY EQUAL TO .. |sigma| unicode:: U+003C3 .. GREEK SMALL LETTER SIGMA .. |Delta| unicode:: U+00394 .. GREEK CAPITAL LETTER DELTA .. |mu| unicode:: U+003BC .. GREEK SMALL LETTER MU .. |beta| unicode:: U+003B2 .. GREEK SMALL LETTER BETA .. |pi| unicode:: U+003C0 .. GREEK SMALL LETTER PI .. |alpha| unicode:: U+003B1 .. GREEK SMALL LETTER ALPHA .. |nu| unicode:: U+003BD .. GREEK SMALL LETTER NU .. |deg| unicode:: U+000B0 .. DEGREE SIGN Atomic Simulation Environment ============================= It is a good idea to quickly skim through the first 5 sections of the `ASE manual `_ 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``. Bulk Al ======= .. experimental a0 = 4.05 A experimental bulk modulus = 76 GPa = (76*10.**9)/(1.602*10.**(-19))/10.**30 = 0.474406 V/A^3 experimental Ec = 3.34 eV/atom .. |a0| replace:: *a*\ :sub:`0` 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 Å. 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 |angst|. 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`_ .. _`Al-fcc-single.py`: attachment: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. 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? .. We use Ecut = 250 ev and MP = 8,8,8 for fcc and MP = 10,10,10 for bcc .. _Al-fcc.py: attachment:Al-fcc.py 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*\ (|infin|) 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 E\ :sub:`0` = E(|a0|). From your fit, calculate the bulk modulus *B* = *M*/(9\ |a0|\ )\ *d*\ :sup:`2`\ *E*/*da*\ :sup:`2` 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? 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**\ :sub:`1` = *a*\ :sub:`bcc` (-0.5, 0.5, 0.5) **b**\ :sub:`2` = *a*\ :sub:`bcc` (0.5, -0.5, 0.5) **b**\ :sub:`3` = *a*\ :sub:`bcc` (0.5, 0.5, -0.5) * Make a qualified starting guess on *a*\ :sub:`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*\ :sub:`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*\ :sub:`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.