Atomic Simulation Environment

for lithium version of Exciting

Purpose: In this tutorial you will learn how to use the simulation tool ASE to setup exciting calculations.

0. Define relevant shell variables and download scripts

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant shell variables are already defined and that the excitingscripts directory has already been downloaded, as specified in Tutorial scripts and environment variables.

From now on the symbol $ will indicate the shell prompt.

1. About ASE

Developed at CAMd, the Atomic Simulation Environment (ASE) provides Python modules for manipulating atoms, analyzing simulations, visualizing results, etc. In ASE, a calculator is a black box that can take input structure information from an input object (e.g., the atoms object) and return properties like the total energy, forces, and stresses.

ASE scripts are written in the general language Python.

2. How to run a Python script

Python is a general-purpose, interpreted high-level programming language whose design philosophy emphasizes code readability. The use of the Python language allows ASE to be used both interactively as well as in scripts. You can start Python interactively by typing in the command line

$ python
Python 2.6.6 (r266:84292, Dec 26 2010, 22:31:48) 
[GCC 4.4.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>

The Python prompt (>>>) appears on the screen. You can now insert your Python command lines and use Ctrl-D to quit Python.

Alternatively, you can collect your Python commands in a script (e.g., named my-python-script.py) and execute it as

$ python my-python-script.py

If the following line

#!/usr/bin/env python

is included as the first line of your script, you can execute directly the script using

$ chmod +x my-python-script.py
$ ./my-python-script.py

3. Features in ASE

Go to the ASE root directory ($EXCITINGASE) and have there a look at the submodules and functionalities. There are many scripts and directories, which contribute to a number of "design goals". Most importantly, one can easily access the codes and make any kind of modification, which represents an advantage in terms of flexibility of the code. In this section, we will demonstrate how to generate input structures, apply constraints, perform force and stress optimization, and visualize the result.

Please note: We are introducing ASE to advanced users, for more details and for a complete information please go to the ASE website. Further useful information can be obtained by watching the ASE mailing list

4. Generating a crystal structure with ASE

There are several ways to build the initial crystal structure

  1. Specify the atomic species, position, and cell parameters manually.
  2. Read an existing structure file of known format, check the available file formats here.
  3. Use the lattice, structure, spacegroup, or surface modules to create an input structure from the avaliable crystal database.
4.1 Build the structure from scratch

Here, we will consider as example the wurtzite structure of bulk GaN. The first step is to create a working directory.

$ mkdir ase-work
$ cd ase-work

Next, save the following code in a script file named GaN.py.

from ase import *
from math import *
 
a = 3.189
c = 2*sqrt(2/3.)*a
atoms = Atoms([Atom('Ga', (0,       0,    0)),
               Atom('Ga', (1/3., 1/3., 1/2.)),
               Atom('N',  (0,       0, 3/8.)),
               Atom('N',  (1/3., 1/3., 7/8.))])
 
cell = [(a,             0, 0),
        (a/2, a*sqrt(3)/2, 0),
        (0,             0, c)]
atoms.set_cell(cell, scale_atoms=True)
 
from ase.visualize import view
view(atoms)

As you can see in this script, first you have to import the default modules in ASE and then create the atoms object. The last two lines enable you to visualize the structure in a prompt out window. To run the python script, type the following command.

$ python GaN.py

Hold the right buttom of the mouse and drag to rotate the structure in the pushed window, the left buttom can be used to select atoms. Information about the object you have selected appears at the bottom of the window.

4.2 Write and read structure

In order to write your structure into a file and read from an existing structure file, you have to use the ase.io module. This can be done by adding the following lines to the script GaN.py.

from ase.io import *
write('GaN.exi', atoms, format='exi')
struct = read('GaN.exi', format='exi')
view(struct)

After the execution of the script, a new structure file GaN.exi is saved in your working directory.

4.3 Structure database in ASE

A list of examples for structures that you can generate with ASE is given in the following.

#------------------------------------------------------------------
# Al bulk packed in the (111) direction (redefine lattice vectors)
from ase.lattice.cubic import FaceCenteredCubic
Al111 = FaceCenteredCubic(directions=[[1,-1,0], [1,1,-2], [1,1,1]],
            size=(2,2,3), 
            symbol='Al', 
            pbc=(1,1,1), 
            latticeconstant=4.05)
 
#------------------------------------------------------------------
# create a compound structure
from ase.lattice.compounds import Zincblende
GaN = Zincblende(directions=[[1,0,0], [0,1,0], [0,0,1]], 
            size=(1,1,1), 
            symbol=('Ga','N'), 
            latticeconstant=4.50)
 
#------------------------------------------------------------------
# complicate skutterudites structure with 32 atoms
from ase.lattice.spacegroup import crystal
skutterudite = crystal(('Co', 'Sb'), 
            basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)], 
            spacegroup=204, 
            cellpar=[9.04, 9.04, 9.04, 90, 90, 90])
 
