26994
Comment:

27033

Deletions are marked like this.  Additions are marked like this. 
Line 205:  Line 205: 
`````````````````````````````````  `````````````````````````````````` 
Line 216:  Line 216: 
calc.SetExecutable('dacapo_2.7.8.run')  
Line 333:  Line 334: 
see the `dacapo pseudopotential page`_.  see the `Pseudopotential library`_ page. 
Line 338:  Line 339: 
``GetPseudoPotential(Z=<elementnumber>)``:  ``GetPseudoPotential(<elementnumber>)``: 
Line 341:  Line 342: 
``SetPseudoPotential(Z=<elementnumber>,path=path)``:  ``SetPseudoPotential(<elementnumber>, path)``: 
Line 569:  Line 570: 
calc.GetFermilevel()  calc.GetFermiLevel() 
Contents
 1 Introduction
 2 Setting up the atomic configuration
 3 Dacapo ASE calculator interface
 3.1 Number of bands and Electronic temperature
 3.2 Planewave and Densitycutoff
 3.3 Changing file names
 3.4 Exchangecorrelation functional
 3.5 Defining a spinpolarized calculation
 3.6 Kpoints
 3.7 Pseudopotentials
 3.8 Chargedensity mixing
 3.9 Eigenvalue solver
 3.10 Dipolecorrection
 3.11 Setting nondefault executable
 3.12 Symmetry
 3.13 Electrostatic decoupling
 3.14 Stress calculation
 3.15 Local density of states
 3.16 Changing convergence parameters
 3.17 dacapo fortran program/python interface communication
 3.18 NetCDF specific methods
 4 Getting information from the calculator
 5 Working with the Dacapo calculator
 6 Minimization and Dynamics using python
 7 Wavefunction and density plots
 8 Wannier orbitals
 9 Nudged Elastic Band
1 Introduction
Dacapo is a total energy program based on density functional theory. It uses a plane wave basis for the valence electronic states and describes the coreelectron interactions with Vanderbilt ultrasoft pseudopotentials. The program performs selfconsistent calculations for both Local Density Approximation (LDA) and various Generalized Gradient Approximation (GGA) exchangecorrelations potentials, using stateofart iterative algorithms. The code may perform molecular dynamics / structural relaxation simultaneous with solving the Schrodinger equations within density functional theory. The program may be compiled for seriel as well as parallel execution and the code has been ported to many hardware platforms.
The manual describes how to use the dacapo program from the Campos Atomistic Simulation enviroment, or for short ASE.
2 Setting up the atomic configuration
The atoms for Dacapo are defined via the ASE ListOfAtoms object.
A ListOfAtoms object is a collection of atoms. The ASE implementation of a ListOfAtoms can be used like this:
from ASE import Atom, ListOfAtoms atoms = ListOfAtoms([Atom('C', (0, 0, 0)), ... Atom('O', (0, 0, 1.0))]) atoms.SetUnitCell((10,10,10))
This will define a CO molecule, with a CO distance of 1A, in a cubic 10Ax10Ax10A unitcell.
For more details of the ListOfAtoms class see listofatoms in the ASE manual.
The atoms defined here can be attached to a Dacapo calculator by using:
atoms.SetCalculator(calculator)
alternatively, one can attach the atoms to a calculator by using:
calculator.SetListOfAtoms(atoms)
The next chapter describes how to define the calculator for dacapo.
3 Dacapo ASE calculator interface
As a minimum for defining a calculation in python you need to specify the atoms as described above, and you need to set the planewave cutoff and the number of bands, this will look like this using the python interface:
from Dacapo import Dacapo calculator = Dacapo() calculator.SetNumberOfBands(6) calculator.SetPlaneWaveCutoff(340) atoms.SetCalculator(calculator) calculator.GetPotentialEnergy()
The rest of this chapter describes the different methods you can use for the calculator object. More detailed explanation of the python construction found above can be found in the ASE Tutorials.
For the most important of the method described one can use keyword arguments in the construction of the Dacapo calculator, the table below list the keywords:
Method  keyword 

SetPlaneWaveCutoff  planewavecutoff 
SetDensityCutoff  densitycutoff 
SetNumberOfBands  nbands 
SetXCFunctional  xc 
SetBZKPoints  kpts 
SetSpinPolarized  spin 
SetNetCDFFile  out 
SetTxtFile  txtout 
DipoleCorrection  dipole 
SetScratch  scratchdir 
The example above can now be written more compact:
calculator = Dacapo( ... nbands=6, ... planewavecutoff=340)
3.1 Number of bands and Electronic temperature
The method 'SetNumberOfBands(nbands)' can be used to define the number of bands for the calculatior. nbands should be a least equal to halv the number of valence electrons.
The occupation statitstics can be set using the method SetOccupationStatistics(method).
dacapo support two methods for occupation statistics:
 FermiDirac:
 FermiDirac occupation statistics. The temperature is controlled by the method SetElectronicTemperature(temp). Default is 0.1 eV.
 MethfesselPaxton:
 Currently, the MethfesselPaxton scheme (PRB 40,3616, 1989) is implemented to 1st order.
3.2 Planewave and Densitycutoff
The method SetPlaneWaveCutoff(planewavecutoff) defines the plavewave cutoff for the calculation. The unit for planewavecutoff is eV.
The method SetDensityCutoff(densitycutoff) sets the density cutoff. By default the two cutoffs are equal, corresponding to using just one grid. For elements like Cu using ultrasoft pseudo potentials the ultrasoft pseudowavefunction are very soft, while the density still requires a fine grid, in this case it can be advantagous to set a larger value for densitycutoff that for planewavecutoff.
3.3 Changing file names
The following two methods can be use to change the file names used by dacapo:
 SetNetCDFFile(filename):
 This will set the NetCDF filename, default is 'out.nc'.
 SetTxtFile(filename)
 This will set the Ascii text file name, default is 'out.txt'. For information on the most inportant data in the Ascii file see dacapo Ascii file.
The method:
SetScratch(scratchdir):
will set the scratch directory used for swapfiles to scratchdir. Default for this directory is /scratch/<username> or if this is not found just /scratch.
3.3.1 Wavefunction of scratch directory
For some applications like NEB calculations, running from a NFS mounted file system, NFS traffic can be reduced by also placing the wavefunctions on the local disk. For version 2.7.8 this can also be done for parallel calculations, because only the master nodes do any reading/writing of the NetCDF input/output files. The following will set do this:
calc.SetExecutable('dacapo_2.7.8.run') calc.infile = '/scratch/<username>/in.nc' calc.SetNetCDFFile('/scratch/<username>/out.nc') calc.SetTxtFile('out.txt')
The txt file must be explicit set, otherwise it will also go to the /scratch directory. In general the wavefunction will be lost doing this, and the calculation must be restarted, not from the dacapo nc file, but from eg. NetCDFTrajectory files.
3.4 Exchangecorrelation functional
You can use the method SetXCFunctional(exc) to select the exchangecorrelation used for the calculation. exc is one of:
 LDA:
 Vosko Wilk Nusair LDAparametrization.
 PW91:
 Perdew Wang 91 GGAparametrization.
 PBE:
 Perdew Burke Ernzerhof GGAparametrization.
 RPBE:
 Revised PBE GGAparametrization, Hammer, Hansen and Nørskov.
The method:
GetXCFunctional
tells which of the above defined exchangecorrelation functionals are currently used.
The method:
GetXCEnergies(xcname)
return the nonselfconsistent energy for the given xcname (PBE,RPBE,PW91,VWN).
3.5 Defining a spinpolarized calculation
One can use the method SetSpinPolarized(True) to tell dacapo that the calculation should be spinpolarized.
One also need to specify an initial magnetic moment for the atoms. You can do this in the construction of the atoms, by using:
>>> from ASE import Atom, ListOfAtoms >>> atoms = ListOfAtoms([Atom('O', (0, 0, 0),magmom=1.0), ... Atom('O', (0, 0, 1.2),magmom=1.0)])
This can also be set using the ListOfAtoms method SetMagneticMoments((1,1)) or setting the SetMagneticMoment method on each atom.
This will define a initial total magnetic moment of 2 Bohr magneton, for the O2 molecule, this being also the selfconsistent result. You could also use, SetMagneticMoments((1,1)), this will define a 'antiferro' magnetic state with total magnetic moment of zero. If this case probably the end result is still the same (2 Bohr magneton) but the number of selfconsistent iterations will increase.
3.6 Kpoints
The method SetBZKPoints(kpts) can be used to define the kpoints in the BZ. In general kpts are a python list of kpoints defined in units of the reciprocal cell. You can however specify a 10x10x10 MonkhorstPack 3dimensional set simple by using SetBZKpoint((10,10,10)).
Below is a list of the special kpoints sets that can be used:
 MonkhorstPack:
 A 3dimensional MonkhorstPack special kpoint set. (see H.J.Monkhorst and J.D.Pack, Physical Review B, vol. 13, page 5188, 1976. You can defined this kpoints set simply using, SetBZKPoints((n1,n2,n3)). This will define a n1 by n2 by n3, Monkhorst Pack special kpoint set. The size along each reciprocal vector given by n1, n2 and n3, respectively.
 ChadiCohen:
A ChadiCohen (see Chadi and Cohen, Phys. Review B., vol 8,page 5747) 1dimensional kpoint can be defined using the ASE Utilities module ChadiCohen:
>>> from ASE.Utilities.ChadiCohen import CC6_1x1 >>> calc.SetBZKPoints(CC6_1x1)
This will define a 6 kpoint 1x1 ChadiCohen special kpoint set. In general the ChadiCohen sets are named CC+'<Nkpoints>'+_+'shape'. For example an 18 k point sq(3)xsq(3) is named 'CC18_sq3xsq3'.
You can use the method GetBZKPoints, to access the kpoints just defined. The method GetIBZKPoints returns the symmetry reduced set of kpoints.
3.7 Pseudopotentials
For an overview of the pseudopotentials for dacapo, see the Pseudopotential library page.
You can change or view the pseudopotentials use for an element using the methods:
 GetPseudoPotential(<elementnumber>):
 View the pseudopotential path for element elementnumber.
 SetPseudoPotential(<elementnumber>, path):
 Set the pseudopotential path for element elementnumber.
3.8 Chargedensity mixing
The method SetChargeMixing(value) can be used to enable or disable charge density mixing.
Using SetChargeMixing(True) the program will use Pulay mixing for the density (default).
You can change the Pulay density mixing coeffficient (default 0.1) and the Pulay mixing history (default 10) using:
calc.SetNetCDFEntry("ChargeMixing",att="Pulay_DensityMixingCoeff",value=0.1) calc.SetNetCDFEntry("ChargeMixing",att="Pulay_MixingHistory",value=10)
Using SetChargeMixing(False) correspond to running a calculation using the HarrisFoulkes functional using the input density. See also the Harris calculation example.
The method SetKerkerPreconditioning(value) can be used to set Kerker preconditioning for the density mixing. For value=True Kerker preconditiong is used, i.e. q is different from zero, see eq. 82 in Kresse/Furthmller: Comp. Mat. Sci. 6 (1996). The value of q is fix to give a damping of 20 of the lowest q vector. For value=False q is zero and mixing is linear (default).
3.9 Eigenvalue solver
The eigenvalue solver can be set using the method SetEigenvalueSolver(method). Here method is one of:
 eigsolve: Block Davidson algorithm (Claus Bendtsen et al). This is the default method.
 rmmdiis: Residual minimization method (RMM), using DIIS (direct inversion in the iterate subspace). The implementation follows closely the algorithm outlined in Kresse and Furthmuller, Comp. Mat. Sci, III.G/III.H
3.10 Dipolecorrection
The method:
DipoleCorrection(MixingParameter=0.2,InitialValue=0.0,AdditiveDipoleField=0)
will add the interslab dipole electrostatic decoupling. Presently, this feature is only active along the third unit cell vector, so you need to choose unit cell accordingly. Corrections for higher electrostatic multipoles are not implemented presently.
The default position of the field discontinuity at the vacuum position farthest from any other atoms on both sides of the slab.
The keyword AdditiveDipoleField, will add a constant electrostatic field along third unit cell vector, corresponding to an external dipole layer. This field is added to the dipole correction field. The field discontinuity follows the position of the dynamical dipole correction.
See also the Dipole correction example and the workfunction pdf document. For more information see the dacapo netcdf manual.
A fixed external field i.e. not a field added on top the the dipole correction, can be set using the NetCDF entry ExternalDipolePotential. This can be done like:
calc.SetNetCDFEntry("ExternalDipolePotential",value=0.5)
corresponding to an external field of 0.5eV/A.
3.11 Setting nondefault executable
 GetJobType:
 SetJobType can be used to view the current executable for dacapo.
 SetJobType(executable=newexecutable):
 SetJobType can be used to change the default dacapo executable dacapo.run to newexecutable.
3.12 Symmetry
 SetSymmetryOff:
 Do not use symmtry then reducing the kpoint set.
 SetSymmetryOn:
 Use symmtry then reducing the kpoint set.
3.13 Electrostatic decoupling
 SetElectrostaticDecoupling(numberofgaussians=3,ecutoff=100):
 Add electrostatic decoupling (Jan Rossmeisl). Decoupling activates the three dimensional electrostatic decoupling. Based on paper by Peter E. Bloechl: JCP 103 page7422 (1995).
3.14 Stress calculation
By default dacapo will not calculate the stress acting on the unitcell, you can use the method:
CalculateStress()
to tell dacapo to calculate the stress.
3.15 Local density of states
Dacapo can calculate the density of states projected onto atomic orbitals.
Warning
The density of states are not implemented for planewave parallel calculations, so the number of nodes must match the number kpoints (after they have been reduce by symmetry). So if you for example have 4 IBZ kpoints, you can use 4 or 2 parallel nodes and still get the projected density of states.
Typically one would do:
atoms = ListOfAtoms(..) calc = Dacapo(..) calc.CalculateAtomicDOS() atoms.SetCalculator(calc) energy = atoms.GetPotentialEnergy() dos = calc.GetLDOS(atoms=[1],angularchannels=["s"]) dos.GetPlot()
This will return a density of states projected onto the sorbital of the first atom (atom numbers starts from 1).
The method:
CalculateAtomicDOS(energywindow=None)
tells the fortran program to calculate the local density of states projected onto atomic orbitals, using energywindow for calculating the energy resolved LDOS in a specific range (must be within the eigenvalue spectrum).
The method:
``GetLDOS(**keywords )``:
returns an instance of AtomProjectedDOSTool (see below).
Keyword for the GetDOS method:
atoms
List of atoms numbers (default is all atoms)
angularchannels
List of angular channel names. The full list of names are: 's', 'p_x', 'p_y', 'p_z','d_zz','dxxyy', 'd_xy', 'd_xz', 'd_yz'. 'p' and 'd' can be used as shorthand for all p and all d channels respectively. (default is all d channels)
spin
List of spins (default all spins)
cutoffradius
'short' or 'infinite'. For cutoffradius = 'short' the integrals are truncated at 1 Angstrom. For cutoffradius = 'infinite' the integrals are not truncated. (default 'infinite')
The resulting DOS are added over the members of the three list (atoms,angularchannels and spin).
GetDOS returns a instance of the class AtomProjectedDOSTool (Thanks to John Kitchen for providing grace plot and band moments methods).
Methods for AtomProjectedDOSTool:
 GetPlot
 returns a GnuPlotAvartar for the energy resolved DOS A parent can be given too combine plot.
 GetIntegratedDOS
 returns the integral up to the fermi energy
 GetData
 returns the energy resolved DOS
 GetBandMoment(moments)
 returns the moments of the projected DOS. GetBandMoment(0) return the integrated DOS. GetBandMoment(1,2) returns the first and second moment of the projected DOS (center and width)
 SaveData(filename)
 saves the data to the filename
 GetGracePlot
 makes a GracePlot (if GracePlot is installed)
Note
The eigen values are not relative to the Fermi level, you can get the fermi level for the calculation using calc.GetFermiLevel()
You add the multicenter orbitals using:
mcenters = [(0,0,0,1.0),(1,0,0,1.0)] calc.SetMultiCenters(mcenters)
this will define an atom1_satom2_s orbital.
The general form is [{(atomnr,l,m,weight)}]
You can retrieve the centers using:
centers = calc.GetMultiCenters()
This will return a list of instances of the class MultiCenterProjectedDOSTool.
This class has the same method as described above, ie. GetBandMoment, `GetCutoff, GetData, GetDescription, GetEFermi, GetGracePlot, GetIntegratedDOS, GetPlot, GetSpin`and SaveData.
You can give the following keywords to SetMultiCenters:
calc.SetMultiCenters(mcenters, energywindow=(15,5), energywidth=0.2, numberenergypoints=250, cutoffradius=1.0)
The direction cosines for which the spherical harmonics are set up are using the next different atom in the list (cyclic) as direction pointer, so the zdirection is chosen along the direction to this next atom. At the moment the rotation matrices is only given in the text file, you can use grep 'MUL: Rmatrix' out_o2.txt to get this information.
Below is a small example, which defines a sorbital on atom 0:
#!/usr/bin/env python from Dacapo import Dacapo from ASE import Atom, ListOfAtoms atoms = ListOfAtoms([Atom('O', (2, 2, 2) ,magmom=1.0), Atom('O', (3.1, 2, 2),magmom=1.0)], cell=(4, 4, 4)) calc = Dacapo(planewavecutoff=350, # in eV nbands=8, # 5 extra empty bands out='o2.nc', spinpol=True ) atoms.SetCalculator(calc) mcenters = [] mcenters.append([(0,0,0,1.0),(1,0,0,0.0)]) # [{(atomnr,l,m,weight)}] calc.SetMultiCenters(mcenters,energywidth=0.1) energy = atoms.GetPotentialEnergy() centers = calc.GetMultiCenters() center1 = centers[0] print 'description:' center1.GetDescription() print 'energy resolved dos:' for d in center1.GetData(): print d[0],d[1]
3.16 Changing convergence parameters
The method:
SetConvergenceParameters(**keywords)
can be used to set convergence parameters for the calculations. The following keywords are defined:
 energy:
Absolute convergence in energy (default 1e5)
 density:
Convergence in density (default 0.0001)
 occupation:
Convergence in occupation numbers (default 0.001)
SetConvergenceParameters(energy=0.00001,density=0.0001,occupation=0.001) correspond to the default settings.
 The method::
 GetConvergenceParameters()
can be used to retrieve the convergence parameters.
3.17 dacapo fortran program/python interface communication
The method:
StayAliveOn()
can be used for normal optimization runs.
This will prevent the dacapo fortran program to terminate and write the full netcdf output file after each call of GetPotentialEnergy. Instead it will wait for python to send new positions. This way of running the program saves IO time and reduce the NFS traffic.
Do not use this mode for NEB calculations, you will end of having 10 fortran programs in memory.
3.18 NetCDF specific methods
Two methods are defined for general reading and writing of NetCDF entries. This provides a method for reading and writing any NetCDF variable defined in the dacapo netcdf manual, so that any settings for the dacapo fortran program, not defined by the standard sets of method from the DFT Calculator or from the extra Dacapo method described above, can be defined. The methods are listed below:
 SetNetCDFEntry(variable,att=None,value=None):
 Define a general netcdf entry for this calculation. If the netcdf entry is allready defined, the value are modified.
 GetNetCDFEntry(variable,att=None):
 Returns the value of NetCDF entry <variable> or NetCDF attribute att, if specified.
The effect of writing:
>>> calc = Dacapo() >>> calc.SetNumberOfBands(10)
and
>>> calc = Dacapo() >>> calc.SetNetCDFEntry('ElectronicBands',att='NumberOfBands',value=10)
is the same. Now you can get the information by either using
>>> print calc.GetNumberOfBands()
or
>>> print calc.GetNetCDFEntry('ElectronicBands',att='NumberOfBands')
4 Getting information from the calculator
Below follow a list of the methods one can use to read information from the calculator:
 GetNumberOBands:
 Return the number of electronic bands.
 GetXCFunctional:
 Return the XCfunctional identifier ('LDA','PBE' ..)
 GetBZKPoints
 Return the Kpoints in the in the Brillouin zone. The Kpoints are in units of the reciprocal space basis.
 GetSpinPolarized:
 Is it a spinpolarized calculation.
 GetMagneticMoment:
 Return the total magnetic moment in units of Bohr magnetons.
 GetIBZKPoints:
 Return kpoints in the irreducible part of the irreducible Brillouin zone. The Kpoints are in units of the reciprocal space basis
 GetIBZKPointWeights
 Weights of the kpoints. The sum of all weights is one.
 GetEigenvalues(kpt=0, spin=0):
 Return the eigenvalues for all bands for the given Kpoint kpt and the given spin spin. Eigenvalues are relative to the Fermilevel.
 GetOccupationNumbers(kpt=0, spin=0):
 Return the occupation numbers for all bands for the given Kpoint kpt and the given spin spin. These numbers fall in the range [0,1] for a spinpolarized calculation, and in the range [0,2] otherwise.
 GetWavefunctionArray(band, kpt=0, spin=0):
Return array of wavefunction values. The wavefunction are returned for the band band, Kpoint kpt and spin spin.
The wavefunction is returned as a threedimensional Numerical array, corresponding to the grid use in the calculation.
 GetWavefunctionArrays(kpt=0, spin=0):
 Return array of wavefunction values for all bands for the given Kpoint kpt and the given spin spin. The wavefunctions are returned as a threedimensional Numerical array.
 GetDensityArray(spin=0):
 Return array of density values for the given spin spin.
5 Working with the Dacapo calculator
5.1 reading from an old calculation
If you want to restart a calculation or read specific information from an old calculation, one have to access the information via the atoms, using:
>>> atoms = Dacapo.ReadAtoms(filename='oldfile')
Now you can get the calculator, attached to the atoms, and get the information, using:
>>> calc = atoms.GetCalculator() >>> print calc.GetEigenValues()
Often ASE methods take the atoms as the argument, using the GetCalculator method to access the relevant information, i.e the VTK methods for 3D plotting:
>>> from ASE.Visualization.VTK import VTKPlotWaveFunction >>> atoms = Dacapo.ReadAtoms(filename='oldfile') >>> VTKPlotWaveFunction(atoms)
6 Minimization and Dynamics using python
Minimization of structures is done using the mimimization algorithms implemented in ASE, for more details see the ASE manual.
7 Wavefunction and density plots
See the example plot wavefunction.
8 Wannier orbitals
Using the ASE Wannier module one can find the most localized set of Wannier orbitals for a dacapo calculation. More details will be added. See the Wannier example.
9 Nudged Elastic Band
See the Nudged Elastic Band method for a general description. See the Nudged Elastic Band example for a dacapo application of this method.