#------------------------------------------------------------------
# create surface with vacuum
from ase.lattice.surface import *
slab = fcc111('Al', size=(2,2,3), a=4.05, vacuum=10.0)
 
#------------------------------------------------------------------
# add an adsorbte on top
from ase.lattice.surface import *
slab = fcc111('Al', size=(2,2,3))
add_adsorbate(slab, 'H', 1.5, 'ontop')
slab.center(vacuum=10.0, axis=2)
 
#------------------------------------------------------------------
# nanotube
from ase.structure import nanotube
cnt = nanotube(3, 3, length=6, bond=1.4, symbol='Si')
 
#------------------------------------------------------------------
# nanoribbons
from ase.structure import graphene_nanoribbon
gnr = graphene_nanoribbon(2, 6, type='zigzag', 
            saturated=True, C_H=1.1, C_C=1.4, 
            vacuum=6.0, 
            magnetic=True, 
            initial_mag=1.12, 
            sheet=False, 
            main_element='C', 
            saturate_element='H', 
            vacc=None)

Save it into a script structure-list.py and start Python interactively by typing

$ python -i structure-list.py

You can check the structures using the view module in ASE. For example, the first structure, corresponding to bulk fcc Al packed in the (111) direction, can be visualized by typing

>>> from ase.visualize import view
>>> view(Al111)

5. Other features

5.1 Constraints

The following types of constraint are avaliable in ASE.

FixCartesian         FixBondLength     FixedMode
FixConstraintSingle  FixAtoms          UnitCellFilter
FixScaled            StrainFilter      FixedPlane 
Filter               FixConstraint     FixedLine 
FixBondLengths       FixInternals

Try to explore the effect of these constraints by having a look to the script $EXCITINGASE/ase/constraints.py. Here are examples of usage for some specific constraints.

#----------------------------------------------------------
# fix positions of all Cu atoms in your structure
fa = FixAtoms(mask=[atom.symbol == 'Cu' for atom in atoms])
atoms.set_constraint(fa)
 
#----------------------------------------------------------
# fix the bond length between the 1st and 2nd atom, 
# as well as the bond legth between the 1st and 3rd atom
fb = FixBondLengths([[0, 1], [0, 2]])
atoms.set_constraint(fb)
 
#----------------------------------------------------------
# combine constraints
atoms.set_constraint([fa, fb])
5.2 Force optimization

In order to understand how optimization is performed in ASE, consider the following example. Simply run this simple script for the water molecule and try to understand what it is doing.

from ase import *
from math import *
import numpy as np
from ase.optimize import *
from ase.calculators.emt import EMT
 
d = 0.9575
t = pi / 180 * 104.51
water = Atoms('H2O',
              positions=[(d, 0, 0),
                         (d * np.cos(t), d * np.sin(t), 0),
                         (0, 0, 0)],
              calculator=EMT())
 
dyn = BFGS(water, trajectory='H2O.traj')
dyn.run(fmax=0.05)

Indeed, the script performs the optimization of the forces acting on the atoms of a water molecule using the simple EMT calculator. Save the script as water.py and run it.

$ python water.py

Optimization trajectories can be visualized using the basic visualization tool ag which is included in ASE. In order to do this, use the command

$ ag H2O.traj

More details on the visualization tolls in ASE can be found in the ASE documentation.

The general Python block for performing the force optimization is the following.

from ase.optimize import *
 
dyn = BFGS(atoms=system, trajectory='qn.traj', logfile='qn.log', restart='qn.pckl')
dyn.run(fmax=0.05)        ### calculation will start from this line

For further details on the keywords used in the ase.optimize Python library, please refer to the ASE.optimize documentation.

6. The exciting calculator

There are actually two ways to construct the exciting calculator.

  1. Using keyword arguments to specify groundstate attributes.
  2. Use paramdict to specify an input file structure with a data structure of dictionaries.
6.1 Using groundstate keywords

Some of the keywords/attributes of the exciting groundstate element can be specified directly:

from ase.calculators.exciting import Exciting
 
calc = Exciting(bin='excitingser', kpts=(4, 4, 3), xctype='GGArevPBE')
6.2 Using the parameter dictionary

The calculator can translate the parameters specified inside the paramdict dictionary to the exciting input format. Below is an example which uses this procedure.

from ase.structure import bulk
from ase.calculators.exciting import Exciting
 
atoms = bulk('Si', 'diamond', a=5.459)
 
calc = Exciting(dir='calc',
           paramdict={"title":{"text()":"Silicon"},
                      "groundstate":{"ngridk":"4 4 4","tforce":"true"},
                      "structureoptimization":{},
                      "properties":{"dos":{},
                      "bandstructure":{"plot1d":{
                      "path":{ "steps":"100",
                      "point":[{"coord":"0.75000   0.50000   0.25000", "label":"W"},
                               {"coord":"0.50000   0.50000   0.50000", "label":"L"},
                               {"coord":"0.00000   0.00000   0.00000", "label":"GAMMA"},
                               {"coord":"0.50000   0.50000   0.00000", "label":"X"},
                               {"coord":"0.75000   0.50000   0.25000", "label":"W"},
                               {"coord":"0.75000   0.37500   0.37500", "label":"K"}]}}}}})
 
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
print 'Energy =', '%15.8f' %energy, ' eV'

This script can be used for the calculation of the total density of states (DOS) and bandstructure of bulk silicon. Save the script as bulk-Si.py and run it.

$ python bulk-Si.py

After running the calculation, you can use exciting templates for visualizing the DOS and the bandstructure.

$ xsltproc --stringparam ID "t///" $EXCITINGVISUAL/xmldos2grace.xsl dos.xml > Silicon_dos.agr
$ xsltproc --stringparam unit "Ha" $EXCITINGVISUAL/xmlband2agr.xsl bandstructure.xml

The total DOS (file Silicon_dos.agr) and bandstructure (Silicon_bandstructure.agr) are generated in the xmgrace format. The result should be similar to the following.

dos.png band.png

7. Additional notes on running exciting with ASE

7.1 When does an exciting calculation starts?

In your script, if you have defined the atoms object and an exciting calculator, the calculation starts immediately when you use the get_total_energy() or get_potential_energy() method from the atoms object, e.g.

from ase.structure import bulk
from ase.calculators.exciting import Exciting
 
atoms = bulk('Si', 'diamond', a=5.459)
 
calc = Exciting(dir='calc', bin='excitingser', kpts=(4, 4, 3))
 
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
print 'Total energy =', energy, 'in eV'
7.2 How to submit a job to a queueing system

Typing some additional lines in your script, you can submit it to the queueing system of a computer cluster.

#!/bin/bash -x
# PBS -l nodes=1:ppn=1
 
cd $PBS_O_WORKDIR
 
### start of job
 
python <<!
from ase.structure import bulk
from ase.calculators.exciting import Exciting
 
atoms = bulk('Si', 'diamond', a=5.459)
 
calc = Exciting(dir='calc', bin='excitingser', kpts=(4, 4, 3))
 
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
print 'Total energy =', energy, 'in eV'

6. Exercises

  • Construct with ASE a supercell containing more than 100 atoms. Then, write the structure into a file with your favorite format.
  • Try some examples from the ASE tutorial page.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License