Software - 岳 岭 - 教师个人主页.pdf

DFTB+ Version 19.1 USER MANUAL 2 Contents 1 Introduction 7 2 Input for DFTB+ 2.1 Main input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Explicit geometry specification . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 GenFormat{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Driver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 SteepestDescent{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 ConjugateGradient{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 gDIIS{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 LBFGS{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.5 SecondDerivatives{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.6 VelocityVerlet{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.7 Socket{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Mixer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 SpinPolarisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 SpinOrbit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.4 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.5 Filling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.6 SlaterKosterFiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.7 KPointsAndWeights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.8 OrbitalPotential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.9 ElectricField . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.10 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.11 DFTB3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.12 Hydrogen corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.13 RangeSeparated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.14 On site corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.15 Differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.16 ForceEvaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 ExcitedState . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.1 Casida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 ParserOptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Parallel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 10 10 10 11 12 12 14 15 15 16 16 22 23 29 31 35 35 38 39 40 42 43 44 49 49 52 53 53 54 54 56 59 59 61 62 3 4 CONTENTS 3 Output of DFTB+ 3.1 band.out . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 detailed.out . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 results.tag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 hamsqrN.dat, oversqr.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 hamrealN.dat, overreal.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 eigenvec.out, eigenvec.bin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 charges.bin / charges.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 md.out . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Electrostatic potential data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Excited state results files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.1 ARPACK.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.2 COEF.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.3 EXC.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.4 SPX.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.5 TDP.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.6 TRA.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.7 TEST_ARPACK.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.8 XCH.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.9 XplusY.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10.10 XREST.DAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 63 63 64 65 65 66 66 67 67 68 68 68 69 69 70 70 70 70 70 71 4 Transport calculations 4.1 Definition of the geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Device{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Contact{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Task = ContactHamiltonian{} . . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Task = UploadContacts{} . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 GreensFunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Solver = TransportOnly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Contour integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Spin-polarised transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Poisson solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7.2 Electrostatic Gates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Parallelisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 TunnelingAndDos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Setting electronic temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Troubleshooting transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Transport Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 73 75 75 76 76 78 79 81 81 82 82 84 87 87 88 88 89 90 90 5 93 93 94 94 MODES 5.1 Input for MODES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Hessian{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 DisplayModes{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CONTENTS 5 6 WAVEPLOT 95 6.1 Input for WAVEPLOT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.1.1 Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.1.2 Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.1.3 ParserOptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7 SETUPGEOM 7.1 8 103 Input for SETUPGEOM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.1.1 Transport{} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.1.2 Code output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 DFTB+ API 107 8.1 Building the library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 8.2 General guidelines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 A The HSD format A.1 Scalars and list of scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Methods and property lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Modifiers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 File inclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Extended format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 110 111 112 112 113 113 B Unit modifiers 117 C Description of the gen format 119 D Atomic spin constants 121 E Slater-Kirkwood dispersion constants 123 F DftD3 dispersion constants 125 G Atomic onsite constants 127 H Description of restart files 129 H.0.1 charges.bin / charges.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 I Publications to cite Index 131 137 6 CONTENTS Chapter 1 Introduction The code DFTB+ is the Fortran 2003 successor of the old DFTB code, which implements the density functional based tight binding approach [15]. The code has been completely rewritten from scratch and extended with various features. The main features of DFTB+ are: • Non-scc and scc calculations (with an expanded range of SCC accelerators) – Cluster/molecular systems – Periodic systems (arbitrary k-point sampling, band structure calculations, etc.) – Open boundary conditions (wire and semi-infinite contacts) • l-shell resolved calculations possible • Spin polarised calculations (both collinear and non-collinear spin) • Geometry and lattice optimisation • Geometry optimisation with constraints (in xyz-coordinates) • Molecular dynamics (NVE, NPH, NVT and NPT ensembles) • Numerical vibrational mode calculations • Finite temperature calculations • Dispersion corrections (van der Waals interactions) • Ability to treat f -electrons • Linear response excited state calculations for cluster/molecular systems • Geometry optimisation and molecular dynamics in singlet and triplet excited states of spinfree molecules • LDA+U/pSIC extensions • L · S coupling • 3rd order on-site corrections (improved H-bonding) 7 8 CHAPTER 1. INTRODUCTION • Long range hybrid corrections • QM/MM coupling with external point charges (smoothing possible) • OpenMP and MPI parallelisation with a range of electronic solvers • Electronic transport calculations (non-equilibrium Green function methodology, transmission calculations) • More general electrostatic boundary conditions via a Poisson equation electrostatic solver • Automatic code validation (autotest system) • New user friendly, extensible input format (HSD or XML) • Dynamic memory allocation • Additional tool for generating cube files for charge distribution, molecular orbitals, etc. (Waveplot) Chapter 2 Input for DFTB+ DFTB+ can read two formats, either XML or the Human-friendly Structured Data format (HSD). If you are not familiar with the HSD format, a detailed description is given in appendix A. The input file for DFTB+ must be named dftb_in.hsd or dftb_in.xml. The input file must be present in the working directory. To prevent ambiguity, the parser refuses to read any input if both files are present. After processing the input, DFTB+ creates a file of the parsed input, either dftb_pin.hsd or dftb_pin.xml. This contains the user input as well as any default values for unspecified options. The file also contains the version number of the current input parser. You should always keep this file, since if you want to exactly repeat your calculation with a later version of DFTB+, it is recommended to use this file instead of the original input. (You must of course rename dftb_pin.hsd into dftb_in.hsd or dftb_pin.xml into dftb_in.xml.) This guarantees that you will obtain the same results, even if the defaults for some non specified options have been changed. The code can also produce dftb_pin.xml from dftb_in.hsd or vice versa if required (see section 2.8). The following sections list properties and options that can be set in DFTB+ input. The first column of each of the tables of options specifies the name of a property. The second column indicates the type of the expected value for that property. The letters “l”, “i”, “r”, “s”, “p”, “m” stand for logical, integer, real, string, property list and method type, respectively. An optional prefixing number specifies how often (if more than once) this type must occur. An appended “+” indicates arbitrary occurrence greater than zero, while “*” allows also for zero occurrence. Alternative types are separated by “|”. Parentheses serve only to delimit groups of settings. Sometimes a property is only interpreted on meeting some condition(s). If this is the case, the the third column gives details of the requirement(s). The fourth column contains the default value for the property. If no default value is specified (“-”), the user is required to assign a value to that property. The description of the properties immediately follows the table. If there is also a more detailed description available for a given keyword somewhere else, the appropriate page number appears in the last column. Some properties are allowed to carry a modifier to alter the provided value (e.g. converting between units). The possible modifiers are listed between brackets ([]) in the detailed description of the property. If the modifier is a conversion factor for a physical unit, only the unit type is indicated (length, energy, force, time, etc.). A list of the allowed physical units can be found in appendix B. 9 10 2.1 CHAPTER 2. INPUT FOR DFTB+ Main input The input file for DFTB+ (dftb_in.hsd/dftb_in.xml) must contain the following property definitions: Name Geometry Hamiltonian Type p|m m Condition Default - Page 10 23 Additionally optional blocks of definitions may be present: Name Driver Options Analysis ExcitedState ParserOptions Parallel Type m p p p p p Condition Default {} {} {} {} {} {} Page 12 54 56 59 61 62 Geometry Specifies the geometry for the system to be calculated. See p. 10. Hamiltonian Configures the Hamiltonian and its options. See p. 23. Driver Specifies a geometry driver for your system. See p. 12. Options Various global options for the run. See p. 54. Analysis Post-run analysis and properties options. See p. 56. ExcitedState Calculations in excited state of the system. See p. 59. ParserOptions Various options affecting the parser only. See p. 61. Parallel Options affecting the MPI-parallel execution. See p. 62. 2.2 Geometry The geometry can be specified either directly by passing the appropriate list of properties or by using the GenFormat{} method. 2.2.1 Explicit geometry specification If the geometry is being specified explicitly, the following properties can be set: Periodic LatticeVectors TypeNames TypesAndCoordinates l 9r s+ (1i3r)+ Periodic = Yes No - Periodic Specifies if the system is periodic in all 3 dimensions or is to be treated as a cluster. If set to Yes, property LatticeVectors{} must be also specified. 2.2. GEOMETRY 11 LatticeVectors [length] The x, y and z components of the three lattice vectors if the system is periodic. TypeNames List of strings with the names of the elements, which appear in your geometry. TypesAndCoordinates [relative|length] For every atom the index of its type in the TypeNames list and its coordinates. If for a periodic system (Periodic = Yes) the modifier relative is specified, the coordinates are interpreted in the coordinate system of the lattice vectors. Example: Geometry of GaAs: Geometry = { TypeNames = { "Ga" "As" } TypesAndCoordinates [Angstrom] = { 1 0.000000 0.000000 0.000000 2 1.356773 1.356773 1.356773 } Periodic = Yes LatticeVectors [Angstrom] = { 2.713546 2.713546 0. 0. 2.713546 2.713546 2.713546 0. 2.713546 } } 2.2.2 GenFormat{} You can use the generic format to specify the geometry (see appendix C). The geometry specification for GaAs would be the following: Geometry = GenFormat { 2 S Ga As 11 0.000000 0.000000 0.000000 22 1.356773 1.356773 1.356773 0.000000 0.000000 0.000000 2.713546 2.713546 0. 0. 2.713546 2.713546 2.713546 0. 2.713546 } It is also possible to include the gen-formatted geometry from a file: Geometry = GenFormat { <<< "geometry.gen" } 12 2.3 CHAPTER 2. INPUT FOR DFTB+ Driver The driver is responsible for changing the geometry of the input structure during the calculation. Currently the following methods are available: {} Static calculation with the input geometry. SteepestDescent{} Geometry optimisation by moving atoms along the acting forces. See p. 12. ConjugateGradient{} Geometry optimisation using the conjugate gradient algorithm. See p. 14. gDIIS{} Geometry optimisation using the modified gDIIS method. See p. 15. LBFGS{} Geometry optimisation using the LBFGS algorithm. See p. 15. SecondDerivatives{} Calculation of the second derivatives of the energy (the Hessian). See p. 16. VelocityVerlet{} Molecular dynamics with the velocity Verlet algorithm. See p. 16. Socket{} Hands over control to an external program via a socket interface. See p. 22. 2.3.1 SteepestDescent{} MovedAtoms MaxForceComponent MaxSteps StepSize OutputPrefix AppendGeometries Constraints LatticeOpt FixAngles FixLengths Isotropic Pressure MaxAtomStep MaxLatticeStep ConvergentForcesOnly (i|s)+ r i r s l (1i3r)* l l 3l l r r r l LatticeOpt = No Periodic = Yes Periodic = Yes, LatticeOpt = Yes FixAngles = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes, LatticeOpt = Yes MovedAtoms 6= None{} Periodic = Yes, LatticeOpt = Yes SCC = Yes 1:-1 1e-4 200 100.0 "geo_end" No {} No No No No No No 0.0 0.2 0.2 Yes MovedAtoms Indices of the atoms which should be moved. The atoms can be specified as a mixture of a list of atoms, ranges of atoms and/or the species of atoms. Index ranges are specified as start:end (without white space as one word!), which inclusively selects all atoms between start and end. MovedAtoms = 1:6 # equivalent to MovedAtoms = { 1 2 3 4 5 6 } Negative indices can be used to count backwards from the last atom (-1 = last atom, -2 = penultimate atom, etc.): 2.3. DRIVER 13 MovedAtoms = 1:-1 # Move all atoms including the last Species names can be used to select all atoms belonging to a given species: MovedAtoms = Ga # select all Ga atoms Various specifiers can be combined together: # Move atoms 1, 2, 3, all Ga atoms, and the last two atoms. MovedAtoms = 1:3 Ga -2:-1 MaxForceComponent [force] Optimisation is stopped, if the force component with the maximal absolute value goes below this threshold. MaxSteps Maximum number of steps after which the optimisation should stop (unless already stopped by achieving convergence). Setting this value as -1 runs a huge() number of iterations. StepSize [time] Step size (δt) along the forces. The displacement δ xi along the ith coordinate is fi given for each atom as δ xi = 2m δt 2 , where fi is the appropriate force component and m is the mass of the atom. OutputPrefix Prefix of the geometry files containing the final structure. AppendGeometries If set to Yes, the geometry file in the XYZ-format will contain all the geometries obtained during the optimisation (instead of containing only the last geometry). Constraints Specifies geometry constraints. For every constraint the serial number of the atom is expected followed by the x, y, z components of a constraint vector. The specified atom is not allowed to move along the constraint vector. If two constraints are defined for the same atom, the atom will only by able to move normal to the the plane containing the two constraining vectors. Example: Constraints = { # Atom one can only move along the z-axis 1 1.0 0.0 0.0 1 0.0 1.0 0.0 } LatticeOpt Allow the lattice vectors to change during optimisation. MovedAtoms can be optionally used with lattice optimisation if the atomic coordinates are to be co-optimised with the lattice.1 FixAngles If optimising the lattice, allow only the lengths of lattice vectors to vary, not the angles between them. For example if your lattice is orthorhombic, this option will maintain that symmetry during optimisation. FixLengths If optimising the lattice with FixAngles = Yes, allow only the lengths of the specified lattice vectors to vary. Example: 1 This is functional but not very efficient at the moment. 14 CHAPTER 2. INPUT FOR DFTB+ Driver = ConjugateGradient { LatticeOpt = Yes FixAngles = Yes # Fix angles between lattice vectors FixLengths = {Yes Yes No} # Allow only lat. vector 3 to change length } Isotropic If optimising the lattice, allow only uniform scaling of the unit cell. This option is incompatible with FixAngles. Pressure [pressure] If optimising the lattice, set the external pressure, leading to a Gibbs free energy of the form G = E + PV − T S being printed as well (the included entropy term is only the contribution from the electrons, therefore this is not the full free energy). MaxAtomStep Sets the maximum possible line search step size for atomic relaxation. MaxLatticeStep Sets the maximum possible line search step size for lattice optimisation. For FixAngles or Isotropic calculations this is as a fraction of the lattice vectors or the volume respectively. ConvergentForcesOnly If using an SCC calculation, this option controls whether the geometry optimisation will prematurely stop (= Yes) if the SCC cycle does not converge at any geometric step. 2.3.2 ConjugateGradient{} MovedAtoms MaxForceComponent MaxSteps OutputPrefix AppendGeometries Constraints LatticeOpt FixAngles Isotropic Pressure MaxAtomStep MaxLatticeStep ConvergentForcesOnly (i|s)+ r i s l (1i3r)* l l l r r r l Periodic = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes MovedAtoms 6= None{} Periodic = Yes, LatticeOpt = Yes SCC = Yes See previous subsection for the description of the properties. 1:-1 1e-4 200 "geo_end" No {} No No No 0.0 0.2 0.2 Yes 2.3. DRIVER 2.3.3 15 gDIIS{} Alpha Generations MovedAtoms MaxForceComponent MaxSteps OutputPrefix AppendGeometries Constraints LatticeOpt FixAngles Isotropic Pressure MaxLatticeStep ConvergentForcesOnly r i (i|s)+ r i s l (1i3r)* l l l r r l Periodic = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes Periodic = Yes, LatticeOpt = Yes SCC = Yes 0.1 8 1:-1 1e-4 200 "geo_end" No {} No No No 0.0 0.2 Yes Specific properties for this method are: Alpha Initial scaling parameter to prevent the iterative space becoming exhausted (this is dynamically adjusted during the run). Generations Number of generations to consider for the mixing. See previous subsection for the description of the other properties.2 2.3.4 LBFGS{} Memory MovedAtoms MaxForceComponent MaxSteps OutputPrefix AppendGeometries Constraints LatticeOpt FixAngles Isotropic Pressure MaxLatticeStep ConvergentForcesOnly i (i|s)+ r i s l (1i3r)* l l l r r l Periodic = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes, LatticeOpt = Yes Periodic = Yes Periodic = Yes, LatticeOpt = Yes SCC = Yes 20 1:-1 1e-4 200 "geo_end" No {} No No No 0.0 0.2 Yes Specific properties for this method are: Memory Number of last steps which are saved and used to calculate the next step via the LBFGS algorithm. The literature recommends that Memory should between 3 and 20 [41]. 2 This approach is distinct from section 2.4.1, but uses a related algorithm based on Ref. [30] and comments from P.R.Briddon. 16 2.3.5 CHAPTER 2. INPUT FOR DFTB+ SecondDerivatives{} Calculates the second derivatives of the energy (currently only using a numerical differentiation of the forces). The derivatives matrix is written out for the i, j and k directions of atoms 1 . . . n as ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ... ∂ xi1 ∂ xi1 ∂ x j1 ∂ xi1 ∂ xk1 ∂ xi1 ∂ xi2 ∂ xi1 ∂ x j2 ∂ xi1 ∂ xk2 ∂ xi1 ∂ xkn ∂ xkn into hessian.out Note: for supercell calculations, the derivatives are obtained at the q = 0 point, irrespective of the k-point sampling used. Important: In order to get accurate results for the second derivatives (and the resulting frequencies) you must set a smaller self-consistent tolerance than the default value in the Hamiltonian{} section. We suggest SCCTolerance = 1e-7 or better. A less accurate tolerance can yield nonphysical vibrational frequencies. Atoms Delta 1:-1 1e-4 i+|m r Atoms Index of the atoms for which to calculate the second derivatives. The atoms can be specified via indices, index ranges and species. (See MovedAtoms in section 2.3.1.) Delta Step size for numerical differentiation of forces to get the second derivatives of the energy with respect to atomic coordinates. 2.3.6 VelocityVerlet{} The code propagates atomic motion using velocity Verlet dynamics with optional thermostats or barostats to control the temperature and/or pressure. Information is printed out during the simulation every MDRestartFrequency steps, and logged in the file md.out (see appendix 3.8). MovedAtoms Steps TimeStep KeepStationary Thermostat OutputPrefix MDRestartFrequency Velocities Barostat ConvergentForcesOnly Xlbomd XlbomdFast Masses (i|s)+ i r l m s i (3r)* m l p p p Periodic = Yes SCC = Yes XlbomdFast not set Xlbomd not set 1:-1 Yes "geo_end" 1 Yes 17 19 20 20 22 MovedAtoms List of atoms to move during the MD. (See more detailed description on page 12.) Steps Number of MD steps to perform. In the case of a thermostat using a TemperatureProfile{} the number of steps is calculated from the profile. 2.3. DRIVER 17 KeepStationary Remove translational motion from the system. TimeStep [time] Time interval between two MD steps. Thermostat Thermostating method for the MD simulation. See p. 17. OutputPrefix Prefix of the geometry files containing the final structure. MDRestartFrequency Interval that the current geometry and velocities are written to the XYZ format geometry file and md.out (see section 3.8). In the case of SCC MD runs, the charge restart information is also written at this interval overriding RestartFrequency (see section 2.5). Velocities [velocity] Specified atomic velocities for all the atoms of the given structure (including “velocities” for any stationary atoms, which are silently ignored). This option can be used to restart an MD run, but make sure the geometry is consistent with the specified velocities. The easiest way to do this is to copy both from the same iteration of the XYZ file produced in the previous run (Note: the velocities printed in the XYZ files are specified in Å/ps, so this should be set in the input). If restarting an SCC MD run, we strongly suggest you use ReadInitialCharges, and in particular read charges for the geometry which you use to restart (iterations at which charges are written to disc are marked in the XYZ file, with the most recent being present in charges.bin or charges.dat depending on the option WriteChargesAsText). Barostat Berendsen method barostat for the MD simulation. See p. 19. ConvergentForcesOnly If using an SCC calculation, this option controls whether the molecular dynamics will prematurely stop (= Yes) if the SCC cycle does not converge at any geometric step. If the option is set to False, forces will be calculated using the non-converged charges and the molecular dynamics continues. In this case you should consider using ForceEvaluation = ’Dynamics’ (or ForceEvaluation = ’DynamicsT0’) in the DFTB block as it gives more accurate forces for non-converged charges. Xlbomd If present, extended Lagrangian type molecular dynamics is applied to speed up the simulation. For further options within the Xlbomd block see p. 20. Masses If present, over-ride the atomic masses from the Slater-Koster files. See p. 22 Thermostat None{} No thermostating during the run, only the initial velocities are set according to either a given temperature or velocities, hence an NVE ensemble should be achieved for a reasonable time step. InitialTemperature r - InitialTemperature [energy] Starting velocities for the MD will be created according the Maxwell-Boltzmann distribution at the specified temperature. This is redundant in the case of specified initial velocities. Andersen{} Andersen thermostat [3] sampling an NVT ensemble. Note: Andersen thermostating has a reputation for suppressing diffusion and also prevents accumulation of data for dynamical properties. 18 Temperature ReselectProbability ReselectIndividually AdaptFillingTemp CHAPTER 2. INPUT FOR DFTB+ No r|m r l l Temperature [energy] Target temperature of the thermostat. It can be either a real value, specifying a constant temperature through the entire run or the TemperatureProfile{} method specifying a changing temperature. (See p. 19.) ReselectProbability Probability for re-selecting velocities from the Maxwell-Boltzmann distribution. ReselectIndividually If Yes, each atomic velocity is re-selected individually with the specified probability. Otherwise all velocities are re-selected simultaneously with the specified probability. AdaptFillingTemp If Yes, the temperature of the electron filling is always set to the current temperature of the thermostat. (The appropriate tag for the temperature of the electron filling is ignored.) Berendsen{} Berendsen thermostat [7] samples something like an NVT ensemble (but without the correct canonical fluctuations, be aware of the “flying ice cube” problem before using this thermostat [22]). Temperature CouplingStrength Timescale AdaptFillingTemp r|m r r l No Timescale not set CouplingStrength not set Temperature [energy] Target temperature of the thermostat. It can be either a real value specifying a constant temperature through the entire run or the TemperatureProfile{} method specifying a changing temperature. (See p. 19.) CouplingStrength Dimensionless coupling strength for the thermostat (given as ∆t/τt in the original paper where ∆t is the MD time step). Alternatively Timescale[time] can be set directly as the characteristic length of time to damp the temperature towards the target temperature. The CouplingStrength and Timescale options are mutually exclusive. AdaptFillingTemp If Yes, the temperature of the electron filling is always set to the current temperature of the thermostat. (The appropriate tag for the temperature of the electron filling is ignored.) NoseHoover{} Nosé-Hoover chain thermostat [34] sampling an NVT ensemble, currently with the chain coupled to all of the atoms in the system. Temperature CouplingStrength ChainLength Order IntegratorSteps Restart AdaptFillingTemp r|m r i i i m l 3 3 1 No 2.3. DRIVER 19 Temperature [energy] Target temperature of the thermostat. It can be either a real value, specifying a constant temperature through the entire run or the TemperatureProfile{} method specifying a changing temperature. (See p. 19, but note that profiles are not well tested with this thermostat yet.) CouplingStrength [Frequency] Frequency of oscillation of the thermostating particles (see section 2.5 of Ref. [34]). This is typically related to the highest vibrational mode frequency of the system. ChainLength Number of particles in the thermostat chain. Order and IntegratorSteps See section 4.3 of Ref. [34]). Restart Specifies the internal state of the thermostat chain, using three keywords (all three must be present): x{}, v{} and g{} containing the internal chain positions, velocities and forces respectively (these are provided in md.out). See also section 2.3.6. AdaptFillingTemp If Yes, the temperature of the electron filling is always set to the current temperature of the thermostat. (The appropriate tag for the temperature of the electron filling is ignored.) TemperatureProfile{} Specifies a temperature profile during molecular dynamics. It takes as argument one or more lines containing the type of the annealing (string), its duration (integer), and the target temperature (real), which should be achieved at the end of the given period. For example: Temperature [Kelvin] = TemperatureProfile { # Temperatures in K constant 1 10.0 # Setting T=10 K for the 0th MD-step linear 500 600.0 # Linearly rising T in 500 steps up to T=600 K constant 2000 600.0 # Constant T through 2000 steps exponential 500 10.0 # Exponential decreasing in 500 steps to T=10 K } The annealing method can be constant, linear or exponential, with the duration of each stage of the anneal specified in steps of the driver containing the thermostat. The starting temperature for each annealing period is the final target temperature of the previous period, with the last step of each stage being at the target temperature. Since the initial stage in the temperature profile has no previous step, the default starting temperature is set to 0, but no actual calculation for that temperature is made. In order to start the calculation from a finite temperature for the first annealing period, a constant profile temperature stage with the duration of only one step should be specified as the first step (see the example). The temperatures of the stages are specified in atomic units, unless the Temperature keyword carries a modifier, as in the example above. Barostat Berendsen barostat [7] samples something like an NPH ensemble for MD (but without the correct fluctuations). Options are provided for either isotropic or cell shape changing pressure control. This can also be used in tandem with a thermostat (p. 17) for the NPT ensemble. If the barostat is used, a partial Gibbs free energy is reported in code output, of the form G = E + PV − T Selectronic . 20 CHAPTER 2. INPUT FOR DFTB+ This does not include structural entropy, only any electronic entropy. For barostated constant energy simulations (no thermostat in use), the conserved quantity is the sum of the kinetic and Gibbs partial energies. Pressure Coupling Timescale Isotropic r r r l Yes Timescale not set Coupling not set Pressure [pressure] Sets the external target pressure. Coupling Coupling strength for the barostat (given as β ∆t/τ p in the original paper where ∆t is the MD time step and β the compressibility, so the coupling is technically dimensioned as reciprocal pressure, but this is usually ignored). Alternatively Timescale[time] can be set directly (β /τ p ) as the characteristic length of time to damp the pressure. Typically, β is assumed to be either the experimental value or ∼ 1, and Ref. [7] chooses the time scale to be around 10–100 fs. The Coupling and Timescale options are mutually exclusive. Isotropic Should isotropic scaling of the unit cell be used, or can the cell shape vary? There is a slight inconsistency between the standard forms of these scalings in the literature, which is reproduced here, in brief the isotropic case scales the cell volume by a factor proportional to the differences in the instantaneous and expected pressures (i.e., the cube of the cell vectors), while the anisotropic case changes the cell vectors proportional to the difference instead. Extended Lagrangian Born-Oppenheimer dynamics For several systems Born-Oppenheimer molecular dynamics simulations can be significantly sped up by using the extended Lagrangian formalism described in Ref. [6]. The XLBOMD integrator can be used in two different modes: • Conventional XLBOMD scheme (Xlbomd): The extended Lagrangian is used to predict the input charge distribution for the next time step, instead of taking charges that were converged for the geometry in the previous time step. The predicted starting charges should then require fewer SCC iterations to converge. • Fast XLBOMD scheme, XlbomdFast (one diagonalisation per time step): The extended Lagrangian is used to predict the population for each time step. This predicted population is then used to build the Hamiltonian, but in contrast to the conventional XLBOMD scheme, there is no self consistent cycle – the forces are calculated immediately after the diagonalisation of the first Hamiltonian. The fast XLBOMD method usually only works for systems without SCC instabilities (e.g. wider gap insulators or molecules without degenerate states). See Ref. [6] for details. The extended Lagrangian dynamics can be activated by specifying either the Xlbomd or the XlbomdFast option block. Both methods offer the following options: IntegrationSteps PreSteps i i 5 0 IntegrationSteps Number of time steps used for determining the population for the next time step. Currently, only integration schemes for 5, 6 or 7 steps are implemented. 2.3. DRIVER 21 PreSteps Number of molecular dynamics time steps before the XLBOMD integration becomes activated. Note: At the step where the XLBOMD integrator becomes active, it is initialised with results of several subsequent converged MD steps, so a further IntegrationSteps + 1 steps will be carried out before the extended Lagrangian dynamics starts to predict the charges and accelerate the calculation. The conventional Xlbomd block has the following specific options in addition to the common XLBOMD settings above: MinSccIterations MaxSccIterations SccTolerance i i r 1 200 1e-5 MinSccIterations Minimum number of SCC iterations to perform at each time step during the extended Lagrangian dynamics. MaxSccIterations Maximum number of SCC iterations to perform at each step in the extended Lagrangian dynamics. If this number of SCC iterations have been reached the forces will be calculated and the MD advances to the next time step. See the note in section 2.4.7 regarding non-convergent k-point sampling. SccTolerance SCC convergence tolerance during the extended Lagrangian dynamics. Once this tolerance has been achieved the SCC cycle will stop and the forces will be calculated. You can use this parameter to override the SccTolerance parameter in the DFTB block for time steps where the extended Lagrangian integrator is active (This way, you can calculate populations with tight SCC tolerance when initialising the XLBOMD integrator, then use a less strict SCC tolerance once the integrator is active). The XlbomdFast block has the following specific options in addition to the common XLBOMD settings above: TransientSteps Scale i r 10 1.0 TransientSteps Enables a smoother transition between Born-Oppenheimer and extended Lagrangian dynamics by carrying out intermediate additional steps with full SCC convergence, during which the converged population and the one predicted by the extended Lagrangian integrator are averaged. Scale Scaling factor for the predicted charge densities ∈ (0, 1]. The optimal value is system dependent. One should take the highest possible value that still produces stable dynamics (good conservation of energy). Example for conventional XLBOMD: Xlbomd { IntegrationSteps = 6 MinSccIterations = 2 22 CHAPTER 2. INPUT FOR DFTB+ MaxSccIterations = 200 SccTolerance = 1e-6 } Fast (SCC-free) XLBOMD with one diagonalisation per time step: XlbomdFast { PreSteps = 5 TransientSteps = 10 IntegrationSteps = 5 Scale = 0.5 } Points to be aware of: • The extended Lagrangian (especially in the fast mode) needs special force evaluation giving more accurate forces for non-convergent charges. Therefore you must set the ForceEvaluation option to ’Dynamics’ (for simulations with finite electronic temperature) or to ’DynamicsT0’ (for simulations at 0 K electronic temperature) in the DFTB block (see p. 54). • The extended Lagrangian implementation only works for the (N,V, E) ensemble so far, so neither thermostats nor barostats are allowed. • The extended Lagrangian implementation currently cannot be used for spin-polarised or spinorbit systems, or when electron filling methods other than Fermi{} filling (see p. 38) are used. Masses Provides values of atomic masses for specified atoms, ranges of atoms or chemical species. This is useful for example to set isotopes for specific atoms in the system. Mass p Any atoms not given specified masses will use the default values from the appropriate homonuclear Slater-Koster file. An example is given below Masses { Mass { Atoms = 1:2 MassPerAtom [amu] = 2.0 } } where Atoms specifies the atom or atoms which each have a mass of MassPerAtom assigned. 2.3.7 Socket{} The code tries to connect to a socket interface to receive control instructions from an external driver code. 2.4. HAMILTONIAN File Prefix Host Port Verbosity Protocol MaxSteps s s s i i m i Host not set Host not set File not set File is set 23 “/tmp/ipi_” for Protocol = i-PI{} 0 i-PI{} 200 File Name of UNIX style file socket to connect to. Prefix Prefix to the file name, in the case of i-PI dialogue, the defaults to the path and file start that i_PI expects. Host Name or ip address of internet host to connect to (“localhost” also possible). Port Port of host to connect to. Verbosity Level of port traffic to document. Protocol Choice of message protocol over the socket connection (only communication with i-PI [8] is currently supported). MaxSteps Number of geometry steps before termination of the DFTB+ instance. Setting this value as -1 runs a huge() number of iterations. Examples First an ip address connection: Driver = Socket { Host = localhost Port = 21012 # port number Verbosity = 0 # minimal verbosity Protocol = i-PI {} # i-PI interface MaxSteps = -1 # Run indefinitely } Then a UNIX socket via the /tmp file system Driver = Socket { File = "dftb" # The protocol defines a default path in this case Protocol = i-PI {} # i-PI interface MaxSteps = 1000 # Terminate this instance after 1000 steps } Please note that this driver changes the default behaviour of some options to remove (usually unneeded) file writing: WriteDetailedOut = No and WriteBandOut = No. 2.4 Hamiltonian Currently only a DFTB Hamiltonian is implemented, so you must set Hamiltonian = DFTB{}. The DFTB{} method may contain the following properties: 24 SCC SCCTolerance MaxSCCIterations EwaldParameter EwaldTolerance ShellResolvedSCC Mixer MaxAngularMomentum Charge SpinPolarisation SpinConstants ShellResolvedSpin SpinOrbit Solver Filling SlaterKosterFiles OldSKInterpolation PolynomialRepulsive KPointsAndWeights OrbitalPotential ReadInitialCharges InitialCharges ElectricField Dispersion HCorrection ThirdOrder ThirdOrderFull RangeSeparated HubbardDerivs OnSiteCorrection Differentiation ForceEvaluation CustomisedHubbards CustomisedOccupations TruncateSKRange CHAPTER 2. INPUT FOR DFTB+ l r i r r l m p r m p l m m m p|m l p|m (4r)+|m m l p p m m l l p p p m s p p p No SCC = Yes 1e-5 SCC = Yes 100 Periodic = Yes SCC = Yes 0.0 Periodic = Yes SCC = Yes 1e-9 SCC = Yes No SCC = Yes Broyden{} 0.0 SCC = Yes {} SpinPolarisation 6= {} ShellResolvedSCC = No No SpinPolarisation 6= Colinear{} {} RelativelyRobust{} Fermi{} No {} Periodic = Yes SpinPolarisation 6= {} {} SCC = Yes No SCC = Yes {} SCC = Yes {} {} SCC = Yes None {} SCC = Yes No SCC = Yes No None ThirdOrder(Full) = Yes SCC = Yes FiniteDiff "Legacy" SCC = Yes SCC = Yes 29 31 34 35 35 38 39 40 42 43 44 49 49 52 53 53 SCC If set to Yes, a self consistent charge (SCC) calculation is made. SCCTolerance Stopping criteria for the SCC. Specifies the tolerance for the maximum difference in any charge between two SCC cycles. MaxSCCIterations Maximal number of SCC cycles to reach convergence. If convergence is not reached after the specified number of steps, the program stops. However in cases where the calculation is not for a static structure (so Driver 6= {}), this behaviour can be overridden using ConvergentForcesOnly. EwaldParameter Sets the dimensionless parameter α in the Ewald electrostatic summation for periodic calculations. This controls the fraction of the Ewald summation occurring in real and reciprocal space. Setting it to zero or negative activates an automatic determination of this parameter (default and recommended behaviour). Setting it positive forces the code to 2.4. HAMILTONIAN 25 use the supplied value. This is useful for very asymmetrical unit cells (typically a slab or nanowire with a huge surrounding vacuum region) since the memory demand of DFTB+ can increase dramatically in these cases (due to storage of a long range real space neighbour list). To determine a suitable value of α for such a cell, you should initially reduce the vacuum region and run a test calculation, looking for the value of the automatically chosen Ewald parameter in the standard output. This is then a suitable choice for the original cell with the large vacuum region. EwaldTolerance Sets the tolerance for Ewald summation in periodic systems. ShellResolvedSCC If set to Yes, all distinct Hubbard U values for the different atomic angular momenta shells are used, when calculating the SCC contributions. Otherwise, the value supplied for the s-shell is used for all angular momenta. Please note, that the old standard DFTB code was not orbitally resolved, so that only the Hubbard U for the s-shell was used. Please check the documentation of the SK-files you intend to use as to whether they are compatible with an orbitally resolved SCC calculation (many of the biological files do not use orbitally resolved charges), before you switch this option to Yes. Even if the Hubbard U values for different shells are the same in the SK-files, this flag would still effect your results, since when it is set to Yes, any charge transfer between atomic shells will change the energy of the system compared to when it is set to set to No. Mixer Mixer type for mixing the charges in an SCC calculation. See p. 29. MaxAngularMomentum Specifies the highest angular momentum for each atom type. All orbitals up to that angular momentum will be included in the calculation. Several main-block elements require d-orbitals, check the documentation of the SK-files you are using to determine if this is necessary. Possible values for the angular momenta are s, p, d, f. Example: MaxAngularMomentum = { Ga = "p" # You can omit the quotes around the As = "p" # orbital name, if you want. } By using the SelectedShells method it is also possible to combine shells from different SlaterKoster files together to treat atoms containing multiple shells with the same angular momentum. The shells to be picked from a particular atom type should be listed without any separating characters. The list of shells of different atom types should be separated by white spaces. Example: # Defining sps* basis for Si and C by combining sp and s shells from # Si and Si2, and C and C2, resp. MaxAngularMomentum = { Si = SelectedShells { "sp" "s" } # Si atom with sps* basis C = SelectedShells { "sp" "s" } # C atom with sps* basis } # Note, that you have to modify the Slater-Koster file definition accordingly SlaterKosterFiles = { Si-Si = "Si-Si.skf" "Si-Si2.skf" "Si2-Si.skf" "Si2-Si2.skf" 26 CHAPTER 2. INPUT FOR DFTB+ Si-C = "Si-C.skf" "Si-C2.skf" "Si2-C.skf" "Si2-C2.skf" C-Si = "C-Si.skf" "C-Si2.skf" "C2-Si.skf" "C2-Si2.skf" C-C = "C-C.skf" "C-C2.skf" "C2-C.skf" "C2-C2.skf" } If for a given atomic type you pick orbitals from more than one species, you must specify an appropriate combinations of file names for the Slater-Koster tables in SlaterKosterFiles{}. For every atom type combination nSK1 × nSk2 Slater-Koster files must be defined, where nSK1 and nSK2 are the number species combined to build up the shells of the two interacting atoms. The file names must be ordered with respect to the interacting species, so that the name for the second interacting species is changed first. Above you see an example, where an extended basis with an s∗ -orbital was generated by introducing the new species "Si2" and "C2", containing the appropriate s∗ -orbital for Si and C, resp., as only orbitals. In the case of multiple Slater-Koster files for a certain interaction, the repulsive data is read from the first specified file (e.g. Si-Si.skf, Si-C.skf, C-Si.skf and C-C.skf in the example above). The repulsive interactions in the other files are ignored. The mass for a certain species is read from the first SK-file for its homo-nuclear interaction. Non-minimal basis Slater-Koster data may be directly defined in the SK-files in future. Charge Total charge of the system in units of the electron charge. Negative values mean an excess of electrons. If the keyword FixedFermiLevel is present (see section 2.4.5), then value specified here will be ignored. SpinPolarisation Specifies if and how the system is spin polarised. See p. 31. SpinConstants Specifies the atom type specific constants needed for the spin polarised calculations, in units of Hartree. See p. 34. SpinOrbit Specifies if the system includes Russel-Saunders coupling. See p. 35 Solver Specifies which solver (eigensolver, Green’s function, etc.) to use with the hamiltonian. See p. 35. Filling Method for occupying the one electron levels with electrons. See p. 38. SlaterKosterFiles Name of the Slater-Koster files for every atom type pair combination. See 39. OldSKInterpolation If set to Yes (strongly discouraged), the look-up tables for the overlap and non-scc Hamiltonian contribution are interpolated with the same algorithm as in the old DFTB code. Please note, that the new method uses a smoother function, is more systematic, and yields better derivatives than the old one. This option is present only for compatibility purposes, and may be removed in the future. PolynomialRepulsive Specifies for each interaction, if the polynomial repulsive function should be used. for every pairwise combination of atoms it should contain a logical value, where Yes stands for the use of a polynomial repulsive function and No for a spline. If a specific pair of species is not specified, the default value No is used. Example: # Use the polynomial repulsive function for Ga-Ga and As-As interactions # in GaAs PolynomialRepulsive = { 2.4. HAMILTONIAN 27 Ga-Ga = Yes Ga-As = No # As-Ga unspecifed, therefore per default set to No As-As = Yes } If you want to apply the same setting for all species pairs, you can specify the appropriate logical value as argument of the SetForAll keyword: # Using polynomial repulsive functions for all interactions in GaAs PolynomialRepulsive = SetForAll { Yes } KPointsAndWeights [relative|absolute] Contains the special k-points to be used for the Brillouinzone integration. See p. 40. For automatically generated k-point grids the modifier should not be set. OrbitalPotential Specifies which (if any) orbitally dependant contributions should be added to the DFTB energy and band structure. See p. 42. ReadInitialCharges If set to Yes the first Hamiltonian is constructed by using the charge information read from the file charges.bin or charges.dat (depending on the option WriteChargesAsText, see section2.5). InitialCharges Specifies initial charges, either for all atoms or for only selected ones. In order to specify it for all atoms, use the keyword AllAtomCharges and supply the charge for every atom as a list of real values: InitialCharges = { AllAtomCharges = { # Specifies charge for each atom in an H2O molecule -0.88081627 # charge for atom 1 (O) 0.44040813 # charge for atom 2 (H1) 0.44040813 # charge for atom 3 (H2) } } Alternatively you can specify charges individually on atoms or species using the AtomCharge keyword. For every AtomCharge declaration you must provide a charge and the list of atoms, which should be initialised to that charge. (You can use the same format for the list of atoms, as described at the MovedAtoms keyword in the section for SteepestDescent): InitialCharges = { # Specifying charge for various species AtomCharge = { Atoms = H ChargePerAtom = 0.44040813 } AtomCharge { Atoms = O ChargePerAtom = -0.88081627 } } 28 CHAPTER 2. INPUT FOR DFTB+ Charge on atoms not appearing in any AtomCharge specification is set to be zero. ElectricField Specifies an external electric field, arising currently from either an applied field or a distribution of electrostatic charges. See p. 43. Dispersion Specifies which kind of dispersion correction to apply. See p. 44. OnSiteCorrection Used to include the on-site matrix elements of Domínguez et al. [10]. See p. 53. Differentiation Specifies how to calculate finite difference derivatives in the force routines. See p. 53. ForceEvaluation Decides which expressions are used to calculate the ground state electronic forces. See p. 54. Note: all methods give the same final forces when the charges are well converged. For non-converged charges however the resulting forces can differ considerably between methods. CustomisedHubbards Enables overriding of the Hubbard U values for given species. If the option OrbitalResolvedScc has been set to Yes, you need to specify one Hubbard U value for each atomic shell in the basis of the given atom type, otherwise only one atomic value is required. For all species not specified in this block, the value(s) found in their respective Slater-Koster files will be used. Warning: This option is for experts only! Overriding values stored in the SK-files may result in inconsistent results. Make sure you understand the consequences when using this option. Example: CustomisedHubbards { Si = 0.32 0.24 } CustomisedOccupations Enables overriding the reference neutral atom electronic occupations of the shells. Note that the atom remains neutral since a corresponding ionic counter charge is implicitly added to its core. This option can be used, for example, to simulate effective doping by setting a slightly larger or smaller number of electrons on certain atoms. Example: CustomisedOccupations{ ReferenceOccupation{ Atoms={1:30} p=2.01 } ReferenceOccupation{ Atoms={31:60} p=1.99 } } 2.4. HAMILTONIAN 29 The example above sets a filling population of +0.01e or -0.01e in the p shell of the corresponding atom indices. When the states are filled up, the electron excess or depletion results in a shift of the Fermi level in the bands. Warning: This option is for experts only! Overriding values stored in the SK-files may result in inconsistent results. Please look at the transport section of the dftb+ recipes to see an example of the correct use of this option. TruncateSKRange Enables overriding of the number of elements to be read from the SlaterKoster parameters, shortening the interaction range of atoms. Warning: This option is for experts only! Overriding values stored in the SK-files may result in inconsistent results. Make sure you understand the consequences when using this option. SKMaxDistance HardCutOff – No r l SKMaxDistance [length] Length at which to cut the overlap and non-SCC interactions for all atoms in the system. If this length is longer than the distances in the Slater-Koster files it will have no effect. HardCutOff The Slater-Koster interpolation DFTB+ produces will smoothly tail off to zero at a short distance beyond the table, which is controlled by OldSKInterpolation. If HardCutOff is set to Yes, then the distance set by SKMaxDistance includes this tail, i.e., no interactions occur beyond that distance. In the case of No this zeroing tail extends slightly beyond the HardCutOff distance. Example: TruncateSKRange = { SKMaxDistance [AA] = 4.0 HardCutOff = Yes } 2.4.1 Mixer DFTB+ currently offers the charge mixing methods Broyden{}, Anderson{}, DIIS{} (DIIS accelerated simple mixer) and Simple{} (simple mixer). Broyden{} Minimises the error function 2 m E = ω02 G(m+1) − G(m) +∑ n=1 ωn2 (n+1) − F (n) n(n+1) − n(n) (m+1) F + G , |F (n+1) − F (n) | |F (n+1) − F (n) | where G(m) is the inverse Jacobian, n(m) and F (m) are the charge and charge difference vector in iteration m. The weights are given by ω0 and ωm , respectively. The latter is calculated as c ωm = √ , F (m) · F (m) c being a constant coefficient [25]. (2.1) 30 CHAPTER 2. INPUT FOR DFTB+ The Broyden{} method can be configured using following properties: MixingParameter InverseJacobiWeight MinimalWeight MaximalWeight WeightFactor 0.2 0.01 1.0 1e5 1e-2 r r r r r MixingParameter Mixing parameter. InverseJacobiWeight Weight for the difference of the inverse Jacobians (ω0 ). MinimalWeight Minimal allowed value for the weighting factors ωm . MaximalWeight Maximal allowed value for ωm . WeightFactor Weighting factor c for the calculation of the weighting factors ωm in (2.1). Note: As the Broyden-mixer stores a copy of the mixed quantity for each SCC iteration at a given geometry, you may consider to choose a different mixer with lower memory requirements, if your system needs density matrix mixing (e.g. DFTB+U), is large and needs a high number of SCCiterations (MaxSCCIteration). Anderson{} Modified Anderson mixer [14]. MixingParameter Generations InitMixingParameter DynMixingParameters DiagonalRescaling r i r (2r)* r 0.05 4 0.01 {} 0.01 MixingParameter Mixing parameter. Generations Number of generations to consider for the mixing. Setting it too high can lead to linearly dependent sets of equation. InitMixingParameter Simple mixing parameter used until the number of iterations is greater or equal to the number of generations. DynMixingParameters Allows specification of different mixing parameters for different levels of convergence during the calculation. These are given as a list of tolerances and the mixing factor to be used below this level of convergence. If the loosest specified tolerance is reached, the appropriate mixing parameter supersedes that specified in MixingParameter. DiagonalRescaling Used to increase the diagonal elements in the system of equations solved by the mixer. This can help to prevent linear dependencies occurring during the mixing process. Setting it to too large a value can prevent convergence. (This factor is defined in a slightly different way from Ref. [14]. See the source code for more details.) Example: 2.4. HAMILTONIAN 31 Mixer = Anderson { MixingParameter = 0.05 Generations = 4 # Now the over-ride the (previously hidden) default old settings InitMixingParameter = 0.01 DynMixingParameters = { 1.0e-2 0.1 # use 0.1 as mixing if more converged that 1.0e-2 1.0e-3 0.3 # again, but 1.0e-3 1.0e-4 0.5 # and the same } DiagonalRescaling = 0.01 } DIIS{} Direct inversion of the iterative space is a general method to acceleration iterative sequences. The current implementation accelerates the simple mix process. InitMixingParameter Generations UseFromStart r i l 0.2 6 Yes MixingParameter Mixing parameter. Generations Number of generations to consider for the mixing. UseFromStart Specifies if DIIS mixing should be done right from the start, or only after the number of SCC-cycles is greater or equal to the number of generations. Simple{} Constructs a linear combination of the current input and output charges as (1 − x)qin + xqout . MixingParameter r 0.05 MixingParameter Coefficient used in the linear combination. 2.4.2 SpinPolarisation In an SCC calculation, the code currently supports three different choices for spin polarisation; non-SCC calculations are not spin polarised. No spin polarisation ({}) No spin polarisation contributions to the energy or band-structure. Colinear{} Colinear spin polarisation in the z direction. The following options can be specified: 32 CHAPTER 2. INPUT FOR DFTB+ UnpairedElectrons RelaxTotalSpin InitialSpins r l p 0 No {} UnpairedElectrons Number of unpaired electrons. This is kept constant during the run, unless the RelaxTotalSpin keywords says otherwise. RelaxTotalSpin If set to Yes, a common Fermi-level is used for both spin channels, so that the total spin polarisation can change during run. In this case, the spin polarisation specified using the UnpairedElectrons keyword is only applied at initialisation. If set to No (default), the initial spin polarisation is kept constant during the entire run. InitialSpins Optional initialisation for spin patterns. If this keyword is present, the default code behaviour is that the initial input charge distribution has a magnetisation of 0. Otherwise if it is not present, the initial input charge distribution has a magnetisation matching the number of UnpairedElectrons keyword. The initial spin distribution for the input charges can be set by specifying the spin polarisation of atoms. For atoms without an explicit specification, a spin polarisation of zero is assumed. The InitialSpins property block must contain either the AllAtomSpins keyword with a list of spin polarisation values for every atom, or one or more AtomSpin blocks giving the spin for a specific group of atoms using the following keywords: Atoms SpinPerAtom (i|s)+ r - Atoms Atoms to specify an initial spin value. The atoms can be specified via indices, index ranges and species. (See MovedAtoms in section 2.3.1.) SpinPerAtom Initial spin polarisation for each atom in this InitialSpins block. For atoms not appearing in any of the SpinPerAtom block, an initial spin polarisation of 0 is set. Example (individual spin setting): SpinPolarisation = Colinear { UnpairedElectrons = 0.0 InitialSpins = { AtomSpin = { Atoms = 1:2 SpinPerAtom = -1.0 } AtomSpin = { Atoms = 3:4 SpinPerAtom = +1.0 } } } Example (setting all spins together): SpinPolarisation = Colinear { UnpairedElectrons = 0.0 InitialSpins = { 2.4. HAMILTONIAN 33 AllAtomSpins = { -1.0 -1.0 1.0 1.0 } # Atoms 1,2: -1.0, atoms 3,4: 1.0 } } NonColinear{} Non-collinear spin polarisation with arbitrary spin polarisation vector on every atom. The only option allowed is to set the initial spin distribution: InitialSpins {} p InitialSpins Initialisation of the x, y and z components of the initial spins for atoms. The default code behaviour is an initial spin polarisation of (0 0 0) for every atom. The initial spin distribution can be set by specifying the spin polarisation vector on all atoms using the AllAtomSpins keyword or by using one or more AtomSpin blocks with the following options: Atoms SpinPerAtom (i|s)+ (3r)+ - Atoms Atoms to specify an initial spin vector. The atoms can be specified via indices, index ranges and species. (See MovedAtoms in section 2.3.1.) SpinPerAtom Initial spin polarisation for each atom in this InitialSpins block. For atoms not appearing in any of the SpinPerAtom block, the vector (0 0 0) is set. Please note, that in contrast to the collinear case, in the non-collinear case you must specify the spin vector (3 components!) for the atoms. Example: SpinPolarisation = NonColinear { InitialSpins = { # Setting the spin for all atoms in the system AllAtomSpins = { 0.408 -0.408 0.816 0.408 -0.408 0.816 -0.408 0.408 -0.816 -0.408 0.408 -0.816 } } } Example: SpinPolarisation = NonColinear { InitialSpins = { AtomSpin = { Atoms = 1:2 SpinPerAtom = 0.408 -0.408 0.816 } AtomSpin = { Atoms = 3:4 34 CHAPTER 2. INPUT FOR DFTB+ SpinPerAtom = -0.408 0.408 -0.816 } } } SpinConstants This environment suplies the atomic constants required for either spin polarised calculations or when evaluating properties which depend on spin interactions (triplet excitations for example). In these cases, for each atomic species in the calculation the spin coupling constants for that atom must be specified. When ShellResolvedSCC = No, an extra keyword in this block controls whether the spin constants are resolved by shell or are identical for all shells: ShellResolvedSpin, defaulting to the same value as ShellResolvedSCC. When shell resolved spin constants are specified, they must be ordered with respect to the pairs of shells they couple, such that the index for the second shell increases faster. For an spd-basis, that gives the following ordering: wss , wsp , wsd , . . . , w ps , w pp , w pd , . . . , wds , wd p , wdd , . . . Example (GGA parameters for H2 O): SpinConstants = { O={ # Wss Wsp Wps Wpp -0.035 -0.030 -0.030 -0.028 } H={ # Wss -0.072 } } Several standard values of atomic spin constants are given in appendix D. Constants calculated with the same density functional as the SK-files should be used. This input block may be moved to the SK-data definition files in the future. When using the SelectedShells method for the keyword MaxAngularMomentum, the spin constants are listed as an array of values running over SK1SK2 . . . in the same order as listed for SlaterKosterFiles. SpinConstants = { # not real values, only an example Si = { # Wss Wsp Wss* -0.035 -0.030 -0.01 # Wps Wpp Wps* -0.030 -0.037 -0.02 # Ws*s Ws*p Ws*s* -0.01 -0.02 -0.01 } 2.4. HAMILTONIAN 35 For cases where ShellResolvedSpin = No, the spin constant for the the highest occupied orbital of each atom should be supplied: Example (GGA parameters for H2 O): SpinConstants = { O={ #Wpp -0.028 } H={ # Wss -0.072 } } 2.4.3 SpinOrbit If present, specifies that L · S coupling should be included for a calculation. Currently spin unpolarised and non-collinear spin polarisation are supported, but not collinear spin polarisation. For every atomic species present in the calculation the spin-orbit coupling constants, ξ , for that atom must be specified for all shells present. The constants must be ordered with respect to the list of shells for the given atom. In the case that the spin-orbit constant for an s orbital has been set to be a non-zero value the code prints a warning. For periodic systems, use of this keyword automatically prevents the folding by inversion normally used in SupercellFolding{}, but manually specified KPointsAndWeights should not be reduced by inversion. Example (GaAs): SpinOrbit = { Ga [eV] = {0.0 0.12 0.0} # s p d shells As [eV] = {0.0 0.32703} # s p shells } The additional option in this block, Dual, sets whether to use a block population for the local spin matrices consistent with the dual populations of Han et al. [20] or the conventional on-site part of the single particle density matrix. The default value of this option is Yes, also giving extra information regarding atomic orbital moments in the detailed output. 2.4.4 Solver Currently the following LAPACK 3.0 [4] eigensolver methods are always available: • QR{} (QR decomposition based solver) • DivideAndConquer{} (this requires about twice the memory of the other solvers) • RelativelyRobust{} (using the subspace form but calculating all states) 36 CHAPTER 2. INPUT FOR DFTB+ • MAGMA{} (Only available for DFTB+ binaries compiled with MAGMA [51, 52, 11] GPU support. WARNING: this is currently an experimental feature, so should be used with care.) None of these solvers need any parameters or properties to be specified. Example: Solver = DivideAndConquer {} For ScaLAPACK enabled compilation, all three solvers are also available for MPI parallel use. If DFTB+ is compiled with the ELSI library also included [54], the additional ELPA, OMM, PEXSI and NTPoly solvers also become available. Note: The ELSI-solvers are not tested with multiple OpenMP-threads. Therefore, DFTB+ will stop with an error, if an ELSI-solver has been selected and the maximal number of allowed threads is greater than one. (You can control the number of allowed OpenMP-threads via the OMP_NUM_THREADS environment variable.) ELPA This is available with either single or two stage solution methods (the second of these should be more efficiently parallel for large problems). Example: Solver = ELPA { Mode = 2 } One caveat for this solver is that the number of parallel groups (see p. 62) must match the number of k-points (times 2 in the case of collinear spin polarisation). Calculations without k-points can use either one or two groups in the case of collinear spin polarisation. OMM This method minimises the single particle density matrix, so does not make band structure information available. It is only stable for insulating grounds states, i.e., systems with a HOMO-LUMO (band) gap. The orbital minimisation method has four options: nInterationsELPA Tolerance Choleskii Sparse i r l l 5 1E-10 Yes No nInterationsELPA Number of initial iterations to be performed with ELPA before the OMM method starts. Tolerance Minimisation tolerance for this solver, larger values are faster by may be less stable. 2.4. HAMILTONIAN 37 Choleskii Whether the overlap is Choleskii factorised before applying OMM. This may increase stability of this method. Sparse Whether the code should use the sparse matrix interface to ELSI solvers. This does not substantially improve memory usage in this case as internally the dense problem is solved with libOMM. PEXSI The PEXSI solver directly calculates the density matrix, so does not make band structure information or Mermi free energy available. The scaling with system size is better than the other solvers d/2+1/2 available in DFTB+, increasing as O(Natom ) where d is the effective dimensionality of the system. Hence for three dimensional structures it will scale as O(N 2 ) for general systems. Poles ProcsPerPole muPoints SymbolicFactorProcs SpectralRadius Sparse 20 1 2 1 10 No i i i i r l Poles number of poles for the complex plane calculation. ProcsPerPole processors used to calculate the inversion at each pole. muPoints number of processors used to search for the Fermi level. SymbolicFactorProcs number of processors to use in evaluating the factorisation pattern of matrices. SpectralRadius [Energy] extension of the complex contour. Sparse Whether the code should use the sparse ELSI matrix interface. NTPoly This method constructs the single particle density matrix via a purification method based on matrix polynomials (hence requires insulating systems). The solver does not make band structure information available, but can be linear scaling in both time and memory depending on settings and system. Currently the solver does not support spin polarisation or k-points. This solver has four options: PurificationMethod Tolerance Truncation Sparse i r r l 2 1E-5 1E-10 No PurificationMethod Allowed choices are 0 for canonical purification, 1 for trace correcting purification, 2 for 4th order trace resetting purification, and 3 for generalised hole-particle canonical purification. 38 CHAPTER 2. INPUT FOR DFTB+ Tolerance Iterative convergence tolerance for this solver, larger values are faster by may be less stable. Truncation Tolerance below which matrix elements in the density matrix are dropped to enforce sparsity. Sparse Whether the code should use the sparse matrix ELSI interface. The default choices of Tolerance and Truncation lead to an accurate, but slow, solutions. Alternatively linear scaling can be achieved at smaller system sizes with a larger choice of these values. Values in the range of 1E-3 and 1E-6 for Tolerance and Truncation may be suitable (but test the quality of the solutions). 2.4.5 Filling There are currently two types of filling supported (see below). Both have common options: Temperature IndependentKFilling FixedFermiLevel r l (1|2)r AdaptFillingTemp = No Periodic = Yes 0.0 No - Temperature [energy] Electron temperature in energy units. This property is ignored for thermostated MD runs, if the AdaptFillingTemp property of the thermostat has been set to Yes (See p. 17). IndependentKFilling Causes the occupation of the eigenstates to be independently determined for each k-point, thus preventing electron transfer between the k-points. Please note that the value for the Fermi level printed out by the code is meaningless in that case, since there is no common Fermi level for all k-points. This option is incompatible with use of the FixedFermiLevel keyword. FixedFermiLevel [energy] Can be used to fix the Fermi-level (total chemical potential, µ) of the electrons in the system. For collinear spin polarisation, values for up and down spin channels are required. Otherwise only a single global chemical potential is required. If this option is present, the total charge and the total spin of the system are not conserved (settings in the options Charge and UnpairedElectrons will be ignored). If a fixed chemical potential is used, the output force related energy includes the contribution to the free energy, −Nµ, hence if differentiated will give the forces and stresses (if periodic). Fermi{} Fills the single particle levels according to a Fermi distribution. When using a finite temperature, the Mermin free energy (which the code prints) should be used instead of the total energy. This is given by E − T S, where the electron entropy S is used. Example: Filling = Fermi { Temperature [K] = 300 } 2.4. HAMILTONIAN 39 MethfesselPaxton{} Produces a Fermi-like distribution but with much lower electron entropy [35]. This is useful for systems that require high electron temperatures (for example when calculating metallic systems). There is an additional option for this type of filling: Order 2 i Order Order of the Methessel-Paxton scheme, the order must be greater than zero, and the 1st order scheme is equivalent to Gaussian filling. Note: Due to the non-monotonic behaviour of the Methfessel-Paxton filling function, the position of the Fermi-level is not necessary unique for a given number of electrons. Therefore, different fillings, band entropies, and Mermin free energies may result, depending which one has been found by the Fermi-level search algorithm. The differences, however, are usually not physically significant. 2.4.6 SlaterKosterFiles There are two different ways to specify the Slater-Koster files for the atom type pairs, explicit specification and using the Type2FileNames{} method. Explicit specification Every pairwise permutation atomic types, connected by a dash, must occur as a property with the name of the corresponding file as an assigned value. Example (GaAs): SlaterKosterFiles = { Ga-Ga = "./Ga-Ga.skf" Ga-As = "./Ga-As.skf" As-Ga = "./As-Ga.skf" As-As = "./As-As.skf" } If you treat shells from different species as shells of one atom by using the SelectedShells{} keyword in the MaxAngularMomentum{} block, you have to specify more than one file name for certain species pairs. (For details see the description about the MaxAngularMomentum{} keyword.) Type2FileNames{} You can use this method to generate the name of the Slater-Koster files automatically using the element names from the input geometry. You have to specify the following properties: Prefix Separator Suffix LowerCaseTypeName s s s l Prefix Prefix before the first type name, usually the path. "" "" "" No 40 CHAPTER 2. INPUT FOR DFTB+ Separator Separator between the type names. Suffix Suffix after the name of the second type, usually extension. LowerCaseTypeName If the name of the types should be converted to lower case. Otherwise they are used in the same way, as they were specified in the geometry input. Example (for producing the same file names as in the previous section): SlaterKosterFiles = Type2FileNames { Prefix = "./" Separator = "-" Suffix = ".skf" LowerCaseTypeName = No } The Type2FileNames method can not be used if an extended basis was defined with the SelectedShells method. 2.4.7 KPointsAndWeights The k-points for the Brillouin-zone integration can either be specified explicitly or using the KLines{} or the SupercellFolding{} methods. If the latter is used the KPointsAndWeights keyword is not allowed to have a modifier. Explicit specification For every k-point four real numbers must be specified: The coordinates of the given k-point followed by its weight. By default, the coordinates are specified in fractions of the reciprocal lattice vectors. If the modifier absolute is set for the KPointsAndWeights keyword, absolute k-point coordinates in atomic units are instead expected. The sum of the k-point weights is automatically normalised by the program. KPointsAndWeights = { # 2x2x2 MP-scheme 0.25 0.25 0.25 1.0 0.25 0.25 -0.25 1.0 0.25 -0.25 0.25 1.0 0.25 -0.25 -0.25 1.0 } SupercellFolding{} This method generates a sampling set containing all the special k-points in the Brillouin zone related to points that would occur in an enlarged supercell repeating of the current unit cell. If two k-points in the BZ are related by inversion, only one (with double weight) is used (in the absence of spin-orbit coupling this is permitted by time reversal symmetry). The SupercellFolding{} method expects 9 2.4. HAMILTONIAN 41 integers and 3 real values as parameters: n11 n12 n13 n21 n22 n23 n31 n32 n33 s1 s2 s3 The integers ni j specify the coefficients used to build the supercell vectors Ai from the original lattice vectors a j : 3 Ai = ∑ ni j a j . j=1 The real values, si , specify the point in the Brillouin-zone of the super lattice, in which the folding should occur. The coordinates must be given in relative coordinates, in the units of the reciprocal lattice vectors of the super lattice. The original l1 × l2 × l3 Monkhorst-Pack sampling [37] for cubic lattices corresponds to a uniform extension of the lattice: l1 0 0 0 l2 0 0 0 l3 s1 s2 s3 where si is 0.0, if li is odd, and si is 0.5 if li is even. For the 2 × 2 × 3 scheme, you would write for example # 2x2x3 MP-scheme according original paper KPointsAndWeights = SupercellFolding { 2 0 0 0 2 0 0 0 3 0.5 0.5 0.0 } To use k-points for hexagonal lattices which are consistent with the erratum to the original paper [38], you should set the shift for the unique “c” direction, s3 , in the same way as in the original scheme. The s1 and s2 shifts should be set to be 0.0 independent of whether l1 and l2 are even or odd. So, for a 2 × 3 × 4 sampling you would have to set # 2x3x4 MP-scheme according modified MP scheme KPointsAndWeights = SupercellFolding { 2 0 0 0 3 0 0 0 4 0.0 0.0 0.5 } It is important to note that DFTB+ does not take the symmetry of your system explicitly into account. For small high symmetric systems with a low number of k-points in the sampling this could eventually lead to unphysical results. (Components of tensor properties–e.g. forces–could be finite, even if they must vanish due to symmetry reasons.) For those cases, you should explicitly specify k-points with the correct symmetry. 42 CHAPTER 2. INPUT FOR DFTB+ KLines{} This method specifies k-points lying along arbitrary lines in the Brillouin zone. This is useful when calculating the band structure for a periodic system. (In that case, the charges should be initialised from the saved charges of a previous calculation with a proper k-sampling. Additionally for SCC calculations the number of SCC cycles should be set to 1, so that only one diagonalisation is done using the initial charges.) The KLines{} method accepts for each line an integer specifying the number of points along the line segment, and 3 real values specifying the end point of the line segment. The line segments do not include their starting points but their end points. The starting point for the first line segment can be set by specifying a (zeroth) segment with only one point and with the desired starting point as end point. The unit of the k-points is determined by any modifier of the KPointsAndWeights property. (Default is relative coordinates.) Example: KPointsAndWeights [relative] = KLines { 1 0.5 0.0 0.0 # Setting (and calculating) starting point 0.5 0.0 0.0 10 0.0 0.0 0.0 # 10 points from 0.5 0.0 0.0 to 0.0 0.0 0.0 10 0.5 0.5 0.5 # 10 points from 0.0 0.0 0.0 to 0.5 0.5 0.5 1 0.0 0.0 0.0 # Setting (and calculating) a new starting point 10 0.5 0.5 0.0 # 10 points from 0.0 0.0 0.0 to 0.5 0.5 0.0 } Note: Since this set of k-points probably does not correctly integrate the Brillouin zone, the default value of MaxSccIterations is set to be 1 in this case. 2.4.8 OrbitalPotential Currently the FLL (fully localised limit) and pSIC [24] (pseudo self interaction correction ) forms of the LDA+U corrections [44] are implemented. These potentials effect the energy of states on designated shells of particular atoms, usually increasing the localisation of states at these sites. The FLL potential lowers the energy of occupied states localised on the specified atomic shells while raising the energy of unoccupied states. The the pSIC potential corrects the local part of the selfinteraction error and so lowers the energy of occupied states (see Ref. [24] for a discussion of the relation between these two potentials, and possible choices for the UJ constant). These particular corrections are most useful for lanthanide/actinide f states and some localised d states of transition metals (Ni3d for example). The Functional option chooses which correction to apply, followed by a list of the specific corrections, listed as an atomic species and the shells on that atom followed by the U − J constant for that block of shells. OrbitalPotential = { Functional = {FLL} Si = { Shells = {1 2} # sp block on the atom UJ = 0.124 } } 2.4. HAMILTONIAN 2.4.9 43 ElectricField This tag contains the specification for an external electric field. Electric fields can only be specified for SCC calculations. You can apply the electric field of point charges3 and/or a homogeneous external field (which may change harmonically in time). The ElectricField block can currently contain either one or more PointCharges blocks and potentially an External block. PointCharges The specification for PointCharges has the following properties: CoordsAndCharges GaussianBlurWidth (4r)+ r Periodic = No 0.0 CoordsAndCharges [length] Contains the coordinates and the charge for each point charge (four real values per point charge). A length modifier can be used to alter the units of the coordinates. The charge must be specified in proton charges. (The charge of an electron is -1.) If you read in a huge number of external charges the parsing time to process this data could be unreasonably long. You can avoid this by including the coordinates and the charges directly from an external file via the DirectRead{} method (see the example in the next paragraph). Please note that when using this method the program will only read the specified number of records from the external file, and ignores any additional data (so do not leave comments in the external file for example). The external file should contain only one record (3 coordinates and 1 charge) per line. GaussianBlurWidth [length] Specifies the half width σ of the Gaussian charge distribution, which is used to delocalise the point charges. The energy of the coulombic interaction EC between the delocalised point charge M with charge QM and the atom A with charge qA is weighted by the error function as hr i qA QM AM EC (A, M) = erf , rAM σ where rAM is the distance between the point charge and the atom. This delocalisation can only be used for non-periodic systems. A length modifier can be used to specify the unit for σ . Example: ElectricField = { # 1st group of charges, with Gaussian delocalisation # We have 100000 charges, therefore we choose the fast reading method. PointCharges = { GaussianBlurWidth [Angstrom] = 3.0 CoordsAndCharges [Angstrom] = DirectRead { Records = 100000 File = "charges.dat" } } # 2nd group of charges, no delocalisation (sigma = 0.0) PointCharges = { 3 Only in calculations with fixed lattice constants. 44 CHAPTER 2. INPUT FOR DFTB+ CoordsAndCharges [Angstrom] = { 3.3 -1.2 0.9 9.2 1.2 -3.4 5.6 -3.3 } } } External Specifies a homogeneous external electric field. In the case of periodic calculations, a saw-tooth potential is currently used, hence it is up to the user to guarantee that there is a vacuum region isolating periodic copies of the system along the applied field direction. We suggest that you place the structure in the ‘middle’ of the unit cell if possible, to reduce the chances of atoms approaching cell boundaries along the direction of the applied electric field. The code will halt if atoms interact with periodic images of the unit cell along the direction of the electric field. The External field keyword has the following options Strength Direction Frequency Phase r 3r r i molecular dynamics used Geometry step offset 0.0 0 Strength [Electric field strength] Specified strength of the applied field. Direction Vector direction of the applied field (the code normalises this vector). In the case of periodic calculations, currently the system must not be continuous in this direction (see above). Frequency [Frequency] If using molecular dynamics, the field can be time varying with this frequency. Phase Initial field phase in units of geometry steps, this is needed if restarting an MD run in an external field to give the offset in phase of the field after the specified number of steps from the old calculation. The applied field is of the form E0 sin(ω∆t(step + phase)) where E0 is the field vector specified by Strength and Direction, ω is the angular Frequency and step is the current MD-step in the simulation, using the MD TimeStep of ∆t (see section 2.3.6). 2.4.10 Dispersion The Dispersion block controls whether DFTB interactions should be empirically corrected for van der Waals interactions, since DFTB (and SCC-DFTB) does not include these effects. Currently, three different dispersion correction schemes are implemented (for the detailed description of the methods see the following subsections): • LennardJones: Dispersion is included via a Lennard-Jones potential between each pair of atoms. The parameters for the potential can either be entered by the user or the program can automatically take the parameters from the Universal Force Field (UFF) [47]. 2.4. HAMILTONIAN 45 • SlaterKirkwood: The dispersion interaction between atoms is taken from a Slater-Kirkwood polarisable atomic model [12]. • DftD3: Dispersion is calculated as in the dftd3 code [18, 19] (see section 2.4.10). Modification hydrogen bond interaction strengths (see section 2.4.12). LennardJones The Lennard-Jones dispersion model in DFTB+ follows the method of Ref. [55], using the following potential:       ri j 6 ri j 12 Ui j (r) = di j −2 + r r Ui j (r) = U0 +U1 r5 +U2 r10 r >= r0 r < r0 where r0 is the distance at which the potential turns from repulsive to attractive. The parameters p di j and ri j are built from atomic parameters di , d j and ri , r j via the geometrical mean (di j = di d j , √ ri j = ri r j ). The parameters U0 , U1 , U2 ensure a smooth functional form at r0 . The parameters ri and di can either be taken from the parameters of the UFF [47] (as in Ref. [55]) or can be specified manually for each species. Example using UFF parameters: Dispersion = LennardJones { Parameters = UFFParameters {} } Example using manually specified parameters: Dispersion = LennardJones { Parameters { H{ Distance [AA] = 2.886 Energy [kcal/mol] = 0.044 } O{ Distance [AA] = 3.500 Energy [kcal/mol] = 0.060 } } } The UFF provides dispersion parameters for nearly every element of the periodic table, therefore it can be used for almost all systems “out of the box”. The parameters are also independent of the atomic coordination number, allowing straight forward geometry relaxation or molecular dynamics (unlike the current implementation of Slater-Kirkwood type dispersion). 46 CHAPTER 2. INPUT FOR DFTB+ SlaterKirkwood A Slater-Kirkwood type dispersion model is also implemented in DFTB+ as described in Ref. [12].4 This model requires atomic polarisation values, van der Waals radii and effective charges for every atom in your system. These parameters are dependent on the coordination of each atom, hence values for different atoms of the same species may vary depending on local environment. You can supply these parameters for the atoms in either of two ways, both using the PolarRadiusCharge tag. The first option is to specify the values within the PolarRadiusCharge environment by providing three real values (polarisability, van der Waals radius, effective charge) for each atom separately. Example: Dispersion = SlaterKirkwood { # Using Angstrom^3 for volume, Angstrom for length and default # unit for charge (note the two separating commas between the units) PolarRadiusCharge [Angstrom^3,Angstrom,] = { # Polar Radius Chrg 0.560000 3.800000 3.150000 # Atom 1: O 0.386000 3.500000 0.800000 # Atom 2: H 0.386000 3.500000 0.800000 # Atom 3: H } } Alternatively you can provide values for each atomic species in your system, but must supply different values for different coordination numbers using the HybridDependentPol{} keyword. The code needs specific parameters for each type of atom in environments with 0, 1, 2, 3, 4 or >5 neighbours. DFTB+ then picks the appropriate values for each atom based on their coordination in the starting geometry. Two atoms are considered to be neighbours if their distance is less than the sum of their covalent radii, hence you need to supply the covalent radii for each atomic species using the CovalentRadius keyword. This is then followed by a HybridPolarisations block containing a list of six values for atomic polarisabilities then six van der Waals radii and finally a single hybridisation independent effective charge for that atomic species. Example: Dispersion = SlaterKirkwood { PolarRadiusCharge = HybridDependentPol { O={ CovalentRadius [Angstrom] = 0.8 HybridPolarisations [Angstrom^3,Angstrom,] = { # Atomic polarisabilities 0-5 van der Waals radii 0-5 chrg 0.560 0.560 0.560 0.560 0.560 0.560 3.8 3.8 3.8 3.8 3.8 3.8 3.15 } } H={ CovalentRadius [Angstrom] = 0.4 HybridPolarisations [Angstrom^3,Angstrom,] = { β α 4 Please note, that Ref. [12] contains two typos: equation (7) should read Cαβ = 2C6 C6 pα pβ , in equation (9) the conβ 6 p2 C +p2 Cα α 6 β 6 tribution from the dispersion should be Edis = − 21 ∑αβ f (Rαβ )C6 (Rαβ )−6 . This option is also currently incompatible with lattice optimisation and the use of barostats. αβ 2.4. HAMILTONIAN 47 # Atomic polarisabilities 0-5 van der Waals radii 0-5 chrg 0.386 0.396 0.400 0.410 0.410 0.410 3.5 3.5 3.5 3.5 3.5 3.5 0.8 } } } } Warning: For both methods of specifying the Slater-Kirkwood dispersion model the code keeps the dispersion parameters fixed for each atom during the entire calculation. Even if the geometry (and therefore the hybridisation) of atoms changes significantly during the calculation, the parameters are unchanged. Therefore if atoms are able to move during your calculation (geometry relaxation or molecular dynamics) you should always check whether the coordination of your atoms has changed during the run. Furthermore, when using the HybridDependentPol{} method we suggest that you first set the StopAfterParsing keyword in the ParserOptions block to Yes (see p. 61) and inspect the generated polarisabilities, radii and charges for every atom in the dftb_pin.hsd file. If fine tuning of the generated values turns out to be necessary, you should replace the hybrid dependent specification in the input file with corrected atom specific values based on dftb_pin.hsd. In order to find suitable parameters for the Slater-Kirkwood model, you should consult Ref. [12] and further references therein. Appendix E contains values which have already been used by some DFTB-users for a few elements. DftD3 The DFT-D3 dispersion correction in DFTB+ is an implementation of the method used in the code ’dftd3’ by Stefan Grimme and coworkers. It is based on the ansatz described in Refs. [18] and [19]. Note: the DFTB+ binary must be compiled with the DFT-D3 library enabled to use this feature. This dispersion correction for DFTB adds a contribution to the general Kohn-Sham-like energy EDFTB-D3 = EDFTB + Edisp with EDFTB being the DFTB total energy and Edisp the dispersion energy. The latter contains twobody and optional three-body contributions: (2) (3) Edisp = Edisp + Edisp The form of the two-body contribution can change depending on the chosen damping factor: • Becke-Johnson damping function: (2) Edisp = − 1 CnAB s ∑ nn AB 2 A6∑ =B n=6,8 rAB + f (R0 ) with AB f (RAB 0 ) = a1 R0 + a2 . • Zero-damping (dispersion at short distances is damped to zero): (2) Edisp = − 1 CnAB s n n f d,n (rAB ) 2 A6∑ rAB =B 48 CHAPTER 2. INPUT FOR DFTB+ with fd,n = 1 −αn 1 + 6(rAB /(sr,n RAB 0 )) In order to adjust the dispersion for various energy functionals, the choice of s6 , s8 and the damping parameters a1 and a2 (for Becke-Johnson-damping) or sr,6 and α6 (for zero damping) are treated as functional-dependent values. All other parameters are fixed based on these parameters. As the DFTB energy functional is largely determined by the underlying parameterisation (the SlaterKoster-files) and the chosen DFTB model (e.g. non-scc, scc, 3rd order, etc.), there are no universal parameter choices which can be used with all settings, but some relevant choices for various parameterisation are given in Appendix F. Note: for the version 6 or earlier of the DFTB+ input parser (see section 2.8) the default values of these parameters are set to be appropriate for DFTB3. But from parser version 7 onwards, no default values are set. Example using adjusted parameters with Becke-Johnson damping: Dispersion = DftD3 { Damping = BeckeJohnson { a1 = 0.5719 a2 = 3.6017 } s6 = 1.0 s8 = 0.5883 } Example using zero-damping: Dispersion = DftD3 { Damping = ZeroDamping { sr6 = 0.7461 alpha6 = 14.0 } s6 = 1.0 s8 = 3.209 } DftD3 optional settings Apart from the functional dependent dispersion parameters, you can also adjust the additional parameters as shown below. The default values for these parameters are taken to be the same as in the dftd3 code. √ Cutoff r 9000 CutoffCN r 40 Threebody l No HHRepulsion l No Cutoff [length] Cutoff distance when calculating two-body interactions. CutoffCN [length] Cutoff distance when calculating three-body interactions. 2.4. HAMILTONIAN 49 Threebody Whether three-body contributions should be included in the dispersion interactions. HHRepulsion Required when calculating the DFTB3-D3H5 [48] modification to D3 dispersion (see section 2.4.12 for details and parameter values). This keyword enables an additional short range repulsion term in all hydrogen–hydrogen pairs [49] which prevents them from approaching too closely together. 2.4.11 DFTB3 If you would like to use what is called “DFTB3” in some publication(s) [17], this group of options include the relevant modifications to the SCC Hamiltonian and energy. To enable the DFTB3 model you will need to set ThirdOrderFull = Yes and damp H–X the interactions (see Section 2.4.12). ThirdOrder If set to Yes the on-site 3rd order correction [53] is switched on. This corrects the SCC-Hamiltonian with the derivatives of the Hubbard U parameters, which you have to specify for every element in HubbardDerivs. This correction only alters the on-site elements and is only maintained for backward compatibility. You should use the full version ThirdOrderFull instead. ThirdOrderFull If set to Yes the full 3rd order correction [17] is switched on. This corrects the SCC-Hamiltonian with the derivatives of the Hubbard U parameters, which you have to specify for every element in HubbardDerivs. HubbardDerivs Derivatives of the Hubbard U for the 3rd order correction (on-site or full). For every element the appropriate parameter (in atomic units) must be specified. If you use shell resolved SCC (with full 3rd order), you must specify a list of derivatives for every element, with one Hubbard U derivative for each shell of the given element. Hamiltonian = DFTB { : ThirdOrder = Yes HubbardDerivs { O = -0.14 H = -0.07 } : } 2.4.12 Hydrogen corrections There are currently two available methods to correct hydrogen interactions (mainly hydrogen bonds) in the HCorrection environment: Damping The Damping method modifies the short range contribution to the SCC interaction between atoms A and B with the damping factor  ζ U +U 2 − Al 2 Bl rAB e 50 CHAPTER 2. INPUT FOR DFTB+ provided that at least one of the two atoms is hydrogen [17, 53]. (UAl and UBl are the Hubbard U values of the two atoms for the l-shell, rAB is the distance between the atoms.) An atom is considered to be a hydrogen-like atom, if its mass (stored in the appropriate homonuclear SK-file) is less than 3.5 amu. The Exponent keyword in this environment sets the parameter ζ for the short range damping: HCorrection = Damping { Exponent = 4.05 } Table 2 of reference [17] gives suggested values of the exponent for different DFTB2 and DFTB3 models applied to light atoms bonded to hydrogen. DFTB3-D3H5 DFTB3-D3H5 [48] is a variant of DFTB3 with additional corrections for non-covalent interactions (dispersion and hydrogen bonds). It consists of a third-order DFTB calculation using the 3OB parameter set, but where the gamma-function damping (Damping method above) is replaced by the H5 correction and an additional D3 dispersion correction in included. This method also includes a repulsive term which is added to prevent unphysically close approach of pairs of hydrogen atoms [49]. Setting the HCorrection environment to H5{} activates this correction for hydrogen bonds [48]. If no additional parameters are provided in the input, suitable values for H-{O,N,S} systems are used (the correction was developed for the DFTB3/3OB model and parameters). HCorrection = H5 {} Note: It was found that DFTB3 overestimates the strength of H-bonds involving the terminal nitrogen of an azide group, and the published results in Ref. [48] were obtained with the H5 correction switched off for these specific atoms. To reproduce this behavior in a system containing nitrogen in several environments, a new atom type with a different name but the same DFTB parameters can be used for specific N atoms to which the correction should not be applied. If you want to specify the parameters manually, H5 accepts following options, corresponding to terms in Ref. [48]: RScaling WScaling H5Scaling r r m 0.714 0.25 RScaling Global scaling factor, sr , when calculating the position of the correcting gaussian functions: r0 = sr (rvdW (X) + rvdW (H)) . WScaling Global scaling factor, sW , when calculating the width of the correcting gaussian functions. The full-width at at half-maximum of the gaussian, w, is normalised to be 1 for a unit value of WScaling: sw (rvdW (X) + rvdW (H)) √ . w= 2 2 ln 2 2.4. HAMILTONIAN 51 H5Scaling Atom type specific scaling pre-factor, kXH , of the correcting gaussian functions when calculating the SCC-interaction: H5 γXH = γXH (rXH − r0 )2 1 + kXH exp − 2w2 !! . You will have to specify one value for each of the chemical species you would like to correct (see the example below). Explicitly setting a negative value (e.g. -1.0) for a given atom type switches off the correction for hydrogen bonds involving that type of atom. In the special cases of N, O or S, if you do not specify a value (and do not disable the contribution by using -1.0), the default value from the reference paper will be used [48]. For any other omitted atom types, the code defaults to a choice of -1.0 (no correction). Hamiltonian = DFTB { : HCorrection = H5 { RScaling = 0.714 WScaling = 0.25 H5Scaling { O = 0.06 N = 0.18 S = 0.21 } } : } Note: The van der Waals radii (rvdW ) of atoms are also required. DFTB+ stores these for most of the periodic table, but for cases that are not available their contribution to this correction are neglected. For a DFTB3-D3H5 calculation, a specific parametrization of the D3 dispersion has to be used. In addition to setting up appropriate values of the D3 parameters, as discussed in Ref. [48], the hydrogen–hydrogen repulsion of Ref. [49] has to also be activated. The complete input is: Hamiltonian = DFTB { : Dispersion = DftD3 { Damping = ZeroDamping { sr6 = 1.25 alpha6 = 29.61 } s6 = 1.0 s8 = 0.49 HHRepulsion = Yes } : } 52 CHAPTER 2. INPUT FOR DFTB+ 2.4.13 RangeSeparated The RangeSeparated keyword specifies the use of a range separated hybrid functional. Currently, only the long-range corrected hybrid functional (LC) [39, 33] is implemented. There, the electrostatic interaction is split up into long and short ranged components according to 1 1 − e−ωr e−ωr = + , r r r with the range-separation parameter ω, which is set in the Slater-Koster files. The option should only be used with corresponding parameter sets created for use with long-range correction. Note: The present release does not yet support long-range corrected excited states calculations with LC-TD-DFTB, and the RangeSeparated keyword can, therefore, not be used in conjunction with the ExcitedState keyword. The RangeSeparated keyword expects either None (default – no use of range-separated hybrid functional) or the LC{} block as value. Latter enables the following option: Screening m Thresholded {} Screening Choice of the screening method. The following choices are allowed: Thresholded {} Screening according to estimated magnitude of terms. This is the recommended choice for speed and accuracy but does not support all of the cases (restarting and spin polarisation). Threshold r 1e-6 CutoffReduction r 0.0 Threashold Threashold, below which elements are considered to be zero. CutoffReduction [length] Reduces the spatial cutoff, beyond which the overlap between atoms is considered to be zero. This can be used as an additional tweak to speed up the LC-calculation, but make sure first, that your results do not change considerably. Default: 0.0 – no reduction, using the cutoff from the SK-files. NeighbourBased Uses a purely neighbour-list based algorithm. This algorithm is usually considerably slower than the Thresholded. CutoffReduction r 0.0 CutoffReduction [length] See description in the Thresholded block. Example for thresholded screening with customised threshold value. RangeSeparated = LC { Screening = Thresholded { Threshold = 1e-5 } } Example for neighbour list based screening with customised cutoff reduction: 2.4. HAMILTONIAN 53 RangeSeparated = LC { Screening = NeighbourBased { CutoffReduction [AA] = 2.0 } } 2.4.14 On site corrections This block enables corrections for on-site matrix elements which improve the description of multicentre integrals [16] leading to, for example, improved hydrogen-bond energies [10]. For each chemical species, the spin-same-spin and spin-different-spin constants should be specified for all combinations of atomic shells. note: the matrix of constants is symmetric and the purely s-with-s entries are zero (the code ignores their value due to symmetry). Example: OnSiteCorrection= { # same spin oxygen Ouu = {0.00000 0.08672 0.08672 -0.00523} # hetero-spin oxygen Oud = {0.00000 0.14969 0.14969 0.03834} # H all zero Huu = {0} Hud = {0} } Some on-site constants are given in appendix G. 2.4.15 Differentiation Calculations of forces currently require the numerical derivatives of the overlap and non-selfconsistent Hamiltonian. This environment controls how these derivatives are evaluated. Note: In earlier DFTB+ versions (up to version 1.2), differentiation was done using finite difference derivatives with a step size of 0.01 atomic units. If you want to reproduce old results, choose the FiniteDiff method and set the step size explicitly to this value. FiniteDiff{} Finite difference derivatives with a specified step size Delta r Delta [length] Step size 1 epsilon /4 54 CHAPTER 2. INPUT FOR DFTB+ Richardson{} Extrapolation of finite difference via Richardson’s deferred approach to the limit (in principle the most accurate of the currently available choices). 2.4.16 ForceEvaluation Chooses the method for evaluating the electronic contribution to the forces. ’traditional’ Uses the “traditional” DFTB-force expression, given for example, in Ref. [13]. ’dynamics’ Force expression from Ref. [6]. This choice should be used if forces are being calculated with non-converged charges (e.g. when doing XLBOMD dynamics). Note: this force expression is only compatible with the Fermi filling (see keyword Filling, p. 38.) ’dynamicsT0’ Simplified dynamic force expression valid for electronic temperature T = 0 K [6]. This choice should be used if forces are calculated with non-converged charges and the electronic temperature is zero (e.g. when doing XLBOMD dynamics at T = 0 K). Note: that XLBOMD calculations (Section 2.3.6) are not able to use the ’traditional’ forces. Example: ForceEvaluation = ’dynamics’ 2.5 Options This block collects some global options for the run. WriteAutotestTag WriteDetailedXML WriteResultsTag WriteDetailedOut RestartFrequency RandomSeed MinimiseMemoryUsage TimingVerbosity ShowFoldedCoords WriteHS WriteRealHS ReadChargesAsText WriteChargesAsText SkipChargeTest l l l l i i l i l l l l l l Driver = {}, SCC = Yes Periodic = Yes ReadInitialCharges = Yes ReadInitialCharges = Yes No No No Yes 20 0 No 0 No No No No No No WriteAutotestTag Turns the creation of the autotest.tag file on and off. (This file can get quite big and is only needed for the autotesting framework.) WriteDetailedXML Turns the creation of the detailed.xml file on and off. (The detailed.xml file is needed among others by the waveplot utility for visualising molecular orbitals.) 2.5. OPTIONS 55 WriteResultsTag Turns the creation of the results.tag file on and off. (That file is used by several utilities processing the results of DFTB+.) For a description of the file format see p. 64. WriteDetailedOut Controls the creation of the file detailed.out (see p. 63). Since this contains the detailed information about the last step of your run, you shouldn’t turn it off without good reasons. RestartFrequency Specifies the interval at which charge restart information should be written to disc for static SCC calculations. Setting it to 0 prevents the storage of restart information. If running an MD calculation, see also section 2.3.6 regarding MDRestartFrequency. RandomSeed Sets the seed for the random number generator. The value 0 causes random initialisation. (This value can be used to reproduce earlier MD calculations by setting the initial seed to the same value.) MinimiseMemoryUsage Tries to minimise memory usage by storing various matrices on disc instead of keeping them in memory. Set it to Yes to reduce the memory requirement for calculations with many k-points or spin polarisation. Note: Currently this option has no effect and you will get a warning if setting it to be Yes. TimingVerbosity Level of information regarding CPU and wall clock timings of sections of the code, higher values becoming more verbose. Setting this parameter to 0 or below supresses any information being printed (default). Setting it to -1 includes all measured timings. ShowFoldedCoords Print coordinates folded back into the central cell, so if an atom moves outside the central cell it will reappear on the opposite side. The default behaviour is to use unfolded coordinates in the output. (Please note, that this option only influences how the coordinates are printed and written, it does not change the way, periodic systems are treated internally.) WriteHS Instructs the program to build the square Hamiltonian and overlap matrices and write them to files. The output files are hamsqrN.dat and oversqr.dat, where N enumerates the spin channels. For a detailed description of the file format see p. 65. Note: If either of the options WriteHS or WriteRealHS are set to Yes, the program only builds the matrices, writes them to disc and then stops immediately. No diagonalisation, no SCCcycles or geometry optimisation steps are carried out. You can use the ReadInitialCharges option to build the Hamiltonian with a previously converged charge distribution. WriteRealHS Instructs the program to build the real space (sparse) Hamiltonian and overlap matrices and write them to files. The output files are hamreal.dat and overreal.dat. For a detailed description of the file format see p. 65. Note: If either of the options WriteHS or WriteRealHS are set to Yes, the program only builds the matrices, writes them to disc and then stops immediately. No diagonalisation, no SCCcycles or geometry optimisation steps are carried out. You can use the ReadInitialCharges option to build the Hamiltonian with a previously converged charge distribution. ReadChargesAsText If No, the program expects the file charges.bin to contain starting charges stored in binary. If Yes, then charges.dat should contain a text file of this data. See section 3.7. WriteChargesAsText If No, the program stores charges in the binary file charges.bin, while if Yes then charges.dat contains text of this data. See section 3.7. 56 CHAPTER 2. INPUT FOR DFTB+ SkipChargeTest If Yes, testing of whether the charges read from file match the total charge (and magnetisation) specified in the DFTB+ input (if relevant) is performed. Skipping this test (setting to No) may be useful if restarting from a charges generated for a similar system with slightly different total charge or magnetisation. Similarly, in the event of serious instabilities in the SCC cycle, the generated charge restart file may fall outside of the check-sum tolerances, hence this option allows a re-start. Finally, in the case of user edited charges.dat file (see section 3.7), the check-sum this option removed the requirement that the checksum values in the file match the charges. 2.6 Analysis This block collects some options to analyse the results of the calculation and/or calculate properties. AtomResolvedEnergies MullikenAnalysis ProjectStates Localise WriteEigenvectors EigenvectorsAsTxt WriteBandOut CalculateForces ElectrostaticPotential l l m m l l l l m No Yes {} {} No No Yes No {} WriteEigenvectors = Yes SCC = Yes 58 AtomResolvedEnergies Specifies whether the contribution of the individual atoms to the total energies should be displayed or not. MullikenAnalysis If Yes, the results of a Mulliken analysis of the system is given. ProjectStates ProjectStates evaluates the Mulliken projection of electronic states onto specific regions of the system being modelled (partial density of states – PDOS). The format of the projected data files is similar to band.out, but the second column is the fraction of the state within that region, instead of its occupation number (for non-collinear and spin-orbit calculations, three additional columns for the magnetisation of the state are also given). Each region for projection is specified within a Region{} block, with the following options Atoms ShellResolved OrbitalResolved Label (i|s)+ l l s No No "regioni" ShellResolved Project onto separate atomic shells of the region. These are taken in order of increasing shell number of the atoms. ShellResolved = Yes is only allowed, if all the selected atoms are of the same type. OrbitalResolved Project onto separate atomic orbitals of the region. These are taken in order of increasing shell number of the atoms. As with ShellResolved, this only allowed, if all the selected atoms are of the same type. 2.6. ANALYSIS 57 Atoms Specification of the atoms over which to make the projection.. Atoms are specified in the same way as MovedAtoms in section 2.3.1.) Label Prefix of the label for the resulting file of data for this region. The default is “regioni.out” where i is the number of the region in the input. In the case that ShellResolved = Yes, the shell index is appended, so that files with names “Label.j.out” are written. For OrbitalResolved = Yes, the shell and then m-value is appended, so that files with names “Label.j.m.out” are written. Examples: ProjectStates = { Region = { # first region Atoms = 23:25 27 # atoms 23, 24, 25 and 27 } Region = { Atoms = N # All nitrogen atoms ShellResolved = Yes # s and p shells separated instead of atomic PDOS Label = "N" # files N.1.out and N.2.out for s and p states } } Localise Convert the single particle states of the calculation to localised orbitals via a unitary transformation. Localised orbitals span the same states as the occupied orbitals, so are equivalent to the usual valence band states, but are more localised in space. Currently only PipekMezey localisation is supported (but not for non-collinear or spin-orbit calculations). Pipek-Mezey [45] localisation transforms the occupied orbitals such that the square of the Mulliken charges for each orbital is maximised. The resulting localised states are output as localOrbs.out and localOrbs.bin following the format given in appendix 3.6 for eigenvec.out and eigenvec.bin. Tolerance MaxIterations r i 1E-4 100 Tolerance Cut off for rotations in the localisation process. MaxIterations Maximum number of total sweeps to perform. For systems with non-gamma-point k-points, no further options are available. Analysis = { Localise = { PipekMezey = { # These are the default options, which are also set if the bracket is left empty. Tolerance = 1.0E-4 MaxIterations = 100 } } } 58 CHAPTER 2. INPUT FOR DFTB+ For molecular and gamma point periodic calculations there are two implementations available, Dense = Yes will use the O(n4 ) scaling conventional algorithm, while Dense = No, uses the default sparse method which may have better scaling properties. Dense SparseTolerances l r+ Dense = No No 1E-1 1E-2 1E-6 1E-12 Dense Selects the conventional method (Yes) using Jacobi sweeps over all orbital pairs or (No) uses the default sparse method. SparseTolerances The sparse method introduces support regions during evaluation to increase performance, and these requires a set of tolerances to determine the regions to be used (these are listed in decreasing order, i.e., with tighter tolerances as the localisation proceeds). WriteEigenvectors Specifies, if eigenvectors should be printed in eigenvec.bin. For a description of the file format see p. 66. EigenvectorsAsTxt If eigenvectors are being written, specifies if a text version of the data should be printed in eigenvec.out. For a description of the file format see p. 66. WriteBandOut Controls the creation of the file band.out which contains the band structure in a more or less human friendly format. CalculateForces If Yes, forces are reported, even if not needed for the actual calculation (e.g. static geometry calculation). ElectrostaticPotential Evaluates the electrostatic potential at specified points in space for SCC calculations. This data is accumulated in a specified text file. OutputFile AppendFile Softening Points Grid s l r (3r)+ m MD or geometry optimisation Grid not set Points not set "ESP.dat" No 1E-6 {} {} OutputFile Text file to store the potential. If external electric fields are present, an additional column gives their values. See p. 67 for a description of the file. AppendFile If running calculations with multiple geometries, should the OutputFile be appended or only contain the last potential information? Softening [length] Modifies the plotted potential to remove the r = 0 divergence of 1/r, by setting ε √ and instead plotting 1/ r2 + ε 2 . Internal potential calculations are unaffected, only the exported data. Points [length] List of cartesian points at which to evaluate the electrostatic field. In the case Periodic = Yes, the modifier "F" may instead be used to specify the points as fractions of the lattice vectors. Grid [length] Specification of a regular 1, 2 or 3 dimensional grid of points. In the case Periodic = Yes, the modifier "F" may instead be used to specify the points as fractions of the lattice vectors. 2.7. EXCITEDSTATE GridPoints Origin Spacing Directions 3i 3r 3r 9r 59 Modifier not F 100 010 001 Spacing Separation between points in each direction. This inherits the modifier for Grid. Origin Location of first point in the grid. This inherits the modifier for Grid. GridPoints Number of points in each of the three direction of the grid (a value of 1 places all points at the Origin of that direction). Directions Set of 3 cartesian vectors along which the grid will become aligned. This can rotate, skew, etc. the grid. The vectors are internally normalised, but must be independent. 2.7 ExcitedState This block collects some options to calculate excited states. Casida 2.7.1 p SCC = Yes {} Casida This tag contains the specifications for a time-dependent DFTB calculation, based on linear response theory [40]. Note: the DFTB+ binary must be compiled with linear response calculations enabled to make use of these features (the ARPACK [32] library or ARPACK-ng [1] is required). The calculation of vertical excitation energies and the corresponding oscillator strengths as well as excited state geometry optimisation can be performed with these options, details of the resulting output files are given in appendix 3.10. Linear response theory is currently implemented only for the SCC-DFTB level of theory and molecular systems.5 Excitations can be calculated for fractional occupations and collinear spin-polarisation, but forces (and hence geometry optimisation or MD) are only available for spin-unpolarised systems with no fractional occupations. The specifications for this block have the following properties: 5 Excitation energies can also be calculated for gamma point periodic systems, but will be incorrect for delocalised excitations or for charge transfer-type excited states. 60 NrOfExcitations StateOfInterest Symmetry EnergyWindow OscillatorWindow WriteTransitions WriteSPTransitions WriteMulliken WriteCoefficients WriteEigenvectors TotalStateCoeffs WriteXplusY WriteTransitionDipole WriteStatusArnoldi TestArnoldi ExcitedStateForces CacheCharges CHAPTER 2. INPUT FOR DFTB+ i i s r r l l l l l l l l l l l l SpinPolarisation = {} WriteCoefficients = Yes CalculateForces = Yes 0 FORTRAN HUGE() -1 No No No No No No No No No No Yes Yes NrOfExcitations Specifies the number of vertical excitation energies to be computed for every symmetry (singlet or triplet). It is recommended that a value slightly greater than the actual number of the states of interest is specified (the eigenvalue solver may not converge to the right roots otherwise). StateOfInterest Specifies the target excited state or states that should be calculated. These are numbered from the first (lowest) excited state as 1, and so on. If the absorption spectrum at a given geometry is required (i.e., a single-point calculation), this parameter should be set to zero (default) and the Driver section (2.3) should be left empty (forces will not be available). A value less than 0 requests that the state with the largest dipole transition moment be found (again a single-point calculation). Symmetry Specifies the spin symmetry of the excited states being computed: “singlet”, “triplet” or “both”. This tag is only applicable for spin restricted calculation. For calculations in the “triplet” or “both” cases, SpinConstants must be supplied (see p. 34). EnergyWindow [energy] Energy range above the last transition at NrOfExcitations to be included in excited state spectrum calculation. OscillatorWindow [Dipole moment] Screening cut-off below which single particle transitions are neglected in excitation spectra calculations. This selects from states above the top of the EnergyWindow (if present). This keyword should not be used if calculating forces or other excited state properties. WriteTransitions If set to Yes, the file TRA.DAT is created. This file contains a description of each requested excited state in terms of its single-particle transitions. WriteSPTransitions If set to Yes, the file SPX.DAT is created, which contains the spectrum at the uncoupled DFTB level (i.e. the single-particle excitations). WriteMulliken If set to Yes, the files XCH.DAT and XREST.DAT are created. The former contains atom-resolved Mulliken (gross) charges for the excited state of interest, the latter the excitedstate dipole moment of the state. 2.8. PARSEROPTIONS 61 WriteCoefficients If set to Yes, the file COEF.DAT is created. This file contains the complex eigenvectors (molecular orbital coefficient) for the excited state of interest. They are derived from the relaxed excited state density matrix. WriteEigenvectors If set to Yes, the file excitedOrbs.bin is created. This file contains the natural orbitals for the specified excited state. TotalStateCoeffs Option to control data from WriteCoefficients or WriteEigenvectors. If set to No the total charge density of the output orbitals corresponds to the change in charge from the ground to excited state. If set to Yes instead it corresponds to the total charge density in the excited state. WriteXplusY If set to Yes, the file XplusY.DAT is created. This file contains the RPA vector (X +Y )IΣ ia for all excited states (c.f., Eqn. (18) in Ref. [23]). WriteTransitionDipole If set to Yes, the file TDP.DAT is created. This file contains the Mulliken transition dipole for each excited state. WriteStatusArnoldi If set to Yes, the file ARPACK.DAT is created, which allows the user to follow the progress of the Arnoldi diagonalisation. TestArnoldi If set to Yes, the file TEST_ARPACK.DAT is created, which gives data on the quality of the resulting eigenstates. ExcitedStateForces If set to Yes, evaluated forces include the contributions from an excited state of interest. By default, it is set to Yes if forces are being calculated (for example in geometry optimisation) and to No otherwise. By setting it explicitly to No, you can calculate the excitations during a molecular dynamics simulation that is being driven by the ground state forces only. CacheCharges If set to No, transition charges are calculated on the fly during the excited states calculation, instead of being cached. This makes the calculation considerably slower, but can help to decrease memory use substantially, if you are short on memory. 2.8 ParserOptions This block contains the options, which are effecting only the behaviour of the HSD/XML parser and are not passed to the main program. ParserVersion WriteHSDInput WriteXMLInput IgnoreUnprocessedNodes StopAfterParsing i l l l l current input version Yes No No No ParserVersion Version number of the input parser, which the input file was written for. If you are using an input file, which was created for an older version of DFTB+, you should set it to the parser version number of that code version. (The parser version number is printed at the beginning of the program run to the standard output.) DFTB+ internally converts the input to its current format. The processed input (written to dftb_pin.hsd) is always in the current format, and the ParserVersion property in it is always set to be the current parser version. 62 CHAPTER 2. INPUT FOR DFTB+ WriteHSDInput Specifies, if the processed input should be written out in HSD format. (You shouldn’t turn it off without really good reasons.) WriteXMLInput Specifies, if the processed input should be written out in XML format. IgnoreUnprocessedNodes By default the code stops if it detects unused or erroneous keywords in the input, which probably indicates error(s) in the input. This dangerous flag suspends these checks. Use only for debugging purposes. StopAfterParsing If set to Yes, the parser stops after processing the input and written out the processed input to the disc. It can be used to make sanity checks on the input without starting an actual calculation. 2.9 Parallel This block contains the options, which are effecting the parallel behaviour of the code. They only take effect, if the code was compiled with MPI-support. Groups UseOmpThreads Blacs i l p 1 .false. {} Groups Number of process groups. Specifying more than one process group enables parallelisation over k-points and spin, as processes in different process groups are working on different k-points and spins at the same time. The number of process groups must be a divisor of the total number of MPI-processes. Default: 1 (all processes work at the same k-point and spin at a given time). Note that transport calculations between contacts are currently incompatible with multiple process groups (see section 4). UseOmpThreads Enables the usage of OpenMP-threads (hybrid MPI/OpenMP-parallelisation). In order to prevent you from accidently running more processes and threads than appropriate for your hardware, this feature is turned off by default. Consequently in this case the MPIparallelised binary will stop if the maximal number of OpenMP-threads is greater than one when DFTB+ is started. (You can usually set the number of maximally allowed OpenMPthreads by setting the OMP_NUM_THREADS environment variable in your shell.) You can enable this option if you wish to run DFTB+ with hybrid parallelisation. You would then typically start fewer MPI-processes than physical cores on each node and also set the number of threads accordingly. This is currently an experimental feature in DFTB+ and is recommended for experienced users only. Blacs Contain BLACS specific settings. Currently only supports BlockSize, which specifies the row and column block size for the block-cyclic distributions (with default size of 32). Example: Parallel { Groups = 2 Blacs { BlockSize = 64 } } Chapter 3 Output of DFTB+ This chapter contains the description of some of the output files of DFTB+ where the output format is not self documenting. Unless indicated otherwise, numbers in the output files are given in atomic units (with Hartree as the energy unit). 3.1 band.out This contains the band energies and occupation of levels in electron volts and electron charge units as columns one and two. The file is printed if WriteBandOut = Yes (see section 2.6). Blocks of numerical results start with a line which labels the k-point and spin channel for the energies. See the DP _ TOOLS package for utilities for converting the data in this file into band-structures and density of states information suitable for plotting. 3.2 detailed.out This file contains details of the total energy and its components, as well as optional information on forces, atomic charges and other properties. It is intended for quick viewing, while values given to more significant figures are available in results.tag. Some of the information available in the file will also depend on the method being used in the calculation. For example, not all electronic solvers make the ground state electronic entropy available, hence only the internal energy would be quoted. Similarly, while the free energy of the system which when differentiated by atomic coordinates or boundary conditions gives the forces or stresses (printed as Force related energy) this is not currently available for some types of non-equilibrium transport calculations. Some of the common energy results printed in this file are: 63 64 CHAPTER 3. OUTPUT OF DFTB+ TS Total Electronic energy Repulsive energy Total energy Extrapolated to 0 Total Mermin free energy Force related energy Gibbs free energy MD Kinetic Energy Total MD Energy Product of the electron entropy and temperature The non-SCC energy plus other contributions to the electronic energy (SCC, spin, . . .) The pairwise contribution to the total energy Sum of electronic energy Estimated zero temperature energy if at finite temperatures U − T S, relevant free energy at finite temperatures Free energy relevant to forces in the system Energy corrected by −pV , i.e. the pressure and volume Kinetic energy of atoms in molecular dynamics Sum of finite temperature electronic, repulsive and atomic kinetic energies Where available the Fermi level µ (i.e. the chemical potential of the electrons in the system) is also printed. For systems with an externally fixed Fermi level (i.e. where the total charge can change), this contribution is included in the Force related energy: ∆E = +qtotal µ, but for calculations with fixed numbers of electrons it is not included in this energy. Note: The total energy reference may not match your required case in some situations, for example a shift with respect to the average electrostatic potential (in periodic cases) or whether the chemical potential should be with respect to the valence band maximum may be needed (see for example the discussion in Ref. [31]). 3.3 results.tag This contains machine readable results labeled with the type and size of the data blocks. The results are given in atomic units and are formatted as: label :type:shape: The variable type is real, complex, integer or logical. The shape information is :ndim: size1 ,size2 ,. . .,sizendim : where ndim is the number of dimensions, organised with the Fortran convention and of size size1 ×size2 ×size2 × . . .. In the special case of scalar variables the shape is :0:. A typical example of mixed scalar and both one and two dimensional results would be similar to: 3.4. HAMSQRN.DAT, OVERSQR.DAT 65 mermin_energy :real:0: -0.672967201447815E+000 total_energy :real:0: -0.672879398682698E+000 forces :real:2:3,3 -0.243590222274811E+000 -0.199780753617099E-001 -0.000000000000000E+000 0.465478448963764E+000 -0.228550455811745E+000 -0.000000000000000E+000 -0.221888226688953E+000 0.248528531173455E+000 -0.000000000000000E+000 gross_atomic_charges:real:1:3 0.171448741143825E+000 -0.254714832621691E+000 0.832660914778645E-001 3.4 hamsqrN.dat, oversqr.dat The files hamsqrN.dat and oversqr.dat contain the square (folded) Hamiltonian and overlap matrices. The number N in the filename hamrealN.dat indicates the spin channel. For spin unpolarised calculation it is 1, for spin polarised calculation it is 1 and 2 for spin-up and spin-down, respectively while for non-collinear spin it is charge, x, y and z for 1, 2, 3 and 4. Spin orbit is not currently supported for this option. Only non-comment lines (lines not starting with "#") are documented: • Flag for signalling if matrix is real (REAL), number of orbitals in the system (NALLORB), number of kpoints (NKPOINT). For non-periodic (cluster) calculations, the number of kpoints is set to 1. • For every k-point: – Number of the k-point. For molecular (non-periodic) calculations only 1 k-point is printed. – The folded matrix for the given k-point. It consists of NALLORB lines × NALLORB columns. If the matrix is not complex (REAL is F), every column contains two numbers (real and imaginary part). The files are produced if requested by WriteHS = Yes (see section 2.5). 3.5 hamrealN.dat, overreal.dat The files hamrealN.dat and overreal.dat contain the real space Hamiltonian and overlap matrices. The number N in the filename hamrealN.dat indicates the spin channel. For spin unpolarised calculation it is 1, for spin polarised calculation it is 1 and 2 for spin-up and spin-down, respectively, while for non-collinear spin it is charge, x, y and z for 1, 2, 3 and 4. Spin orbit is not currently supported for this option. Note: The sparse format contains only the "lower triangle" of the real space matrix. For more details about the format and how to obtain the upper triangle elements, see reference [5]. Also note, that for periodic systems the sparse format is based on the folded coordinates of the atoms, resulting in translation vectors (ICELL) which look surprising at first glance. Only non-comment lines (lines not starting with "#") are documented: 66 CHAPTER 3. OUTPUT OF DFTB+ • Number of atoms in the system (NATOM) • For every atom: – Atom number (IATOM), number of neighbours including the atom itself (NNEIGH), number of orbitals on the atom (NORB) • For every neighbour of every atom: – Atom number (IATOM1), neighbour number (INEIGH), corresponding image atom to the neighbour in the central cell (IATOM2F), coefficients of the translation vector between the neighbour and its corresponding image (ICELL(1), ICELL(2), ICELL(3)). Between the coordinates of the neighbour rINEIGH and the image atom rIATOM2F the relation 3 rINEIGH = rIATOM2F + ∑ ICELL(i) ai i=1 holds, where ai are the lattice vectors of the supercell. – The corresponding part of the sparse matrix. The data block consists of NORB(IAT1) lines and NORB(IAT2F) columns. The files are produced if requested by WriteRealHS = Yes (see section 2.5). 3.6 eigenvec.out, eigenvec.bin These files contain the eigenvectors from the Hamiltonian, stored either as plain text (eigenvec.out) or in the native binary format of your system (eigenvec.bin). The plain text format file eigenvec.out contains a list of the values of the components of each eigenvector for the basis functions of each atom. The atom number in the geometry, its chemical type and the particular basis function are listed, followed by the relevant value from the current eigenvector and then the Mulliken population for that basis function for that level. The particular eigenvector, k-point and spin channel are listed at the start of each set of eigenvector data. In the case of non-collinear spin, the format is generalised for spinor wavefunctions. Complex coefficients for both the up and down parts of the spinors are given (instead of single eigenvector coefficient) followed by four values – total charge, then (x, y, z) magnetisation. The binary format file eigenvec.bin contains the (unique) runId of the DFTB+ simulation which produced the output followed by the values of the eigenvectors. The eigenvector data is ordered so that the individual components of the current eigenvector are stored, with subsequent eigenvectors for that k-point following sequentially. All k-points for the current spin channel are printed in this order, followed by the data for a second channel if spin polarised. The files are produced if requested by setting WriteEigenvectors = Yes, with EigenvectorsAsTxt being also required to produce the plain text file (see section 2.6 for details). 3.7 charges.bin / charges.dat The file charges.bin contains the orbitally-resolved charges for each atom. In later versions of DFTB+ this format includes a check sum for the total charge and magnetisation. In the case 3.8. MD.OUT 67 of orbital potentials (p. 42) the file also contains extra population information for the occupation matrices. This file is produced as part of the mechanism to restart SCC calculations, see sections 2.5 and 2.3.6. Equivalent data can also be present in the file charges.dat, but stored as plain text. The options WriteChargesAsText and ReadChargesAsText control which cases are generated and read respectively. Appendix H contains details of the contents of the file. 3.8 md.out This file is only produced for VelocityVerlet{} calculations (See p. 16). It contains a log of information generated during MD calculations, and appended every MDRestartFrequency steps. In the case of small numbers of atoms and long MD simulations it may be useful to set WriteDetailedOut to No and examine the information stored in this file instead. 3.9 Electrostatic potential data The output from evaluating the electrostatic potential. The first line consists of a comment mark followed by a logical variable as to whether there is an external electric field (or not), followed by 3 values for any regular grid pattern present in the system and the total number of points. If the data is gridded, the next four lines contain the origin and grid separation vectors in Ångstroms. The next line is a comment, then the locations and the potential experience for a positive charge due to the internal field plus optionally the external field (from point charges or homogeneous electric fields). Values are given in Volts. In the case of gridded data, the location field is omitted. For an example with a regular grid # T 1 1 11 # 0.000000000000E+00 -0.200000000000E+01 -0.200000000000E+01 # 0.200000000000E+01 0.000000000000E+00 0.000000000000E+00 # 0.000000000000E+00 0.200000000000E+01 0.000000000000E+00 # 0.000000000000E+00 0.000000000000E+00 0.200000000000E+01 # Internal (V) External (V) 0.173386318927E-10 0.314737193575E+00 In the case where there is no regular grid: # T 0 0 01 # Location (AA) Internal (V) External (V) 0.0000E+00 -0.2000E+01 -0.2000E+01 0.173386318927E-10 0.314737193575E+00 In the case where data is generated for multiple geometry steps, this is also shown in the label: # F 1 5 5 25 # 0.000000000000E+00 -0.200000000000E+01 -0.200000000000E+01 # 0.100000000000E+01 0.000000000000E+00 0.000000000000E+00 68 CHAPTER 3. OUTPUT OF DFTB+ # 0.000000000000E+00 0.100000000000E+01 0.000000000000E+00 # 0.000000000000E+00 0.000000000000E+00 0.100000000000E+01 # Internal (V) Geo 0 0.215249473376E-01 . . # Internal (V) Geo 10 0.215815549672E-01 . . 3.10 Excited state results files Several files are produced during excited state calculations depending on the particular settings from section 2.7. Note: in the case of degeneracies, the oscillator strengths depend on arbitrary phase choices made by the ground state eigensolver. Only the sum over the degenerate contributions is well defined for most single particle transition properties, and label ordering of states may change if changing eigensolver or platform. For the excited state, properties like the intensities for individual excitations in degenerate manifolds again depend on phase choices made by both the ground and excited eigensolvers. 3.10.1 ARPACK.DAT Internal details of the ARPACK solution vectors, see the ARPACK documentation [32] for details. 3.10.2 COEF.DAT Data on the projection of this specific excited state onto the ground state orbitals. For the specific exited state, the (complex) decomposition of its single particle states onto the ground state single particle levels, together with its fractional contribution to the full excited state are given. General format: 3.10. EXCITED STATE RESULTS FILES 69 TF 1 1.9999926523 2.0000000000 -0.1944475716 0.0000000000 -0.1196876988 0.0000000000 .... -0.1196876988 0.0000000000 -0.1944475703 0.0000000000 .... . . . 2 1.9999866161 2.0000000000 -0.2400145188 0.0000000000 -0.1767827333 0.0000000000 .... Legacy flags level 1, fraction of total WF, 2.0 real then imaginary projection of level 1 onto ground state 1, then ground state 2, etc. level 2 real then imaginary projection of state 2 . . . 3.10.3 EXC.DAT Excitations data including the energies, oscillator strength, dominant Kohn-Sham transitions and the symmetry. Example first few transitions for C4 H4 : w [eV] Osc.Str. Transition Weight KS [eV] Sym. ========================================= 5.551 5.592 0.5143882 0.0000000 11 -> 10 -> 12 12 1.000 1.000 4.207 5.592 S S Two examples of singlet transitions with energies of 5.551 and 5.592 eV. The first is dipole allowed, the second not. In both cases they are transitions primarily (weight of 1.000) to single particle state 12, and are of singlet character (“S”). In the case of spin-polarised calculations, an additional column of values are given instead of the symmetry, showing the level of spin contamination in the state (labelled as D), with typically states where a magnitude of less than 0.5 is usually considered reliable [16]. 3.10.4 SPX.DAT Single particle excitations (SPX) for transitions between filled and empty single particle states of the ground state. These are given in increasing single particle energy and show the oscillator strength and index of the Kohn-Sham-like states that are involved. # w [eV] Osc.Str. Transition ============================ 70 CHAPTER 3. OUTPUT OF DFTB+ 1 2 3 4 5 6 5.403 5.403 5.403 5.403 6.531 6.531 0.2337689 0.2337689 0.2337689 0.2337689 0.0000000 0.0000000 15 14 15 14 13 12 -> -> -> -> -> -> 16 16 17 17 16 16 3.10.5 TDP.DAT Detail of the magnitude and direction of the transition dipole from the ground to excited states. 3.10.6 TRA.DAT Decomposition of the transition from the ground state to the excited states. The energy and spin symmetry are given together with the contributions from each of the single particle transitions. 3.10.7 TEST_ARPACK.DAT Tests on the quality of the eigenvalues and vectors returned by ARPACK. For the ith eigen-pair, the eigenvalue deviation corresponds to the deviation from (hxi |H|xi i − εi ), The eigen-vector deviation is a measure of rotation of the vector under the action of the matrix: |(H|xi i − εi |xi i)|2 , the normalisation deviation is hxi |xi i − 1 and finally largest failure in orthogonality to other eigenvectors is given. Example: State Ei deviation 1 -0.19428903E-15 2 0.27755576E-16 3 -0.12490009E-15 3.10.8 Evec deviation 0.80601119E-15 0.85748374E-15 0.88607302E-15 Norm deviation 0.19984014E-14 0.48849813E-14 0.88817842E-15 Max non-orthog 0.95562226E-15 0.36924443E-15 0.60384195E-15 XCH.DAT Charges on atoms in the specified excited state. The top line contains the symmetry (Singlet or Triplet) and the number of the excited state. The next line is the number of atoms in the structure followed by some header text. Then on subsequent lines the number of each atom in the structure and its charge are printed. 3.10.9 XplusY.DAT Expert file with the RPA (X +Y )IΣ ia data for all the calculated excited states. Line 1: number of single particle excitations and the number of calculated excited states Line 2: Level number 1, nature of the state (S, T, U or D) then excitation energy (in Hartree) Line 3: expansion in the KS single particle transitions . . . 3.10. EXCITED STATE RESULTS FILES Line 2: Level number 2, nature of the state (S, T, U or D) then excitation energy (in Hartree) 3.10.10 XREST.DAT Dipole moment of the specified excited state in units of Debye. 71 72 CHAPTER 3. OUTPUT OF DFTB+ Chapter 4 Transport calculations Non-equilibrium Green’s function (NEGF) calculations are now available with DFTB+. Within this formalism it is possible to treat quantum mechanical systems with open boundary conditions, i.e. systems connected to external reservoirs and therefore quantum transport. A new specific Transport{} block has been added to specify the geometry of such transport problems. Additional solvers have been added to the Solver section to either fully solve the open boundary problem (using the keyword GreensFunction{}) or the transmission through the system (TransportOnly). Finally a real-space Poisson solver is available for self consistent charge calculations and electrostatic gates (within a new section, Electrostatics, using the keyword Poisson{}). 4.1 Definition of the geometry The input geometry for transport calculations can be a little tricky. In comparison to cluster or supercell boundary conditions, the geometry for transport calculation must also contain information about the contacts (external reservoirs). The contacting leads (or surfaces) are actually semi-infinite structures, supporting travelling waves. Unlike finite structures, where mo stationary current is possible, travelling waves can only exist in such open systems. The simulation is therefore partitioned into a device region and two or more contact regions. Rules to build a valid input geometry: 1. All device atoms must come first in the structure. 2. Each contact must comprise of two subsequent unit cells, called principal layers (PLs). The two PLs together give all information about the contact structure and in the following are referred generally as a “contact”. 3. A PL is a unit cell of the contacting lead that has interactions only with its nearest neighbour PLs in tight-binding terms (i.e. the Slater-Koster interactions only extend into immediately neighbouring PLs). 4. The ordering of the atoms within the two PLs of a contact must be consistent, in the sense that the two PLs must be exact periodic replicas of each other: If each PL comprise N atoms, the ith atom in the first PL must have a corresponding identical atom, i + N, in the second PL which is related by translation to the position of atom i. 5. The first PL in a contact should be always the one which is closer to the device region. 73 74 CHAPTER 4. TRANSPORT CALCULATIONS 6. All blocks should be contiguous in the structure and each atom must belong to one and only one region. 7. The geometry can be defined as a cluster or a supercell. In the first case is it understood that the contacts are one-dimensional wire leads. 8. If a structure is defined as being a supercell, only the lattice vectors that are transverse to the transport direction are meaningful. The periodicity specified along the transport direction is treated as a dummy vector (but must be present). 9. For each contact the periodicity along the transport direction is actually deduced from the separation between the two PLs (using the coordinate difference r(i + N) − r(i)). We refer to this vector as aligned along the contact direction. 10. All lattice vectors (including the contact direction vector) must be aligned parallel to one of the Cartesian axes x, y or z. In practice only rectangular cells are allowed in transport calculations at present. An example of a non-periodic device with contacts attached is shown in Figure 4.1. Contact 2nd PL Contact vector v Contact 1st PL Contact 1st PL Contact 2nd PL Contact Contact 2nd PL 1st PL Device 1st PL Device 2nd PL Device 3rd PL Contact vector v Figure 4.1: Example of a valid 3 contact device with principal layers marked. Note: The code does not currently check: if the device regions are consistently defined (rules 1 and 6); if the PL defined are really PLs (rule 3); or if the first PL defined is really the one closest to the device (rule 5). The code does check rules 4, 8, 9 and 10. The check for rules 4 and 9 is performed on the atomic coordinates, such that R2i+N = R1i + v ∀i ∈ PL (4.1) where R2i are atomic coordinates of atoms in the second PL of the contact, R1i are atomic coordinates of atoms in the first PL and v is the contact lattice vector. The equality is verified within an accuracy that can be set by the user (see below for PLShiftTolerance). Please take care when building structures and to cross-check them. Also consider looking at the examples distributed with the code. The input structure is often the first suspect when there are problems in transport calculations. 4.2. TRANSPORT 4.2 75 Transport The Transport section collects together the information needed whenever open boundary conditions are used. It contains the description of the partitioning of the system into a device and the contact regions and additional information needed to calculate the required self-energies associated with the contacts. The transport block contains the following properties: Device Contact Task UploadContacts p p m 75 76 76 An example transport geometry specification looks like: Transport { Device { AtomRange = 1 8 } Contact { Id = "source" AtomRange = 9 24 } Contact { Id = "drain" AtomRange = 25 40 } } Where the associated atomic geometries follow the rules of Section 4.1. In this specific example, there is only one principal layer in the device, but each contact contains two principle layers (atoms 9–16 and 17–24 in the “source” contact, atoms 25–32 and 33–40 in the “drain” contact). 4.2.1 Device{} The Device blocks contains the following properties: AtomRange FirstLayerAtoms 2i i+ 1 75 AtomRange defines the first and last atom of the device region. FirstLayerAtoms defines the first atom of PLs in the device region. By default there is only one layer (the entire device region). Alternatively the user can manually reorder and group the atoms in the structure into distinct layers for more efficient Green’s function calculations. The device layers, unlike the contact PLs, do not need to represent unit cell repetitions. If the device geometry has specified principal layers, these must be ordered in such a way that all the atoms within each of the layer are contiguous in space and adjacent layers must be placed next to each other in the structure. This ensures that the constructed hamiltonian and overlap are block tri-diagonal. Refer to [43] for a description of the efficient iterative Green’s function algorithm that can then be applied. 76 4.2.2 CHAPTER 4. TRANSPORT CALCULATIONS Contact{} The contact block contains the following properties: Id AtomRange PLShiftTolerance Temperature FermiLevel Potential WideBand LevelSpacing s 2i r r r r l r 1E-5 0.0 WideBand = Yes 0.0 No 0.735 The sections Device and Contact are used to define the atomic range of each region. The user can also assign a label (Id) to each contact that can be used later for cross referencing. In the section Contact the user can add a keyword that specifies the accuracy for the internal check of the PLs (tolerance for rule 4 of structures, i.e. that accuracy to which (4.1) must be satisfied). Id Assign a text label to the contact (must be 50 or fewer characters). AtomRange Defines the first and last atom of the device region. Note the contacts should be defined such that the atoms included in the range are in continuous increasing order in the structure. PLShiftTolerance [length] Used to set the absolute accuracy used to check principal layer (PL) consistency (see above). The default is 10−5 atomic units. Please be aware that using a large values may hide errors due to an inconsistent definition of the contacts, therefore it should not be modified. Temperature [energy] Specifies the electronic temperature of the contact (see a more detailed discussion after the section Änalysis¨). FermiLevel [energy] Specifies the Fermi level of the contact. Potential [energy] Specifies any additional electrostatic potential applied to the contact. The natural units of this quantity are a potential energy. WideBand Use the wide band approximation for the contact. If set to Yes, the surface Green’s function of the contact is not explicitly calculated but is instead assumed to be local and constant according to a specified density of states. LevelSpacing [energy] Specifies the inverse of the density of states (DOS) per atom to be used for the Wide Band approximation. As an example, the DOS of gold at the Fermi level is 0.05 eV−1 atom−1 , which corresponds to an energy spacing of 20 eV ≈0.735 Hartree (the default value). 4.2.3 Task = ContactHamiltonian{} The Task option is used to define which type of calculation should be performed. Before performing a transport calculation it is necessary to compute some equilibrium properties of the contacts by running a periodic boundary condition DFTB calculation. This necessary step must be carried out separately for each contact and can be done by specifying a Task=ContactHamiltonian block, as in the following example to calculate the source case. 4.2. TRANSPORT 77 Task = ContactHamiltonian { ContactId = source ContactSeparation [Angstrom] = 50.0 } When Task=ContactHamiltonian the following options can be defined ContactId ContactSeparation s r 1e3 ContactId Id label of the contact to be calculated. ContactSeparation [length] Dummy separation in transverse direction (see the following explanation). The contact calculation computes the bulk Hamiltonian, self-consistent charges (if SCC) and Fermi level for each contact. This is a usual DFTB+calculation for which appropriate parameters must be included in the input file. For supercell structures the calculation of the contact is performed using corresponding supercells in which the transverse lattice vectors are those specified in the Geometry tag and the lattice vector along the contact direction is deduced from the PL separations (rule 9). If the structure is defined as a cluster, the contact calculation is performed for a supercell in which the contact is treated as one-dimensional periodic wire with a surrounding vacuum region. However, since DFTB+does not support pure one- and two-dimensional calculations, dummy lattice vectors are defined for the two remaining directions. The default value for these lattice vectors is 1000 a.u. (527 Å), which should guarantee sufficient wire to wire distances to avoid Coulomb interactions. The user can specify an alternative contact separation using the keyword ContactSeparation placed in the ContactHamiltonian block. Each contact computation produces one output file called shiftcont_ContactId.dat which storing energy shifts and Mulliken charges that must be present in the working folder in all subsequent transport calculations. Note that during the contact calculation you will need to perform a k-point integration in the contact direction (as the contacts are semi-infinite). Whenever the system is defined as a cluster, DFTB+ will automatically extract the periodicity vectors from the geometry such that the first reciprocal vector will correspond to the contact direction. Therefore you must specify a k-point sampling for the periodic calculation by sampling along the first reciprocal lattice vector. As an example, if the structure is defined as a cluster (i.e., 1-dimensional wire leads), the source contact calculation will have an input file similar to: ... Task = ContactHamiltonian { ContactId = source } ... Hamiltonian = DFTB { ... KpointsAndWeights = SupercellFolding { 8 0 0 # sampling points here regardless of the transport direction 0 1 0 0 0 1 0.5 0.0 0.0 78 CHAPTER 4. TRANSPORT CALCULATIONS } } On the other hand, if your structure is defined as a supercell (as an example, a molecule with bulk contacts) and the transport direction is along the y direction, your the source contact calculation will have an input file similar to: ... Task = ContactHamiltonian { ContactId = source } ... Hamiltonian = DFTB { ... KpointsAndWeights = SupercellFolding { 4 0 0 # points in periodic direction 0 8 0 # points in transport direction 0 0 4 # points in periodic direction 0.5 0.5 0.5 } } This could seem confusing, but the underlining reasons is that in the cluster calculation the reciprocal lattice is set up by the code itself, while in the periodic calculation is set up by the user who can chose any arbitrary direction. Refer to the transport cookbook and to the distributed examples for further clarification. 4.2.4 Task = UploadContacts{} After the contact calculations, it is possible to perform actual transport calculations. This is activated simply specifying Task = UploadContacts, without additional options (Note if no task is specified, DFTB+ assumes UploadContacts is the required task in the transport block). In order to set up a proper transport calculation the user should also define the contacts’ Fermi levels (as printed in the files previously produced by using ContactHamiltonian tasks) and any required potential shift for each contact. Transport { Device { AtomRange = 1 8 } Contact { Id = "source" AtomRange = 9 24 FermiLevel [eV] = -8.4123 Potential = 0.0 } Contact { Id = "drain" 4.3. GREENSFUNCTION 79 AtomRange = 25 40 FermiLevel [eV] = -8.4123 Potential = 1.0 } Task = UploadContacts } Note: During the transport calculation you will not need to set up the k-point integration when the structure is defined as a cluster, just as in a regular DFTB+ calculation. For supercell calculations, integration perpendicular to the transport direction will need to be accurate, but the sampling grid can in the transport direction itself can have only a single value. In the special case where your device is a supercell but also wire-like, with a vacuum region lateral to its transport direction, the Gamma-point can be chosen: KPointsAndWeights = { 0 0 0 1.0 } 4.3 GreensFunction For calculations in open systems, instead of calculating the eigenstates of the system, a Green’s function method is used to obtain the density matrix of the system. The Green’s function (GF) solver can also be used for conventional supercell/cluster boundary conditions if required. In order to activate Green’s function calculations the user must define the keyword Solver = GreensFunction in the Hamiltonian section. The GF solver, either under equilibrium (no bias applied) or under nonequilibrium conditions, builds the density-matrix of the device region such that it is consistent with any contacts that are present. Strictly speaking the GF does not solve for the eigenstates of the system, however it logically substitutes the traditional construction of the density matrix from the eigenstates of the system, as would be obtained after the diagonalisation step. The usual DFTB+ self-consistent calculations can be driven using the GF solver. The following table gives the important parameters of the solver: Name Delta ContourPoints LowestEnergy FermiCutoff EnclosedPoles RealAxisStep RealAxisPoints SaveSurfaceGFs ReadSurfaceGFs FirstLayerAtoms FermiLevel LocalCurrents Type r 2i r i i r r l l i+ r l Condition RealAxisPoints=undefined RealAxisStep=undefined Transport = undefined Transport = undefined Default 1E-5 20 20 -2.0 10 3 6.65E-4 Page Yes No 1 No Note: For efficient GF calculation the device region must be partitioned into layers whose fundamental property is to interact with nearest-neighbour layers only (see section 4.1). 80 CHAPTER 4. TRANSPORT CALCULATIONS Delta [energy] A small positive imaginary delta used in the GF definition and required for the x contour integration. ContourPoints The number of points along the complex contour integration of the GF along the segments C and L (see contour integration in section 4.5). LowestEnergy [energy] The initial energy from which the integration starts. It should be low enough to ensure that all the electronic states are correctly included in the integration. The default is -2.0 Hartree (see contour integration). FermiCutoff Integer number setting the Fermi distribution cutoff in units of kT . It is read only if the Fermi distribution temperature is greater than 0 (see contour integration). EnclosedPoles The number of Poles enclosed in the contour. It is meaningful only in finite temperature calculations (see contour integration). RealAxisStep [energy] The energy step along the real axis integration for non-equilibrium calculations. Note: RealAxisStep and RealAxisPoints cannot both be defined at the same time. RealAxisPoints The number of points along the real axis integration needed in non-equilibrium calculations. The default depends on the electronic temperature and bias. Note: RealAxisStep and RealAxisPoints cannot both be defined at the same time. SaveSurfaceGFs As the SCC cycle usually needs to repeat the calculation of the Green’s function at given energy points and as the surface Green functions do not change during the SCC cycle, this flag allows for saving the surface Green functions to disk and so save computational time on every SCC cycle after the first. ReadSurfaceGFs Loads the surface Green’s function from a file at the the first SCC cycle. Note that this operation only makes sense if the energy integration points are identical to the calculation used to generate the surface Green’s function files. The code does not verify whether this condition is fulfilled. In general there is no need to modify the defaults for ReadSurfaceGFs and SaveSurfaceGFs. FirstLayerAtoms As described in Device block. Can be specified only if no Transport block exists. Note: the Green solver can be used also to calculate the density matrix when there are no open boundary conditions, for example to take advantage of the iterative scheme in quasi1d systems. In this case, a Transport block is not defined and therefore FirstlayerAtoms{} should be given in the GreensFunction block. Also, the Fermi level of the system must be known and provided to fill up the electronic states. FermiLevel [energy] Required to set the Fermi level used by the Green’s solver to fill up the electronic states, unless already specified by there being contacts already present. LocalCurrents if set to Yes, local bond-currents are computed using the non-equilibrium density matrix. This task is currently limited to non-periodic systems. The output is placed in a file lcurrent_u.dat (or lcurrent_d.dat depending on spin). The files are arranged in a table in order of increasing neighbour distance, Atom(i) x y z nNeighbours j1 Ii, j1 j2 Ii, j2 j3 Ii, j3 ... This file can be processed using the small code flux provided in tools/transport that helps in building plots for jmol. GreensFunction section example: 4.4. SOLVER = TRANSPORTONLY 81 Solver = GreensFunction { FirstLayerAtoms = 1 61 92 145 Delta [eV] = 1E-4 ContourPoints = 20 20 RealAxisPoints = 55 LowestEnergy [eV] = -60.0 FermiCutoff = 10 EnclosedPoles = 3 } 4.4 Solver = TransportOnly The GreensFunction block is used to solve the full self-consistent NEGF transport problem. However, the block TunnelingAndDos{} within the Analysis block (see below), can be used to calculate the transmission function according to the Landauer formula, without solving for the full density matrix. This can be applied even for calculations where the density matrix and charge densities are not computed. Similarly, in these cases the electrostatics block should be omitted (i.e. the Electrostatics = Poisson). The keyword Solver = TransportOnly is used to jump straight to the post-SCC analysis. Note: This option cannot be used together with calculations which require forces, including geometry relaxations or md calculations. 4.5 Contour integration Much of the computational work for transport is in the integration of the energy resolved density matrix, as represented via the NEGF matrix. The integration is efficiently performed with a complex contour integration and a real axis integration, as shown in Figure 4.2 and discussed in references [2, 42, 43]. All integrations are performed with Gaussian quadratures and the number Figure 4.2: Contour integration in the complex plane for the Green’s functions. The crosses represent poles of either Gr or the Fermi function. of points must be specified manually. The complex contour integration is subdivided into two sections: the first section is the arc of a circle, C , that can be computed with a few integration points (default 20); the second section is a line that intersects the contour and runs parallel to the real axis at a distance that depends on the number of poles of the Fermi function that are enclosed within the contour. Usually a good choice for the number of poles is between 3 and 5 (the default is 3). The poles are placed at the complex points zm = EF + i(2m + 1)πkB T and therefore are separated from 82 CHAPTER 4. TRANSPORT CALCULATIONS each other by 2πkB T , where kB is Boltzmann’s constant. At T = 300 K this corresponds to a separation of 156 meV. It should be noted that, as the temperature decreases, the separation between poles reduces. This makes the contour integration harder as it needs to walk across two singularities. At very low temperatures, T = 10 K, the separation is 5.2 meV. Below this temperature, the contour integration is treated as T = 0 in order to avoid numerical inaccuracies. The integration along the segment L extends up to Re [z] = EF +nkB T , where n is an integer number specified by the keyword FermiCutoff and has a default value of 10. In the limit T = 0 K the poles collapse into a non-analytic branch cut and the contour needs to be changed such that the second section of the complex contour becomes the arc of circle closing on the real axis. Finally, the real axis integration extends between the lowest and highest chemical potentials. The number of quadrature points should depend on the bias itself and can be set using RealAxisPoints or RealAxisStep. The default value is 1 pt/0.018 eV (actually 1500 pt/1 Hartree). In finite temperature calculations the segment is extended to include the Fermi cutoff by nkB T on both sides (µ1 − nkB T, µ2 + nkB T ). In this case the number of quadrature points are increased by assuming the same point density defined in the range (µ1 , µ2 ). Example: for a bias of 0.2 V, the default number of points is 0.2 · 1500/27.21139 = 11. At T = 300 K the interval is increased by 0.26 eV on both sides, therefore 0.26 · 1500/27.21139 = 14.33 which is truncated to 14 points, leading to a total of 38 points along the real axis. The use of the keyword RealAxisStep is usually more convenient because it ensures a consistent real axis integrations during, for example, a bias sweep. Note: The GF solver can be used also for calculations other than the transport context. In cases where the position of the Fermi Energy is known with good accuracy, the density matrix solver based on the GF can be used to compute the electronic properties of clusters and supercells. The recursive algorithm may be an efficient solution to large problems having an elongated one dimensional shape. 4.6 Spin-polarised transport Spin-polarised transport is not available yet. Collinear spin transport will be available soon, as all the needed machinery has been implemented and is undergoing debugging and testing. 4.7 Poisson solver The Poisson solver is a fundamental part of charge self-consistent non-equilibrium transport calculations and must be declared whenever an SCC NEGF calculation is performed using Electrostatics = Poisson. Under non-equilibrium conditions the self-consistent potential of the KS equations cannot be solved using the efficient γ-functional, but instead requires the definition of appropriate boundary conditions for the potentials imposed by the contacts. However, since the γ-functional is functionally related to a pure Hartree potential, it can be obtained in real space by solving a Poisson solver. The Poisson equation is solved in a box with hexahedral prism shape. This restriction is imposed by the Poisson solver being employed. This restricts calculations of supercell structures to orthorhombic super-lattices. An additional restriction is that the box sides must be aligned with the Cartesian axes, x, y, z. 4.7. POISSON SOLVER Name Verbosity PoissonBox BoxExtension MinimalGrid PoissonAccuracy AtomDensityTolerance AtomDensityCutoff CutoffCheck NumericalNorm SavePotential PoissonAccuracy MaxPoissonIterations BuildBulkPotential ReadOldBulkPotential OverrideDefaultBC OverrideBulkBC BoundaryRegion Gate MaxParallelNodes RecomputeAfterDensity PoissonThickness 83 Type i 3r r 3r r r r l l l r i l l m m m m m l r Condition Default 51 0.0 0.5 0.5 0.5 1E-7 1E-6 14.0 Yes No No 1E-6 60 Yes Yes none{} none{} global{} none{} none{} No Page 84 84 84 84 84 87 87 contacts = 1 Verbosity This parameter controls the level and amount of output messages and takes values ranging from 1 to 100. PoissonBox [length] Dimension of the Poisson box along the directions x, y and z. BoxExtension [length] With this value it is possible to tune the position of the box interface between the device and contacts. By default (BoxExtension=0.0) the boundary is placed at the midpoint between the last device atom and the first contact atom. MinimalGrid [length] The minimal requested grid spacing along x, y and z. The actual grid spacing chosen by the multigrid solver will be lower than this. AtomDensityTolerance In order to calculate the potential, the Mulliken charges are projected on the real space grid. This parameter defines the cutoff after which the charge is considered to vanish (i.e., the spatial extension of the projected charge). The default is 1E-6 e. Note that the contacts must be at least twice the length over which a projected Mulliken charge extends. If this conditions is not fulfilled and CutoffCheck is set to Yes, the code will exit with an error message. Setting this parameter to a lower value could allow shorter contacts to be defined in some cases. However this could lead to error in the potential and hence to spurious reflections, therefore it should be left at its default value (or changed very carefully). AtomDensityCutoff [length] Defines the atomic radius cutoff. This is an alternative to AtomDensityTolerance and directly specifies the distance over which charge density associated with an atom is considered to vanish. CutoffCheck If set to No, consistency between contact length and charge extension is not verified (see AtomDensityTolerance and/or AtomDensityCutoff). The default is Yes. As with AtomDensityTolerance, this parameter should not be changed unless you know exactly what you’re doing. 84 CHAPTER 4. TRANSPORT CALCULATIONS SavePotential Save the electrostatic potential to the file potential.dat and the charge density to charge_density.dat. Additional files Xvector.dat, Yvector.dat, Zvector.dat and box.dat are also created. These files can be converted to a cube file that can be visualised in jmol. See section 4.13 about transport tools. PoissonAccuracy Defines the accuracy for the approximate solution of the Poisson equation (default value 10−6 ). MaxPoissonIterations Defines the maximum number of iterations allowed for the solver. RecomputeAfterDensity When set to Yes, Poisson’s equation is solved again after the density matrix is created in order to make the electrostatic energy consistent with the newly updated charges. In transport calculations it is set to No by default in order to avoid the extra time spent on the Poisson step. This does not affect the SCC loop or other calculations apart from the electronic energy and forces. PoissonThickness In the special case of a single contact (cases like the end of semi-infinite wires or surfaces of crystals), the thickness of the Poisson box normal to the surface of the contact can be set with this command. Note: The Poisson box can be specified using the keyword PoissonBox. In calculations where two contacts face each other along the same axis, setting the box-size along this axis will has no effect (the code adjusts to the correct size internally). The PoissonBox keyword is redundant (and should not be specified) when the system is periodic, since in this case the box geometry is taken from the supercell lattice vectors. Numerical error in the potential will results in spurious discontinuities at the contact-device interfaces. The default tolerances should be sufficient to avoid this in most cases. Below is a a typical example of the whole Poisson block specification. Some of the keywords are described in the next subsections. Electrostatics = Poisson { PoissonBox [Angstrom] = 20.0 20.0 20.0 MinimalGrid [Angstrom] = 0.3 0.3 0.3 SavePotential = No BuildBulkPotential = Yes ReadOldBulkPotential = No BoundaryRegion = Global {} PoissonAccuracy = 1E-7 Gate = Planar{ GateLength_l [Angstrom] = 10.0 GateLength_t [Angstrom] = 20.0 GateDistance [Angstrom] = 7.0 GatePotential [eV] = 1.0 } } 4.7.1 Boundary Conditions The Poisson equation is solved imposing boundary conditions (BC) on the potential at the six faces of the Poisson Box. In transport calculations for non-supercell geometries comprising two contacts placed along the same axis, the BCs are chosen as follows: 4.7. POISSON SOLVER 85 Dirichlet fixed potentials on the two contact faces with values defined by the applied potentials (see UploadContacts). Neumann zero normal field on the remaining 4 lateral box faces In periodic supercells the BCs are: Dirichlet (fixed potentials) on the two contact faces with values defined by the applied potentials (see UploadContacts) and Periodic on the remaining 4 lateral box faces. In some specific cases Neumann BCs can be set on one contact. In order to do so it is necessary to use OverrideDefaultBC (see below). The device and contact potentials should smoothly join at the interface. In order to achieve this the code computes the bulk potential of each contact and uses the result as a BC on the contact face of the Poisson box. This is useful when the contact potential is not uniform due to charge rearrangements. Any externally applied contact potential (Potential) is added to the bulk potential. The user can deactivate this calculation by setting the keyword BuildBulkPotential to No. Note: The bulk potential is computed on a special box that has “lateral” sizes copied from the device box, and has the size of one PL along the contact direction. The BCs are—so to speak—inherited from the device region. In particular: 1. Along the contact direction periodic BCs are imposed on both faces. 2. On the other four faces the BCs are copied from the device region (supercell or cluster). 3. The user can override this setting using OverrideBulkBC (see below). 4. When all four faces inherit Neumann BC (default for the device region), these are ALL internally changed to Dirichlet, because the solver cannot handle this situation as it gives rise to a singular matrix. BuildBulkPotential (default: Yes) is used to calculate the electrostatic potential of the contacts and the result is used as a Dirichlet boundary condition on the contact face (superimposed to the contact potential). ReadOldBulkPotential Read a previously computed bulk potential from hard-disk. BoundaryRegion Specifies how the Dirichlet boundary conditions are treated on each contact face of the Poisson box. It can be Global, Square or Circle. Global means that the BC is applied to the entire face of the box, whereas the other keywords imply that the Dirichlet BC are applied on a cross-section projected on the contact face. This is useful for instance when handling nanowire contacts, for which it is not really correct to impose a constant potential on the whole face of the Poisson box. BufferLength [length] can be used to set the size of the boundary region beyond the atomistic size which is determined as the minimal circle or square containing all atoms of the contact cross-section. Name BufferLength Example: Type r Condition Default 9.0 Page 86 CHAPTER 4. TRANSPORT CALCULATIONS BoundaryRegion = Circle { BufferLength [Angstrom] = 3.0 } In some special cases it might be necessary to override the default BCs applied by the code to the Poisson equation. Currently this can be done using the keywords: OverrideDefaultBC and OverrideBulkBC. In the special case of a single contact, the boundary condition on the other side of the box to that contact is automatically over-ridden to be of Neumann type (but can still then be over-ridden with OverrideDefaultBC). OverrideDefaultBC block is used to override the BCs described above. It can be used to force Dirichlet or Neumann BCs along some specified directions or on one of the four lateral faces of the Poisson box. Boundaries is used to specify on which face different BCs must be imposed. Assuming contacts are aligned along z, the keyword can be set to be any of xmin, xmax, x, ymin, ymax or y. OverrideDefaultBC = Dirichlet { Boundaries = xmin } For instance, setting a Dirichlet BC on Boundaries = xmin imposes φ (x, y, z) = 0 on the face placed at x = xmin , while boundaries = xmax would impose φ (x, y, z) = 0 on the face placed at x = xmax . When Dirichlet needs to be forced on both faces it is possible to use either boundaries = xmin,xmax or simply boundaries = x. The same syntax can be used to impose conditions on more faces, using boundaries = x,y or boundaries = x,ymin. A similar strategy can be used to impose different boundary conditions on the contacts. For instance, a Neumann BC can be set on one contact face by using OverrideDefaultBC = Neumann { Boundaries = zmin } Note that the user should know which face of the Poisson Box corresponds to the desired contact. Furthermore, if the user sets Neumann at all contacts the Poisson solver will not converge (singular matrix) unless the Dirichlet condition is imposed somewhere else (e.g., a gate potential is present). It is also possible to override the default BCs when computing the bulk potential. OverrideBulkBC block is used to override bulk BC usually copied from the device region. Boundaries has the same meaning and syntax as in OverrideDefaultBC. For example by choosing OverrideBulkBC = Neumann { Boundaries = x, y } 4.8. PARALLELISATIONS 4.7.2 87 Electrostatic Gates The option Gate can be used to specify an electrostatic gate. Currently the available gate types are Planar and Cylindrical. There are some restrictions as the planar gate must be placed with its face parallel to the x-z plane, i.e., the gate direction must be along y. At the same time the transport direction should be along the z-axis (i.e. perpendicular to the gate). The latter is not really a restriction but it gives meaning to “longitudinal” (l) and “transverse” (t) in the geometrical definitions of the gate lengths. Example: Gate = Planar { GateLength_l [Angstrom] = 20.0 GateLength_t [Angstrom] = 20.0 GateDistance [Angstrom] = 7.0 GatePotential [eV] = 1.0 } Gate = Cylindrical { GateLength [Angstrom] = 10.0 GateRadius [Angstrom] = 7.0 GatePotential [eV] = 1.0 } The various options for the gates have the following meanings: GateLength_l [length] Sets the gate length along the transport direction (always assumed to be z). The gate is centred in the middle of the device region. GateLength_t [length] Sets the gate extent transverse to the transport direction (assumed to be x). The gate is centred in the middle of the device region. GateDistance [length] Sets the distance of the gate from the centre axis of the device region. GatePotential [energy] Sets the potential applied to the gate. GateRadius [length] For a cylindrical gate, sets the distance of the gate from the centre axis or gate radius. Name GateLenth_l GateLenth_t GateDistance GatePotential GateRadius Type r r r r r Condition Default 0.0 0.0 0.0 0.0 0.0 Page Note that the gate option has not be tested thoroughly and may still contain bugs. Please report any problems you encounter to the developers. 4.8 Parallelisations The code has been parallelised in two main parts. The Non-equilibrium Green’s functions are computed by distributing the energy points along the contour and real axis calculations. Contour 88 CHAPTER 4. TRANSPORT CALCULATIONS and real axis integrations are independent and separately distributed. Load balancing has to be taken care of by the user. For instance if ContourPoints = {20 20} (i.e. 40 in total) and RealAxisPoints = 60, by setting 10 MPI nodes, each node will handle 4 points along the contour and 6 points along the real axis. Mixed OpenMP/MPI calculations are possible. When compiling DFTB+ the user should link against threaded BLAS/MKL, rather than sequential. Numerical experiments show that best performance on multicore CPUs is generally obtained by running independent MPI processes on physical sockets and exploiting OpenMP multithreading within each socket. For instance NEGF can exploit threaded matrix-matrix products. The user can experiment by setting the environment variable OMP_NUM_THREADS. The Poisson solver itself has not been parallelised yet. Currently the assembly of the charge density on the real-space grid and the projection of the potential onto the atoms has been parallelised. Since the gathering of the charge density on each node can easily hit communication bottlenecks, the user can use the parameter MaxParallelNodes to control distributions of these calculations. The default is MaxParallelNodes=1, this can be increased until speedups are observed. MaxParallelNodes 1 i Note: At the moment, calculations with more than one process group (see section 2.9) are not supported. 4.9 Analysis The Analysis block is used to specify post-scf calculations such as tunnelling or projected DOS. Analysis{ TunnelingAndDOS{ EnergyRange [eV] = {-5.0 -3.0} EnergyStep [eV] = 0.02 } } 4.10 TunnelingAndDos This method block can be specified in Analysis 56 and it is used to calculate the transmission by means of the Caroli formula, the current by means of the Landauer formula and the density of states from the spectral function. This block can only be specified if an open boundary conditions system has been defined in Transport (see p.75). EnergyRange EnergyStep TerminalCurrents ContactTemperature Region WriteTunn WriteLDOS 2r r p Nr p l l Telec 89 Yes Yes 4.11. SETTING ELECTRONIC TEMPERATURE 89 EnergyRange [energy] Contains the energy range over which the transmission function and local density of states are computed. EnergyStep [energy] Is the energy sampling step for evaluating properties. TerminalCurrents{} in multi-terminal configurations is used to define the terminal across which current must be computed. The terminal pairs are defined by using the keyword EmitterCollector, for example: TerminalCurrents{ EmitterCollector = {"source" "drain"} EmitterCollector = {"source" "gate"} } The block TerminalCurrents may be omitted since the code automatically sets all possible independent combinations for the terminal currents. For example in a 4-contact calculations the currents are 1–2, 1–3, 1–4, 2–3, 2–4 and 3–4. ContactTemperature [energy] Specifies the electronic temperature for the contacts used in the calculation of currents. It expects an array of real values, one per contact, which following the order the contacts are listed in Transport. Region{} This block defines atomic ranges or orbitals on to which the local density of states is projected. The definition in the block follow the same syntax as a DFTB+ calculation without transport (see section 2.6). WriteTunn The transmission coefficients are written also to a separate file for quick reference. If set to No, the transmission coefficient are only written to DFTB+ output files (detailed.out and detailed.xml, autotest.tag). WriteLDOS same as above, but for the density of states. 4.11 Setting electronic temperature In the current state of the code the electronic temperature of the system can be set in different places. One place is within the Hamiltonian block in the Filling section. The Temperature specified here is effective for the whole device and applies to all contact calculations as well. Note that during contact calculations the temperature is read from the Filling section and not from the contact section. During contact calculations only Fermi filling can be used. When a temperature is specified in the Contact section of the Transport block, it overrides the system temperature specified in filling. Note that the present electronic behaviour for transport is going to be changed soon. It is also possible to specify a (different) temperature in the section TunnelingAndDOS, within the block Analysis. The latter applies only in the calculation of currents (integration of the transmission function). Although slightly inconsistent, in some cases it is useful to be able to set a somewhat larger temperature in the calculation of the density than used for property calculations (e.g. currents), as this helps the convergence of the self-consistent loop. 90 4.12 CHAPTER 4. TRANSPORT CALCULATIONS Troubleshooting transport The DFTB+ transport machinery is designed to calculate transport in structures with a large number of atoms. To take full advantage of the iterative algorithm, be sure that the system is correctly partitioned into Principal Layers, as described in section 2.6. Be aware that an incorrect partitioning will lead to wrong results. If you are not completely confident, you can run a calculation on a test system with and without principal layer partitioning in the device region and the results should be the same. On some systems, a Segmentation Fault error could occur while running relatively large structures. This can happen because the stack memory limit on your system has been exceeded (the Intel compiler for example can show this behaviour). You can troubleshoot this by setting a higher limit for the stack memory. In bash you can remove the stack memory limitation with the command line ulimit -s unlimited. 4.13 Transport Tools Some tools useful for transport calculations can be found in tools/misc/transport. buildwire This tool can be used to create a one-dimensional nanowire, ready for transport calculations. A Principal Layer must be defined as a gen file, complete with supercell information. The code needs as input the number of PLs in the device region and the direction of the device. The resulting geometry will include 2 PLs for each contact and the specified number of PL repeats in the device region. flux This can be used to visualise the local bond currents in a junctions. The code reads the output files lcurrents.dat and writes out a script for jmol with arrows of different length/thickness for the currents. makecube This program can be used to convert the real-space potential.dat or charge_density.dat files computed on the Poisson box to a cube file that can be plotted using jmol. makecube potential.dat [-r refpot] [-b boxfile xfile yfile zfile] Options: -r refpot provides a reference potential that is subtracted from potential.dat For instance it is possible to subtract the equilibrium potential from the bias cases. -b The code reads by default the files box.dat, X,Y,Zvector.dat, but different filenames can be given with this flag Once the cube file has been created it can be read into jmol and visualised using the following script, script colormap128.jmol load "structure.xyz" isosurface pl1 fullplane plane {-1.2 0.8 0 0} color range all colorscheme ’user’ ’potential.cube’ 4.13. TRANSPORT TOOLS 91 Notice that a 128 palette colour map is provided in the tool folder. Also note that the structure should be converted to xyz to be read into jmol. 92 CHAPTER 4. TRANSPORT CALCULATIONS Chapter 5 MODES The MODES program calculates vibrational modes using data created by DFTB+. 5.1 Input for MODES The input file for MODES must be named modes_in.hsd and should be a Human-friendly Structured Data (HSD) formatted file (see Appendix A). The program can read the input in XML instead of HSD format if the input file is modes_in.xml. The input file must be present in the working directory. As with DFTB+ to prevent ambiguity, the parser refuses to read in any input if two types of input file are present. The table below contains the list of the properties, which must occur in the input file modes_in.hsd: Name Geometry Hessian SlaterKosterFiles Type p|m p p|m Condition Default {} - Page 10 94 Default 1:-1 No No No No Page 94 Additionally optional definitions may be present: Name DisplayModes Atoms WriteHSDInput WriteXMLInput RemoveTranslation RemoveRotation Type p i+|m l l l l Condition Geometry Specifies the geometry for the system to be calculated. See p. 10. Hessian Contains the second derivatives matrix of the system energy with respect to atomic positions. See p. 94. SlaterKosterFiles Name of the Slater-Koster files for every atom type pair combination. See p. 39. DisplayModes Optional settings to plot the eigenmodes of the vibrations. See p. 94. 93 94 CHAPTER 5. MODES Atoms Optional list of atoms, ranges of atoms and/or the species of atoms for which the Hessian has been supplied. This must be equivalent to the setting you used for MovedAtoms in your DFTB+ input when generating the Hessian. WriteHSDInput Specifies, if the processed input should be written out in HSD format. (You shouldn’t turn it off without good reason.) WriteXMLInput Specifies, if the processed input should be written out in XML format. RemoveTranslation Explicitly set the 3 translational modes of the system to be at 0 frequency. RemoveRotation Explicitly set the rotation modes of the system to be at 0 frequency. Note, for periodic systems, this is usually incorrect (if used for a molecule full inside the central cell, it may be harmless). 5.1.1 Hessian{} Contains the second derivatives of the energy supplied by DFTB+, see p. 16 for details of the options to generate this data. The derivatives matrix must be stored as the following order: For the i, j and k directions of atoms 1 . . . n as ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ∂ 2E ... ∂ xi1 ∂ xi1 ∂ x j1 ∂ xi1 ∂ xk1 ∂ xi1 ∂ xi2 ∂ xi1 ∂ x j2 ∂ xi1 ∂ xk2 ∂ xi1 ∂ xkn ∂ xkn Note: for supercell calculations, the modes are currently obtained at the q = 0 point, irrespective of the k-point sampling used. 5.1.2 DisplayModes{} Allows the eigenvectors of the system to be plotted out if present PlotModes Animate XMakeMol i+|m l l 1:-1 Yes Yes PlotModes Specifies list of which eigenmodes should be plotted as xyz files. Remember that there are 3N modes for the system (including translation and rotation). Animate Produce separate animation files for each mode or a single file multiple modes where the mode vectors are marked for each atom. XMakeMol Adapt xyz format output for XMakeMol dialect xyz files. Chapter 6 WAVEPLOT The WAVEPLOT program is a tool for the visualisation of molecular orbitals. Based on the files created by a calculation performed by DFTB+ it is capable of producing three dimensional information about the charge distribution. The information is stored as cube files, which can be visualised with many common graphical tools (e.g. VMD or J MOL). The user controls WAVEPLOT through an input file, choosing which orbitals and charge distributions should be plotted for which spatial region. Since the information about the shape of the basis functions is usually not contained in the Slater-Koster files, the coefficients and exponents for the Slater type orbitals must be entered by the user as part of the input file. The WAVEPLOT tool offers the following plotting capabilities: • Total charge density. • Total spin polarisation. • Difference between the total charge density and the density obtained by the superposition of the neutral atomic densities (visualisation of the charge shift). • Electron density for individual levels. • Real and imaginary part of the wavefunctions for individual levels. 6.1 Input for WAVEPLOT The input file for WAVEPLOT must be named waveplot_in.hsd and should be a Human-friendly Structured Data (HSD) formatted file (see Appendix A) The program can read the input in XML instead of HSD format if the input file is waveplot_in.xml. The input file must be present in the working directory. To prevent ambiguity, the parser refuses to read in any input if both files are present. The table below contains the list of the properties, which must occur in the input file waveplot_in.hsd: 95 96 CHAPTER 6. WAVEPLOT Name Options DetailedXML EigenvecBin GroundState Basis Type p s s s p Condition Default Yes - Page 96 100 Options Contains the options for WAVEPLOT. See p. 96. DetailedXML Specifies the name of the file, which contains the detailed XML output of the DFTB+ calculation (presumably detailed.xml). EigenvecBin Specifies the name of the file, which contains the eigenvectors in binary format (presumably eigenvec.bin). GroundState Read ground or excited state occupation data from the detailed XML output. Basis Contains the definition of the Slater-type orbitals which were used as basis in the DFTB+ calculation. At the moment, due to technical reasons this information has to be entered by the user per hand. In a later stage, it will be presumably read in by WAVEPLOT automatically. See p. 100. Additionally optional definitions may also be present: Name ParserOptions 6.1.1 Type p Condition Default {} Page 102 Options This property contains the options (as a list of properties), which the user can set, in order to influence the behaviour of WAVEPLOT. The following properties can be specified: PlottedRegion NrOfPoints PlottedKPoints PlottedLevels PlottedSpins TotalChargeDensity TotalSpinPolarisation TotalChargeDifference TotalAtomicDensity ChargeDensity RealComponent ImagComponent FoldAtomsToUnitCell FillBoxWithAtoms NrOfCachedGrids Verbose RepeatBox ShiftGrid p|m 3i i+|m i+|m i+|m l l l l l l l l l i l 3i l periodic system complex wavefunction periodic system No No No No No No No No No -1 No {1 1 1} Yes 98 6.1. INPUT FOR WAVEPLOT 97 PlottedRegion Regulates the region which should be plotted. See p. 98. NrOfPoints Specifies the resolution of the equidistant grid on which the various quantities should be calculated. The three integers represent the number of points along the three vectors of the parallelepiped specifying the plotted region. The number of all calculated grid points is the product of the three integers. Example: NrOfPoints = { 50 50 50 } # 125 000 grid points PlottedKPoints The list of integers specified here represent the k-points, in which the molecular orbitals should be plotted. The first k-point in the original DFTB+ calculation is represented by "1". The order of the specified k-points does not matter. You can also use the specification of the form from:to to specify ranges. (For more details on range specification, see the MovedAtoms keyword in the DFTB+ manual.) The actual list of molecular orbitals to plot is obtained by intersecting the specifications for PlottedKPoints, PlottedLevels and PlottedSpins. The option is ignored if the original calculation was not periodic. Example: PlottedKPoints = 1 3 5 # The 1st, 3rd and 5th k-point is plotted PlottedLevels The list of integers specified here represent the states, which should be plotted. The first (lowest lying) state in the original DFTB+ calculation is represented by "1". The order of the specified states does not matter. You can also use the specification of the form from:to to specify ranges. (For more details on range specification, see the MovedAtoms keyword in the DFTB+ manual.) The actual list of molecular orbitals to plot is obtained by intersecting the specifications for PlottedKPoints, PlottedLevels and PlottedSpins. Example: PlottedLevels = 1:-1 # All levels plotted PlottedSpins The list of integers specified here represent the spins, for which the molecular orbitals should be plotted. The first spin in the original DFTB+ calculation is represented by "1". The order of the specified spins does not matter. You can also use the specification of the form from:to to specify ranges. (For more details on range specification, see the MovedAtoms keyword in the DFTB+ manual.) The actual list of molecular orbitals to plot is obtained by intersecting the specifications for PlottedKPoints, PlottedLevels and PlottedSpins. Example: PlottedSpins = 1 2 # Both spin-up and spin-down plotted ChargeDensity If true, the absolute square of the wavefunction is plotted for the selected molecular orbitals. RealComponent If true, the real component of the wavefunction is plotted for the selected molecular orbitals. ImagComponent If true, the imaginary component of the wavefunction is plotted for the selected molecular orbitals. This option is only parsed, if the wavefunctions in the DFTB+ calculation were complex. 98 CHAPTER 6. WAVEPLOT TotalChargeDensity If true, the total charge density of the system is plotted. TotalSpinPolarisation If true, the total spin polarisation of the system (difference of the spin up and spin down densities) is plotted. This option is only interpreted if the processed DFTB+ calculation was spin polarised. TotalChargeDifference If true, the difference between the total charge density and the charge density obtained by superposing the neutral atomic densities is plotted. TotalAtomicDensity If true, the superposed neutral atomic densities are plotted. FoldAtomsToUnitCell If true, the atoms are folded into the parallelepiped unit cell of the crystal. FillBoxWithAtoms If true, the geometry is extended by those periodic images of the original atoms, which falls in the plotted region or on its borders. It sets FoldAtomsToUnitCell to Yes. NrOfCachedGrids Specifies how many grids should be cached at the same time. The value -1 stands for as many as necessary to be as fast as possible. Since the plotted grids could eventually become quite big, you should set it to some positive non-zero value if you experience memory problems. Example: NrOfCachedGrids = 5 # Maximal 5 cached grids RepeatBox The three integers specify how often the plotted region should be repeated in the generated cube files. Since repeating the grid is not connected with any extra calculations, this is a cheap way to visualise a big portion of a solid. You want probably set the FillBoxWithAtoms option to Yes to have the atoms also repeated (otherwise only the plotted function is repeated). In order to obtain the correct picture, you should set the plotted region to be an integer multiple of the unit cell of the crystal. Please note, that the phase of the wavefunctions in the repeated cells will be incorrect, except in the Γ-point. Example: RepeatBox = { 2 2 2 } # Visualising a 2x2x2 supercell ShiftGrid Whether the grid should be shifted, so that the specified origin lies in the middle of a cell and the grid fills out the specified plotted region symmetrically. The default is Yes. If set to No, the specified grid origin will be at the edge of a cell. Verbose If true, some extra messages are printed out during the calculation. PlottedRegion Specifies the region, which should be included in the plot. You can specify it explicitly (as property list), or let WAVEPLOT specify it automatically using either the UnitCell{} or the OptimalCuboid{} methods. 6.1. INPUT FOR WAVEPLOT 99 Explicit specification Specifies origin and box size explicitly. Origin Box - 3r 9r Origin [length] Specifies the xyz coordinates of the origin as three real values. Box [length] Specifies the three vectors which span the parallelepiped of the plotted region. The vectors are specified sequentially (a1x a1y a1z a2x a2y a2z a3x a3y a3z ). You are allowed to specify an arbitrary parallelepiped with nonzero volume here. Please note, however, that some visualisation tools only handles cube files with cuboid boxes correctly. Example: PlottedRegion = { Origin = { 0.0 0.0 0.0 } Box [Angstrom] = { 12.5 12.5 -12.5 12.5 -12.5 12.5 -12.5 12.5 12.5 } } UnitCell{} For the periodic geometries, this method specifies the plotted region to be spanned by the three lattice vectors of the crystal. The origin is set to (0 0 0). For cluster geometries, the smallest cuboid containing all atoms is constructed. For a cluster geometry the UnitCell{} object may have the following property: MinEdgeLength r 1.0 MinEdgeLength [length] Minimal side length of the cuboid, representing the plotted region. This helps to avoid cuboids with vanishing edge lengths (as it would be the case for a linear molecule). Example: PlottedRegion = UnitCell { MinEdgeLength [Bohr] = 2.0 } OptimalCuboid{} Specifies the plotted region as a cuboid, which contains all the atoms and enough space around them, that no wavefunctions are leaking out of the cuboid. This object does not have any parameters. Example: PlottedRegion = OptimalCuboid {} 100 6.1.2 CHAPTER 6. WAVEPLOT Basis The basis definition is done by specifying the following properties: Resolution ElementName1 ElementName2 .. . r p p - 101 101 Resolution Specifies the grid distance used for discretising the radial wavefunctions. Setting it too small, causes a long initialisation time for WAVEPLOT. Setting it too high causes a very coarse grid with bad mapping and inaccurate charges. Values around 0.01 seem to work fine. (Units must be in Bohr.) ElementName1 Basis for the first atom type. The name of this property is the name of that atom type. ElementName2 Basis for the second atom type. The name of this property is the name of that atom type. Before describing the properties (and their sub-properties) in detail, the full basis definition for carbon (sp) and hydrogen (s) should be presented as example: Basis = { Resolution = 0.01 C={ # Basis of the C atom AtomicNumber = 6 Orbital = { # 2s orbital AngularMomentum = 0 Occupation = 2 Cutoff = 4.9 Exponents = { 6.00000 3.00000 1.50000 } Coefficients = { 1.050334389886e+01 2.215522018905e+01 9.629635264776e+00 -4.827678012828e+01 -5.196013014531e+00 -2.748085126682e+01 3.072783267234e+01 -1.007000163584e+01 8.295975876552e-01 } } Orbital = { # 2p orbital AngularMomentum = 1 Occupation = 2 Cutoff = 5.0 Exponents = { 6.00000 3.00000 1.50000 } Coefficients = { -2.659093036065e+00 -6.650979229497e+00 -1.092292307510e+01 2.190230021657e+00 -9.376762008640e+00 -5.865448605778e-01 8.208019468802e+00 -2.735743196612e+00 2.279582669709e-01 } } } H={ # Basis for the H atom 6.1. INPUT FOR WAVEPLOT 101 AtomicNumber = 1 Orbital = { # 1s orbital AngularMomentum = 0 Occupation = 1 Cutoff = 4.2 Exponents = { 2.00000 1.00000 } Coefficients = { 1.374518455977e+01 1.151435212242e+01 2.072671588012e+00 -1.059020844305e+01 3.160957468828e+00 -2.382382105798e-01 } } } } Basis for an atom type The actual basis for every atom type is specified as a property with the name of that type: AtomicNumber Orbital .. . - i p 101 AtomicNumber The atomic number of the species. This is not needed in the actual calculations, but for creating proper cube-files. Orbital Contains the parameters of the orbitals. For every orbital a separate Orbital block must be created. See below. Orbital For every orbital there is an orbital block which specifies the radial wavefunction. Thereby the following properties must be used: AngularMomentum Occupation Cutoff Exponents Coefficients - i r r r+ r+ AngularMomentum Angular momentum of the current orbital. (s – 0, p – 1, d – 2, f – 3) Occupation Occupation of the orbital in the neutral ground state. (Needed to obtain the right superposed atomic densities.) Cutoff Cutoff for the wave function. You should choose a value, where the value of 4πr2 |R(r)|2 drops below 10−4 or 10−5 . R(r) is the radial part of the wave function. If you do not have the possibility to visualise the radial wave function, take the half of the longest distance, for which an overlap interaction exists in the appropriate homonuclear Slater-Koster file. (Value must be entered in Bohr.) Exponents The radial wave function with angular momentum l has the form: nexp npow Rl (r) = ∑ ∑ ci j rl+ j−1 e−αi r i=1 j=1 (6.1) 102 CHAPTER 6. WAVEPLOT This property defines the multiplication factors in the exponent (αi ). Coefficients This property contains the coefficients ci j as defined in equation (6.1). The sequence of the coefficients must be as follows: c11 c12 . . . c1npow c21 c22 . . . c2npow . . . 6.1.3 ParserOptions This block contains the options, which are effecting only the behaviour of the HSD/XML parser and are not passed to the main program. IgnoreUnprocessedNodes StopAfterParsing l l No No IgnoreUnprocessedNodes By default the code stops if it detects unused or erroneous keywords in the input. This dangerous flag suspends these checks. Use only for debugging purposes. StopAfterParsing If set to Yes, the parser stops after processing the input and written out the processed input to the disc. It can be used to make sanity checks on the input without starting an actual calculation. Chapter 7 SETUPGEOM The program utility SETUPGEOM can help in preparing the input geometry for transport calculations, following the rules specified in the Transport section. Starting from a geometry that can be the output of a previous relaxation step or any other building step, SETUPGEOM can be used to specify the system partitioning into contacts and device regions and reorder the atom numbers such that the device is placed before the contacts. Additionally the tool reorders the atoms of the two PLs of each contact or can create the second PL if only one is specified. Finally, the device region is reordered and partitioned into PLs for more efficient Green’s function calculations. A practical example is discussed in DFTB+ recipes. 7.1 Input for SETUPGEOM The input of the code must be named setup_in.hsd and should be a Human-friendly Structured Data (HSD) formatted file (see Appendix A). The file is similar to the DFTB+ input, where just 2 sections are needed. The table below contains the list of the properties, that must occur in the input file: Name Geometry Transport Type p|m p Condition Default {} Page 10 Geometry Specifies the geometry for the system to be calculated. See p. 10. 7.1.1 Transport{} The transport block must specify the atoms in each contact. An example of the Transport section is reported below: Transport { Contact { Id = "source" Atoms [zeroBased] = {9:24 56:78} ContactVector = 0.0 0.0 3.78 SpecifiedPLs = 2 103 104 CHAPTER 7. SETUPGEOM } Contact { Id = "drain" Atoms [zeroBased] = {81:100} ContactVector = 0.0 0.0 3.78 SpecifiedPLs = 2 } Task = SetupGeometry{ SlaterKosterFiles = type2names{ ... } TruncateSKRange = { SKMaxDistance [AA] = 5.0 HardCutOff = Yes } } } Id Atoms SpecifiedPLs ContactVector TruncateSKRange SlaterKosterFiles s li i 3r p p 1 or 2 2 - Id Specifies a unique contact name. Atoms Sets the list of atoms belonging to the named contact. This list can be easily obtained using some external atomic manipulation tool like for instance Jmol. NOTE the modifier [zeroBased] specifing that the defined atom numbers starts from 0 rather than 1. Use [oneBased] or no modifier for normal numbering starting from 1. In general follows the rules for MovedAtoms in section 2.3.1. SpecifiedPLs Specifies the number of PLs given for the named contact. If this value is 2 (default) the total number of atoms in the contact are divided by 2 and the 2nd PL is reordered according to the 1st with the help of ContactVector. If this value is 1, the correct numbers of PLs are created according to the interaction cutoff distance. ContactVector Sets the translation vector connecting the 2 PLs along the transport direction. Since contact must be aligned to a cartesian axis, so must be this vector. Different contact can be in different directions. Also notice that the vector must be specified along the positive axis direction. TruncateSKRange This section is the same as that described in section 2.4. SlaterKosterFiles This section is the same as that described in section 2.4. The SK files are used to compute the cutoff distance. 7.1. INPUT FOR SETUPGEOM 7.1.2 105 Code output The code writes two files, processed.gen and transport.hsd. The first file is the processed geometry, reordered according to the needs of transport calculations. NOTE that coordinates are folded to the unit cell such that all fractional coordinates are in the range 0 to 1. The structure is first translated such that all absolute coordinates have posive values. This step is important in order to take properly into account the periodic images. The file transport.hsd contains the details of the geometry partitioning for transport, as described in the Transport section and that can be included in the input file. For convenience this file also contains the block TruncateSKRange in order to make the Hamiltonian consistent with the MaxSKCutoff set in there. 106 CHAPTER 7. SETUPGEOM Chapter 8 DFTB+ API You can compile DFTB+ into a library and access some of its functionality via an API. Currently the API offers high level access only: you can set the current geometry and extract energy and forces for that geometry. 8.1 Building the library In order to compile the DFTB+ library, issue make api After the compilation, the library (libdftb+.a) can be found in the api/mm directory of your build tree. You will find also all mod-files here. If you issue make install_api the library and the modfiles will be installed in the appropriate sub-directories under the install target directory. You can test the library functionality by issuing make test_api 8.2 General guidelines Although the DFTB+ library contains nearly all internal routines of the DFTB+ code, you should access the code functionality only via the provided API and not by calling internal routines directly. We aim to keep the API stable over time, but the internal routines themselves can change without notice between releases. The API version can be found in the API_VERSION file in the api/mm folder. We use semantic versioning, a change in the major (first) version number indicates backwards incompatible changes, while changes in the minor (second) version number indicate backwards compatible extensions of the API. When using the API, we suggest that ParserVersion should be set in order to ensure that you can maintain backwards compatibility with later versions of DFTB+. 107 108 CHAPTER 8. DFTB+ API The Fortran interface is documented in the source code file api/mm/mmapi.F90, while api/mm/capi.F90 gives the bindings for calling from C. DFTB+ uses atomic units internally, hence exchanged values should be in this unit system (however HSD formatted data can carry unit modifiers, see examples of input parsing for details). Appendix A The HSD format The Human-friendly Structured Data (HSD) format is a structured input format, which can be bijectively mapped onto a subset of the XML-language. Its simplified structure and notation should make it a more convenient user interface than reading and writing XML tags. This section contains a brief overview of the most important aspects of this format. An input file in the HSD format consists basically of property assignments of the form Property = value where the value value was assigned to the property Property. The value must be one of the following types (detailed description of each follows later on): • Scalar, such as – integer – real – logical – string • list of scalars • method • list of further property assignments An unquoted hash mark (#) is interpreted as the start of a comment, everything after it, up to the end of the current line, is ignored by the parser (hash marks inside of quotes are taken as literals not comments): # Entire line with comment Prop1 = "hell#oo" # Note, that the first hashmark is quoted! The name of the properties, the methods and the logical values are case insensitive, so the assignments 109 110 APPENDIX A. THE HSD FORMAT Prop1 = 12 prOP1 = 12 Prop2 = Yes Prop2 = YES are pairwise identical. Quoted strings (specified either as a value for a property or as a file name), however, are case sensitive. Due to technical issues, the maximal line length is currently limited to 1024 characters. Lines longer than this are chopped without warning. If a property, which should only appear once, is defined more than once, the parser uses the first definition and ignores all the other occurrences. Thus specifying a property in the input a second time, does not override the first definition. (For advanced use the HSD syntax also offers the possibility of conditional overriding/extending of previous definitions. For more details see A.6.) A.1 Scalars and list of scalars The following examples demonstrate the assignments with scalar types: SomeInt = 1 SomeInt2 = -3 SomeRealFixedForm = 3.453 SomeRealExpForm = 2.12e-45 Logical1 = Yes Logical2 = no SomeString = "this is a string value" As showed above, real numbers can be entered in either fixed or exponential form. The value for logical properties can be either Yes or No (case insensitive). Strings should always be enclosed in quotation marks, to make sure that they are treated as one string and that they are not interpreted by the parser: String1 = "quoted string" String2 = this value is actually a list of 9 strings # list of strings! String3 = "Method { ;" # This is a string assignment String4 = Method { # This is syntactically incorrect, since # it tries to assign a method to String4 A list of scalars is created by sequentially writing the scalars separated by one or more spaces: PlottedLevels = 1 2 3 Origin = 0.0 0.0 0.0 ConfirmItTwice = Yes Yes SpeciesNames = "Ga" "As" The assignments statements are usually terminated by the end of the line. If the list of the assigned values goes over several lines, it must be entered between curly (brace) brackets. In that case, instead of the line end, the closing bracket will signal the end of the assignment. It is allowed to put a list of scalars in curly brackets, even if it is only one line long. A.2. METHODS AND PROPERTY LISTS 111 PlottedLevels = { 123 } Origin = { 0.0 0.0 0.0 } Short = { 1 2 3 } If you want to put more than one assignment in a line, you have to separate them with a semi-colon: Variable = 12; Variable2 = 3.0 If a property should be defined as empty, either the empty list must be assigned to it or it must be defined as an empty assignment terminated by a semi-colon: EmptyProperty = {} EmptyProperty2 = ; Please note, that explicitly specifying a property to be empty is not the same as not having specified it at all. In the latter case, the parser substitutes the default value for that property (if there is a default for it), while in the first case it interprets the property to be empty. If a property without default value is not specified, the parser stops with an appropriate error message. A.2 Methods and property lists Besides the scalar values and the list of scalars, the right hand side of an assignment may also contain a method, which itself may contain one or more scalar values or further property assignments as parameters: Diagonaliser = LapackDAC {} # Method without further params PlottedLevels = Range { 1 3 } # Range needs two scalar params PlottedRegion = UnitCell { # UnitCell needs a property list MinEdgeLength = 1.0 # as parameter SomeOtherProperty = Yes } The first assignment above is an example, where the method on the right hand side does not need any parameters specified. Please note, that even if no parameters are required, the opening and closing brackets after the method are mandatory. If the brackets are missing, the parser interprets the value as a string. In the second assignment, the method Range needs only two integers as parameters, while for the method UnitCell several properties must be specified. A method may contain either nothing or scalars or property assignments, but never scalars and property assignments together. So the following assignment would be invalid: InvalidSpecif = SomeMethod { 123 Property1 = 12 "Some strings here" } 112 APPENDIX A. THE HSD FORMAT Very often a value for the property is represented by a list of further property assignments (as above, but without naming an explicit method beforehand). In that case, the property assignments must be put between curly brackets (property list): Options = { SubOption1 = 12 Suboption2 = "string" } A.3 Modifiers Each property may carry a modifier, which changes the interpretation of the assigned value: LatticeConstant [Angstrom] = 12.23 Here, the property LatticeConstant possesses the Angstrom modifier, so the specified value will be interpreted to be in Ångström instead of the default length unit. Specifying a modifier for a property which is not allowed to carry one leads to parsing error. The syntax of the HSD format also allows methods (used as values on the right hand side of an assignment) to carry modifiers, but this is usually not used in the current input structures. Sometimes, the assigned value to a property contains several values with different units, so that more than one modifiers can be specified. In that case, the modifiers must be separated by a comma. VolumeAndChargePerElement [Angstrom^3,au] = { 1.2 0.3 # first element 4.2 0.1 # second element } You have to specify either no modifier or all modifiers. If you want specify the default units for some of the quantities, you can omit the name of the appropriate modifier, but you must include the separating comma: # Specifying the default unit for the charge. Note the separating comma! VolumeAndChargePerElement [Angstrom^3,] = { 1.2 0.3 # first element 4.2 0.1 # second element } Specifying not enough or too many modifiers leads to parser error. A.4 File inclusion It is possible to include files in an HSD-formatted input by using the <<< and <<+ operators. The former includes the specified file as raw text without parsing it, while latter parses the included text: A.5. PROCESSING 113 Geometry = GenFormat { <<< "geo_start.gen" } Basis = { <<+ "File_containing_the_property_definitons_for_the_basis" } The file included with the <<+ operator must be a valid HSD document in itself. A.5 Processing After having parsed and processed the input file, the parser writes out the processed input to a separate file in HSD format. This file contains the internal representation for all properties, which can be specified by the user. In particular, all default values are explicitly set and all automatic definitions (e.g. ranges) are converted to their internal representations. Assuming the following example as input # Lattice contant specified in Angstrom. # Internal representation uses Bohr, so it will be converted. LatticeConstant [Angstrom] = 12.0 # This property is not set, as its commented out, so the # default value will be set for this (let’s assume, it’s Yes) #DoAProperJob = No # Plotted levels specified as a range with parameter 1:3. # This will be replaced by an explicit listing of the levels PlottedLevels = { 1:3 } the parsed and processed input (written to a special file) should look something like LatticeConstant = 22.676713499923075 DoAProperJob = Yes PlottedLevels = { 123 } If you want to reproduce your calculation later, you should use this processed input. It should give you identical results, even if the default setting for some properties had been changed in the code. Since the HSD format is mapped by the parser internally to an XML tree, most codes using this format allow (or hopefully will allow) to dumping out of the processed input in the XML format as well, and to use that later as an input, instead of the HSD formatted input. This is practical for people preferring to work with XML or if the input should be automatically generated by a script. A.6 Extended format As stated earlier, if a property, which should be defined only once, occurs more than once in the input, the parser uses per default the first definition and ignores all the others. Sometimes this is not 114 APPENDIX A. THE HSD FORMAT the desired behaviour, therefore, the HSD format also offers the possibility to override properties that were set earlier. This feature can be very useful for scripts which are generate HSD input based on some user provided template. By just appending a few lines to the end of the user provided input the scripts can make sure that certain properties are set correctly. Thus, the script can modify the user input, without having to parse it at all. The parser builds internally an XML DOM-tree from the HSD input. For every property or method name an XML tag with the same name (but lowercased) is created, which will contain the value of the property or the method. If the value contains further properties or methods, new XML tags are created inside the original one. Shortly, the HSD input is mapped on a tree, whereas the assignment and the containment (equal sign and curly brace) are turned into a parent-child relationships.1 As an example an HSD input and the corresponding XML-representation is given below: Level0Elem1 = 1 Level0Elem2 = { 1 2 3 } Level0Elem3 = { Level1Elem1 = 12 Level1Elem2 = Level2Elem1 { Level3Elem1 = "abcd" Level3Elem2 = { Level4Elem1 = 12 } } } 11 2 312"abcd"12 By prefixing property and method names, the default behaviour of the parser can be overridden. Instead of creating a new tag (on the current encapsulation level) with the appropriate name, it will look for the first occurrence of the given tag and will process that one. Depending of the prefix character, the tag is processed in the following ways: +: If the tag exists already, it’s value is modified, otherwise the parser stops. ?: If the tag exists already, it’s value is modified, otherwise the parser ignores the prefixed HSD construct. *: If the tag exists already, it’s value is modified, otherwise it is created (and then it’s value is modified). /: If the tag does not exist yet, it is created and modified, otherwise the prefixed HSD construct is ignored. !: The tag is newly created and modified. If it exists already, the old occurrence is deleted first. The way the value of the tag is going to be modified, is ruled by the constructs inside the prefixed property or method name. If the parser finds non prefixed constructs here, the appropriate tags are just added, otherwise the behaviour is determined by the rules above, just acting one level deeper in the tree. The following examples should make this a little bit more clear. 1 In the internal tree representation of the HSD input there is no difference between properties and methods, both are just elements capable to contain some value or further elements. The differentiation in the HSD input is artificial and is only for human readability (equal sign after property names, curly brace after method names), A.6. EXTENDED FORMAT 115 • Changing the value of Level0Elem1 to 3. If the element does not exist, it should be created with the value 3. !Level0Elem1 = 3 • Changing the value of Level0Elem3/Level1Elem1 to 21 (the slash indicates the parent-child relationship). If the element does not exist, stop with an error message: # Make sure the containing element exists. If yes, go inside, otherwise die. +Level0Elem3 = { # Set the value of Level1Elem1 or die, if it does not exist. +Level1Elem1 = 21 } Please note, that each tag in the path must be prefixed. Using the following construct instead of the original one # Not prefixed, so it creates a new tag with empty value Level0Elem3 = { # The new tag doesn’t contain anything, so the parser stops here +Level1Elem1 = 21 } would end with an error message. Since Level0Elem1 is not prefixed here, a tag is created for it with an empty value (no children). It does not matter, whether the tag already existed before or not, a new tag is created and appended as the last element (last child) to the current block. Then the parser is processing its value. Due to the +Level1Elem1 directive it is looking for a child tag . Since the tag was newly created, it does not contain any children, so the parser stops with an error message. • Create a new tag Level1Elem3 inside Level0Elem3 with some special value. If the tag already exists, replace it. # Modifing the children of Level0Elem3 or dying if not present +Level0Elem3 = { # Replacing or if not existent creating Level1Elem3 !Level1Elem3 = NewBlock { NewValue1 = 12 } This example also shows, that the value for the new property can be any arbitrary complex HSD construct. • Provide a default value "string" for Level0Elem3/Level1Elem2/Level2Elem1/Level3Elem1. If the tag is already present do not change its value. # Modify Level0Elem3 or create it if non-existent *Level0Elem3 = { # Modify Level1Elem2 and Level2Elem1 or create them if non-existent *Level1Elem2 = *Level2Elem1 { # Create Level3Elem1 if non-existent with special value. /Level3Elem1 = "string" } } 116 APPENDIX A. THE HSD FORMAT • If Level0Elem3/Level1Elem2 has the value Level2Elem1, make sure that Level3Elem1 in it exists, and has "" as value. If Level1Elem2 has a different value, do not change anything. # If Level0Elem3 is present, process it, otherwise skip this block ?Level0Elem3 = { # The same for the next two containers ?Level1Elem2 = ?Level2Elem1 { # Create or replace Level3Elem1 !Level3Elem1 = "" } } Appendix B Unit modifiers The DFTB+ code uses internally atomic units (with Hartree as the energy unit). The value of every numerical property in the input is interpreted to be in atomic units (au), unless the property carries a modifier. The allowed modifiers and the corresponding conversion factors are given below.1 (The modifiers are case insensitive). Length: Angstrom, AA (for Ångström) Meter, m pm Bohr, au 0.188972598857892E+01 0.188972598857892E+11 0.188972598857892E-01 1.000000000000000E+00 Mass: amu au da dalton 0.182288848492937E+04 1.000000000000000E+00 0.182288848492937E+04 0.182288848492937E+04 Volume: Angstrom∧ 3, AA∧ 3 meter∧ 3, m∧ 3 pm∧ 3 bohr∧ 3, au 0.674833303710415E+01 0.674833303710415E+31 0.674833303710415E-05 1.000000000000000E+00 Energy: Rydberg, Ry Electronvolt, eV kcal/mol Kelvin, K cmˆ-1 Joule, J Hartree, Ha, au 0.500000000000000E+00 0.367493245336341E-01 0.159466838598749E-02 0.316681534524639E-05 0.455633507361033E-05 0.229371256497309E+18 1.000000000000000E+00 1 The conversion factors listed here were calculated with double precision on i686-linux architecture. Depending on your architecture, the values used there may deviate slightly. 117 118 APPENDIX B. UNIT MODIFIERS Force: eV/Angstrom, eV/AA Joule/meter, J/m Hartree/Bohr, Ha/Bohr, au 0.194469064593167E-01 0.121378050512919E+08 1.000000000000000E+00 Time: femtosecond, fs picosecond, ps second, s au 0.413413733365614E+02 0.413413733365614E+05 0.413413733365614E+17 1.000000000000000E+00 Charge: Coulomb, C au, e 0.624150947960772E+19 1.000000000000000E+00 Velocity: au m/s pm/fs AA/ps 1.000000000000000E+00 0.457102857516272E-06 0.457102857516272E-03 0.457102857516272E-04 Pressure: Pa au 0.339893208050290E-13 1.000000000000000E+00 Frequency: Hz THz cmˆ-1 au 0.241888432650500E-16 0.241888432650500E-04 0.725163330219952E-06 1.000000000000000E+00 Electric field strength: v/m au 0.194469063788953E-11 1.000000000000000E+00 Dipole moment: CoulombMeter,Cm Debye au 0.117947426715764E+30 0.393430238326893E+00 1.000000000000000E+00 Appendix C Description of the gen format The general (gen) format can be used to describe clusters and supercells. It is based on the xyz format introduced with xmol, and extended to periodic structures. Unlike some earlier implementations of gen, the format should not include any neighbour mapping information. The first line of the file contains the number of atoms, n, followed by the type of geometry. C for cluster (non-periodic), S for supercell in Cartesian coordinates or F for supercell in fractions of the lattice vectors. The supercells are periodic in 3 dimensions. The second line contains the chemical symbols of the elements present separated by one or more spaces. The following n lines contain a list of the atoms. The first number is the atom number in the structure (this is currently ignored by the program). The second number is the chemical type from the list of symbols on line 2. Then follow the coordinates. For S and C format, these are x, y, z in Å, but for F they are fractions of the three lattice vectors. If the structure is a supercell, the next line after the atomic coordinates contains the coordinate origin in Å (this is ignored by the parser). The last three lines are the supercell vectors in Å. These four lines are not present for clusters. Example: Geometry of GaAs with 2 atoms in the fractional supercell format 2 F # This is a comment Ga As 11 0.0 0.0 0.0 22 0.25 0.25 0.25 0.000000 0.000000 0.000000 2.713546 2.713546 0. 0. 2.713546 2.713546 2.713546 0. 2.713546 Note The DFTB+ input parser as well as the dptools utilities will ignore any lines starting with the # comment mark. 119 120 APPENDIX C. DESCRIPTION OF THE GEN FORMAT Appendix D Atomic spin constants These are suggested values for some atomic spin constants (W values) as given in reference [29], only the first two decimal places of the finite spin constants are numerically significant. These constants may eventually be included in the Slater-Koster files directly. Check the documentation of the Slater-Koster files required for a calculation to decide whether to use the LDA or PBE-GGA spin constants. W H C N O Si S Fe (3d 7 4s1 ) Ni s s p s p s p s p d s p d s p d s p d s -0.064 -0.028 -0.024 -0.030 -0.026 -0.032 -0.028 -0.018 -0.013 0.000 -0.019 -0.016 0.000 -0.013 -0.009 -0.003 -0.009 -0.009 -0.003 LDA p -0.024 -0.022 -0.026 -0.025 -0.028 -0.027 -0.013 -0.012 0.000 -0.016 -0.014 0.000 -0.009 -0.011 -0.001 -0.009 -0.010 -0.001 d 0.000 0.000 -0.019 0.000 0.000 -0.010 -0.003 -0.001 -0.015 -0.003 -0.001 -0.017 121 s -0.072 -0.031 -0.025 -0.033 -0.027 -0.035 -0.030 -0.020 -0.015 0.002 -0.021 -0.017 0.000 -0.016 -0.012 -0.003 -0.016 -0.012 -0.003 PBE p d -0.025 -0.023 -0.027 -0.026 -0.030 -0.028 -0.015 -0.014 0.002 -0.017 -0.016 0.000 -0.012 -0.029 -0.001 -0.012 -0.022 -0.001 0.000 0.000 -0.032 0.000 0.000 -0.080 -0.003 -0.001 -0.015 -0.003 -0.001 -0.018 122 APPENDIX D. ATOMIC SPIN CONSTANTS Appendix E Slater-Kirkwood dispersion constants The following table contains recommended dispersion constants for some elements with the SlaterKirkwood dispersion model (see p. 46). The values have been tested for biological systems, C, N, O and H predominantly for DNA [12]. If you would like to calculate different systems or you’re looking for other elements, check references [36] and [26]. The values of the atomic polarisabilities and cutoffs are given for zero, one, two, three, four and more than four neighbours. Element O N C H P S Polarisability [Å3 ] 0.560 0.560 0.000 0.000 0.000 0.000 1.030 1.030 1.090 1.090 1.090 1.090 1.382 1.382 1.382 1.064 1.064 1.064 0.386 0.386 0.000 0.000 0.000 0.000 1.600 1.600 1.600 1.600 1.600 1.600 3.000 3.000 3.000 3.000 3.000 3.000 123 Cutoff [Å] 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.5 3.5 3.5 3.5 3.5 3.5 4.7 4.7 4.7 4.7 4.7 4.7 4.7 4.7 4.7 4.7 4.7 4.7 Chrg 3.15 2.82 2.50 0.80 4.50 4.80 Note PO4 only S, not SO2 124 APPENDIX E. SLATER-KIRKWOOD DISPERSION CONSTANTS Appendix F DftD3 dispersion constants These are suggested dispersion values for some of the DFTB parameterizations for use with the DftD3 model (see p. 47). The table below gives the old defaults for DFTB3, those for which the 3OB parameters were fitted at the halogen correction stage, along with choices for the various OB2 parameterizations for range separated calculations. Becke-Johnson damping a1 a2 s6 s8 old default 0.5719 3.6017 1.0 0.5883 3OB 0.746 4.191 1.0 3.209 OB2 (base) 0.717 2.565 1.0 0.011 OB2 (shift) 0.816 2.057 1.0 0.010 OB2 (split) 0.497 3.622 1.0 0.010 Note: for the H5 corrected model, the DftD3 zero-damping parameters given on page 51 should be used. 125 126 APPENDIX F. DFTD3 DISPERSION CONSTANTS Appendix G Atomic onsite constants These are suggested values for some atomic on-site correction constants as given in reference [16], see p. 53. W H↑↑ H↑↓ C↑↑ C↑↓ N↑↑ N↑↓ O↑↑ O↑↓ S↑↑ S↑↓ Ti↑↑ Ti↑↓ Au↑↑ Au↑↓ s s s p s p s p s p s p s p s p d s p d s p d s p d s p d s p d s 0.000 0.000 0.00000 0.04973 0.00000 0.10512 0.00000 0.06816 0.00000 0.12770 0.00000 0.08672 0.00000 0.14969 0.00000 0.07501 0.00398 0.00000 0.11653 0.03915 0.00000 0.02659 -0.00587 0.00000 0.06881 0.01239 0.00000 0.03752 0.00073 0.00000 0.06928 0.01339 PBE p d 0.04973 -0.01203 0.10512 0.02643 0.06816 -0.00879 0.12770 0.03246 0.08672 -0.00523 0.14969 0.03834 0.07501 0.00310 0.01100 0.11653 0.03058 0.04979 0.02659 -0.01297 -0.00523 0.06881 0.01640 0.01144 0.03752 -0.00505 -0.00002 0.06928 0.01677 0.01228 0.00398 0.01100 -0.01792 0.03915 0.04979 0.01582 -0.00587 -0.00523 -0.00750 0.01239 0.01144 0.02604 0.00073 -0.00002 0.00531 0.01339 0.01228 0.02519 127 128 APPENDIX G. ATOMIC ONSITE CONSTANTS Appendix H Description of restart files H.0.1 charges.bin / charges.dat Initial charges and the current orbital charges are stored in these files. Both contain the same information, but charges.bin is stored as unformatted binary data, while charges.dat is a text file. The first line of the file is: version tBlockCharges tImaginaryBlock nSpin CheckSum Where version is currently 3, tBlockCharges and tImaginaryBlock are logical variables as to whether real and imaginary block charges are present. nSpin is the number of spin channels (1, 2 or 4 for spin free, collinear and non-collinear) and checksum is the totals for the charges in each spin channel. The subsequent nAtom × nSpin lines contain the individual orbital occupations for each atom in spin channel 1 (then 2 . . . 4, if present). If tBlockCharges is true, then the on-site block charges for each atom and spin channel are stored, followed by the imaginary part if tImaginaryBlock is true. Examples of the contents of charges.dat are given below for an H2 O molecule in the yz aligned with its dipole along y. Using the mio-1-1 Slater-Koster set, this file would contain: 3FF 1 6.5926151655316767 0.70369241723366482 0.70369241723466003 8.0000000000000018 0.0000000000000000 0.0000000000000000 0.0000000000000000 When OrbitalResolved = No. So, this is version 3 of the format, without block charges and it is spin free with 8 electrons. The electronic charges are grouped into the lowest atomic orbitals of each atom in this case. There is some small numerical noise in some of these these values (< 10−14 ). With OrbitalResolved = Yes, the oxygen has 1.7 2s electrons and 4.83 2p orbitals (electrons listed in the lowest labelled state in each shell). 3FF 1 1.7335403452609417 0.71592615825295036 0.71592615825154060 8.0000000000000018 4.8346073382345685 0.0000000000000000 0.0000000000000000 While for a pseudo-SIC calculation, where the net spin is 0: 3TF 2 1.7566193972978825 0.66337902366501833 0.66337902366479040 8.0000000000000018 1.7147230821039328 0.0000000000000000 1.2018994732683752 129 2.0000000000000013 130 APPENDIX H. DESCRIPTION OF RESTART FILES 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 1.7566193972978825 -0.28455491415632111 7.9592308124161928E-014 -7.7482799007142168E-027 -0.28455491415632111 1.7147230821039328 5.6922736516833719E-014 4.1776281206242259E-026 7.9592308124161928E-014 5.6922736516833719E-014 1.2018994732683752 -4.6749196043609904E-016 -7.7482799007142168E-027 4.1776281206242259E-026 -4.6749196043609904E-016 2.0000000000000013 0.66337902366501833 0.66337902366479040 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 Note the 0 blocks for the spin polarisation in channel 2, and also that the diagonal of the block charges matches the orbital charges for the atoms. Appendix I Publications to cite The following publications should be considered for citation, if you are publishing any results calculated by using DFTB+. DFTB+ code non-SCC DFTB SCC DFTB Collinear spin polarisation Non-collinear spin polarisation Spin orbit coupling QM/MM coupling (external charges) Van der Waals interaction (dispersion) DFTB+U 3rd order corrections linear-response TD-DFTB [5] [46], [50] [13] [28] [27] [27] [9], [21] [12] [24] [53] [40] 131 132 APPENDIX I. PUBLICATIONS TO CITE Bibliography [1] https://github.com/opencollab/arpack-ng. 59 [2] Introducing Molecular Electronics, chapter Tight-Binding DFT for Molecular Electronics (gDFTB), pages 153–184. Lecture Notes in Physics. Springer, 2005. 81 [3] H. C. Andersen. Molecular dynamics at constant pressure and/or temperature. J. Chem. Phys., 72:2384, 1980. 17 [4] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. 35 [5] B. Aradi, B. Hourahine, and Th. Frauenheim. DFTB+, a sparse matrix-based implementation of the DFTB method. J. Phys. Chem. A, 111(26):5678, 2007. 65, 131 [6] B. Aradi, A. M. N. Niklasson, and T. Frauenheim. Extended lagrangian density functional tight-binding molecular dynamics for molecules and solids. J. of Chem. Theory Comput., 11:3357–3363, 2015. 20, 54 [7] H. J. C. Berendsen, J. P. M. Postma, W. F. Van Gunsteren, A. Dinola, and J. R. Haak. Molecular-dynamics with coupling to an external bath. J. Chem. Phys., 81(8):3684–3690, 1984. 18, 19, 20 [8] M. Ceriotti, J. More, and D. E. Manolopoulos. i-pi: A python interface for ab initio path integral molecular dynamics simulations. Computer Phys. Comm., 185:1019–1026, 2014. 23 [9] Q. Cui, M. Elstner, T. Frauenheim, E. Kaxiras, and M. Karplus. Combined self-consistent charge density functional tight-binding (SCC-DFTB) and CHARMM. J. Phys. Chem. B, 105:569, 2001. 131 [10] A. Domínguez, T. A. Niehaus, and T. Frauenheim. Accurate hydrogen bond energies within the density functional tight binding method. The Journal of Physical Chemistry A, 119(14):3535–3544, 2015. 28, 53 [11] Jack Dongarra, Mark Gates, Azzam Haidar, Jakub Kurzak, Piotr Luszczek, Stanimire Tomov, and Ichitaro Yamazaki. Accelerating numerical dense linear algebra calculations with gpus. Numerical Computations with GPUs, pages 1–26, 2014. 36 [12] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, and E. Kaxiras. Hydrogen bonding and stacking interactions of nucleic acid base pairs: a density-functional-theory based treatment. J. Chem. Phys., 114:5149, 2001. 45, 46, 47, 123, 131 133 134 BIBLIOGRAPHY [13] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B, 58:7260, 1998. 54, 131 [14] V. Eyert. A comparative study on methods for convergence acceleration of iterative vector sequences. J. Comp. Phys., 124:271, 1996. 30 [15] T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Kohler, M. Amkreutz, M. Sternberg, Z. Hajnal, A. Di Carlo, and S. Suhai. Atomistic simulations of complex materials: groundstate and excited-state properties. J. Phys. Cond. Matter, 14(11):3015–3047, Mar 2002. 7 [16] Adriel Dominguez Garcia. Density functional approaches for the interaction of metal oxides with small molecules. PhD thesis, Universität Bremen, 2014. http://elib.suub.unibremen.de/edocs/00103868-1.pdf. 53, 69, 127 [17] M. Gaus, Q. Cui, and M. Elstner. DFTB3: Extension of the Self-Consistent-Charge DensityFunctional Tight-Binding Method (SCC-DFTB). J. of Chem. Theory Comput., 7:931–948, 2011. 49, 50 [18] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements HPu. J. Chem. Phys., 132:154104, 2010. 45, 47 [19] S. Grimme, S. Ehrlich, and L. Goerigk. Effect of the damping function in dispersion corrected density functional theory. J. Chem. Phys., 32:1456–1465, 2011. 45, 47 [20] M. J. Han, T. Ozaki, and J. Yu. O(N) LDA+U electronic structure calculation method based on the nonorthogonal pseudoatomic orbital basis. Phys. Rev. B, 73:045110, 2006. 35 [21] W. Han, M. Elstner, K. J. Jalkanen, T. Frauenheim, and S. Suhai. Hybrid SCC-DFTB/molecular mechanical studies of H-bonded systems and of N-acetyl-(L-Ala)n N’-Methylamide helices in water solution. Int. J. Quant. Chem., 78:459, 2000. 131 [22] S. C. Harvey, R. K. Z. Tan, and T. E. Cheatham. The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition. J. Comp. Chem., 19(7):726– 740, 1998. 18 [23] D. Heringer, T. A. Niehaus, M. Wanko, and T. Frauenheim. Analytical excited state forces for the time-dependent density-functional tight-binding method. J. Comp. Chem., 28(16):2589, 2007. 61 [24] B. Hourahine, S. Sanna, B. Aradi, C. Köhler, T. Niehaus, and Th. Frauenheim. Self-interaction and strong correlation in DFTB. J. Phys. Chem. A, 111(26):5671, 2007. 42, 131 [25] D. D. Johnson. Modified Broyden’s method for accelerating convergence in self consistent calculations. Phys. Rev. B, 38:12807, 2003. 29 [26] Y. K. Kang and M. S. Jhon. Additivity of atomic static polarizabilities and dispersion coefficients. Theoretica Chimica Acta, 61:41, 1982. 123 [27] C. Köhler, T. Frauenheim, B. Hourahine, G. Seifert, and M. Sternberg. Treatment of collinear and noncollinear electron spin within an approximate density functional based method. J. Phys. Chem. A, 111(26):5622, 2007. 131 [28] C. Köhler, G. Seifert, and T. Frauenheim. Density-functional based calculations for Fe(n), (n≤32). Chem. Phys., 309:23, 2005. 131 BIBLIOGRAPHY 135 [29] Christof Köhler. Berücksichtigung von Spinpolarisationseffekten in einem dichtefunktionalbasierten Ansatz. PhD thesis, Department Physik der Fakultät fur Naturwissenschaften an der Universität Paderborn, 2004. http://ubdata.uni-paderborn.de/ediss/06/2004/koehler/. 121 [30] A. Kovalenko, S. Ten-no, and F. Hirata. Solution of three-dimensional reference interaction site model and hypernetted chain equations for simple point charge water by modified method of direct inversion in iterative subspace. J. Comp. Chem., 20:928–936, 1999. 15 [31] Stephan Lany and Alex Zunger. Accurate prediction of defect properties in density functional supercell calculations. Modelling and Simulation in Materials Science and Engineering, 17(8):084002, Nov 2009. 64 [32] R. B. Lehoucq, D. C. Sorensen, and C. Yang. Arpack users guide: Solution of large scale eigenvalue problems by implicitly restarted arnoldi methods., 1997. 59, 68 [33] V. Lutsker, B. Aradi, and T. A. Niehaus. Implementation and benchmark of a long-range corrected functional in the densityfunctional based tight-binding method. J. Chem. Phys., 143:184107, 2015. 52 [34] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein. Explicit reversible integrators for extended systems dynamics. Molecular Phys., 87:1117–1157, 1996. 18, 19 [35] M. Methfessel and A. T. Paxton. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B, 40:3616, 1989. 39 [36] K. J. Miller. Additivity methods in molecular polarizability. J. Am. Chem. Soc., 112:8533, 1990. 123 [37] H. J. Monkhorst and J. D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B, 13:5188, 1976. 41 [38] H. J. Monkhorst and J. D. Pack. "special points for Brillouin-zone integrations"–a reply. Phys. Rev. B, 16:1748, 1977. 41 [39] T. A. Niehaus and F. Della Sala. Range separated functionalsin the density functional based tight-binding method: Formalism. Phys. Status Solidi B, 249(2):237–244, 2012. 52 [40] T. A. Niehaus, S. Suhai, F. Della Sala, P Lugli, M. Elstner, G. Seifert, and Th. Frauenheim. Tight-binding approach to time-dependent density-functional response theory. Phys. Rev. B, 63:085108, 2001. 59, 131 [41] Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, New York, NY, USA, second edition, 2006. 15 [42] A. Pecchia and A. Di Carlo. Atomistic theory of transport in organic and inorganic nanostructures. Rep. Prog. Phys., 67:1497, 2004. 81 [43] A. Pecchia, G. G Penazzi, L. Salvucci, and A. Di Carlo. Non-equilibrium Green’s functions in density functional tight binding: method and applications. New J. of Physics, 10(6):065022, 2008. 75, 81 [44] A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein. Correlated metals and the LDA+U method. Phys. Rev. B, 67:153106–4, 2003. 42 136 BIBLIOGRAPHY [45] J. Pipek and P. G. Mezey. A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys., 90:4916, 1989. 57 [46] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner. Construction of tightbinding-like potentials on the basis of density-functional theory: Application to carbon. Phys. Rev. B, 51:12947, 1995. 131 [47] A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard III, and W. M. Skiff. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc., 114:10024–10035, 1992. 44, 45 [48] Jan Řezáč. Empirical Self-Consistent Correction for the Description of Hydrogen Bonds in DFTB3. J. of Chem. Theory Comput., 13(10):4804–4817, 2017. 49, 50, 51 [49] Jan Řezáč and Pavel Hobza. Advanced corrections of hydrogen bonding and dispersion for semiempirical quantum mechanical methods. J. of Chem. Theory Comput., 8:141–151, 2012. 49, 50, 51 [50] G. Seifert, D. Porezag, and T. Frauenheim. Calculations of molecules, clusters, and solids with a simplified LCAO-DFT-LDA scheme. Int. J. Quant. Chem., 58:185, 1996. 131 [51] Stanimire Tomov, Jack Dongarra, and Marc Baboulin. Towards dense linear algebra for hybrid GPU accelerated manycore systems. Parallel Computing, 36(5-6):232–240, June 2010. 36 [52] Stanimire Tomov, Rajib Nath, Hatem Ltaief, and Jack Dongarra. Dense linear algebra solvers for multicore with GPU accelerators. In Proc. of the IEEE IPDPS’10, pages 1–8, Atlanta, GA, April 19-23 2010. IEEE Computer Society. DOI: 10.1109/IPDPSW.2010.5470941. 36 [53] Y. Yang, H. Yu, D. York, Q. Cui, and M. Elstner. Extension of the self-consistent-charge density-functional tight-binding method: Third-order expansion of the density functional theory total energy and introduction of a modified effective coulomb interaction. J. Phys. Chem. A, 111:10861, 2007. 49, 50, 131 [54] V. W.-z. Yu, F. Corsetti, A. García, W. P. Huhn, M. Jacquelin, W. Jia, B. Lange, L. Lin, J. Lu, W. Mi, A. Seifitokaldani, Á Vázquez-Mayagoitia, C. Yang, H. Yang, and V. Blum. ELSI: A unified software interface for Kohn-Sham electronic structure solvers. Computer Phys. Comm., 222:267–285, 2018. 36 [55] L. Zhechkov, Th. Heine, S. Patchkovskii, G. Seifert, and H. A. Duarte. An efficient a posteriori treatment for dispersion interaction in density-functional-based tight binding. J. of Chem. Theory Comput., 1:841–847, 2005. 45 Index AdaptFillingTemp, 18 AllAtomCharges, 27 AllAtomSpins, 32, 33 Alpha, 15 Analysis, 10, 88 AngularMomentum, 101 Animate, 94 AppendFile, 58 AppendGeometries, 12, 14, 15 ARPACK.DAT, 68 Atom list, 12 AtomCharge, 27 AtomDensityCutoff, 83 AtomDensityTolerance, 83 AtomicNumber, 101 AtomRange, 75, 76 AtomResolvedEnergies, 56 Atoms, 16, 22, 32, 33, 56, 93, 104 AtomSpin, 32, 33 band structure calculation, 42 band.out, 63 Barostat, 16 Basis, 96 Blacs, 62 Boundary Conditions, 84 BoundaryRegion, 83 Box, 99 BoxExtension, 83 Brillouin-zone sampling, 40 BufferLength, 85 BuildBulkPotential, 83, 85 CacheCharges, 60 CalculateForces, 56 Casida, 59 ChainLength, 18 Charge, 24 ChargeDensity, 96 charges.bin, 66 Choleskii, 36 Circle, 85 COEF.DAT, 68 Coefficients, 101 Constraints, 12, 14, 15 Contact, 75 Contact{}, 76 ContactHamiltonian, 77 ContactId, 77 ContactSeparation, 77 ContactTemperature, 88 ContactVector, 104 ContourPoints, 79 ConvergentForcesOnly, 12, 14–16, 24 CoordsAndCharges, 43 Coupling, 20 CouplingStrength, 18 CovalentRadius, 46 CustomisedHubbards, 24 CustomisedOccupations, 24, 28 Cutoff, 48, 101 CutoffCheck, 83 CutoffCN, 48 CutoffReduction, 52 Cylindrical, 87 Delta, 16, 53, 79 Dense, 58 detailed.out, 63 DetailedXML, 96 Device, 75 Device{}, 75 DFTB+U, 42 DFTB3, 49 DiagonalRescaling, 30 Differentiation, 24 Direction, 44 Directions, 59 DirectRead{}, 43 Dispersion, 24 DisplayModes, 93 Driver, 10 DynMixingParameters, 30 eigenvec.bin, 66 137 138 eigenvec.out, 66 EigenvecBin, 96 EigenvectorsAsTxt, 56 ElectricField, 24 Electrostatic Gates, 87 ElectrostaticPotential, 56 Electrostatics, 73 ELPA, 36 EnclosedPoles, 79 EnergyRange, 88 EnergyStep, 88 EnergyWindow, 60 ESP.dat, 67 EwaldParameter, 24 EwaldTolerance, 24 EXC.DAT, 69 ExcitedState, 10, 52 ExcitedStateForces, 60 Exponents, 101 External, 43, 44 FermiCutoff, 79, 82 FermiLevel, 76, 79 File, 23 FillBoxWithAtoms, 96 Filling, 24 FirstLayerAtoms, 75, 79 FixAngles, 12, 14, 15 FixedFermiLevel, 38 FixLengths, 12 FoldAtomsToUnitCell, 96 ForceEvaluation, 24 Frequency, 44 g{}, 19 Gate, 83, 87 GateDistance, 87 GateLenth_l, 87 GateLenth_t, 87 GatePotential, 87 GateRadius, 87 GaussianBlurWidth, 43 Generations, 15, 30, 31 Geometry, 10, 93, 103 Global, 85 GreensFunction, 79 GreensFunction{}, 73 Grid, 58 GridPoints, 59 GroundState, 96 INDEX Groups, 62 H5Scaling, 50 Hamiltonian, 10 hamreal.dat, 65 hamsqr.dat, 65 HardCutOff, 29 HCorrection, 24 Hessian, 12, 16, 94 Hessian, 93 HHRepulsion, 48 Host, 23 HubbardDerivs, 24 HybridDependentPol{}, 46 HybridPolarisations, 46 i-PI{}, 23 Id, 76, 104 IgnoreUnprocessedNodes, 61, 102 ImagComponent, 96 IndependentKFilling, 38 InitialCharges, 24 InitialSpins, 32, 33 InitialTemperature, 17 InitMixingParameter, 30, 31 IntegrationSteps, 20 IntegratorSteps, 18 InverseJacobiWeight, 30 Isotropic, 12, 14, 15, 20 KeepStationary, 16 KPointsAndWeights, 24 Label, 56 LatticeOpt, 12, 14, 15 LatticeVectors, 10 LBFGS{}, 15 LevelSpacing, 76 List of atoms, 12 LocalCurrents, 79 Localise, 56 localOrbs.bin, 57 localOrbs.out, 57 LowerCaseTypeName, 39 LowestEnergy, 79 Mass, 22 Masses, 16 MassPerAtom, 22 MaxAngularMomentum, 24 MaxAtomStep, 12, 14 INDEX MaxForceComponent, 12, 14, 15 MaximalWeight, 30 MaxIterations, 57 MaxLatticeStep, 12, 14, 15 MaxParallelNodes, 83, 88 MaxPoissonIterations, 83 MaxSCCIterations, 24 MaxSccIterations, 21, 42 MaxSKCutoff, 105 MaxSteps, 12, 14, 15, 23 md.out, 67 MDRestartFrequency, 16, 67 Memory, 15 MinEdgeLength, 99 MinimalGrid, 83 MinimalWeight, 30 MinimiseMemoryUsage, 54 MinSccIterations, 21 Mixer, 24 MixingParameter, 30 Monkhorst-Pack scheme, 41 MovedAtoms, 12, 14–16 MullikenAnalysis, 56 muPoints, 37 nInterationsELPA, 36 No, 29 NPH ensemble, 19 NPT ensemble, 19 NrOfCachedGrids, 96 NrOfExcitations, 60 NrOfPoints, 96 NTPoly, 36 NumericalNorm, 83 NVE ensemble, 17 NVT ensemble, 17, 18 Occupation, 101 OldSKInterpolation, 24, 29 OMM, 36 OnSiteCorrection, 24 Options, 10, 96 Orbital, 101 OrbitalPotential, 24 OrbitalResolved, 56 Order, 18, 39 Origin, 59, 99 OscillatorWindow, 60 OutputFile, 58 OutputPrefix, 12, 14–16 139 overreal.dat, 65 OverrideBulkBC, 83 OverrideDefaultBC, 83 oversqr.dat, 65 Parallel, 10 Parallelisations, 87 ParserOptions, 10, 96 ParserVersion, 61, 107 Periodic, 10 PEXSI, 36 Phase, 44 PipekMezey, 57 Planar, 87 PlotModes, 94 PlottedKPoints, 96 PlottedLevels, 96 PlottedRegion, 96 PlottedSpins, 96 PLShiftTolerance, 76 PointCharges, 43 Points, 58 Poisson{}, 73 PoissonAccuracy, 83 PoissonBox, 83 PoissonThickness, 83 Poles, 37 PolynomialRepulsive, 24 Port, 23 Potential, 76, 85 Prefix, 23, 39 Pressure, 12, 14, 15, 20 PreSteps, 20 ProcsPerPole, 37 ProjectStates, 56 Protocol, 23 PurificationMethod, 37 RandomSeed, 54 RangeSeparated, 24, 52 ReadChargesAsText, 54 ReadInitialCharges, 24 ReadOldBulkPotential, 83 ReadSurfaceGFs, 79, 80 RealAxisPoints, 79, 80, 82 RealAxisStep, 79, 80, 82 RealComponent, 96 RecomputeAfterDensity, 83 Region, 88 RelaxTotalSpin, 32 140 RemoveRotation, 93 RemoveTranslation, 93 RepeatBox, 96 ReselectIndividually, 18 ReselectProbability, 18 Resolution, 100 Restart, 18 RestartFrequency, 54 results.tag, 64 RScaling, 50 SavePotential, 83 SaveSurfaceGFs, 79, 80 Scale, 21 SCC, 17, 24 SCCTolerance, 24 SccTolerance, 21 Screening, 52 sec:Damp X-H, 49 sec:DFTB3-D3H5, 50 SelectedShells, 25 Separator, 39 ShellResolved, 56 ShellResolvedSCC, 24 ShellResolvedSpin, 24, 34 ShiftGrid, 96 ShowFoldedCoords, 54 SkipChargeTest, 54 SKMaxDistance, 29 SlaterKosterFiles, 24, 93, 104 Softening, 58 Solver, 24, 73 Spacing, 59 Sparse, 36, 37 SparseTolerances, 58 SpecifiedPLs, 104 SpectralRadius, 37 SpinConstants, 24, 60 SpinOrbit, 24 SpinPerAtom, 32, 33 SpinPolarisation, 24 SPX.DAT, 69 Square, 85 state resolved Mulliken population, 66 StateOfInterest, 60 Steps, 16 StepSize, 12 StopAfterParsing, 61, 102 Strength, 44 Suffix, 39 INDEX SymbolicFactorProcs, 37 Symmetry, 60 Task, 75–77 Task = ContactHamiltonian{}, 76 Task = UploadContacts, 78 Task = UploadContacts{}, 78 TDP.DAT, 70 Temperature, 18, 38, 76 TemperatureProfile{}, 16, 18, 19 TerminalCurrents, 88 TEST_ARPACK.DAT, 70 TestArnoldi, 60 Thermostat, 16 ThirdOrder, 24 ThirdOrderFull, 24 Threebody, 48 Threshold, 52 Timescale, 18, 20 TimeStep, 16, 44 TimingVerbosity, 54 Tolerance, 36, 37, 57 TotalAtomicDensity, 96 TotalChargeDensity, 96 TotalChargeDifference, 96 TotalSpinPolarisation, 96 TotalStateCoeffs, 60 TRA.DAT, 70 TransientSteps, 21 Transport, 75, 80, 103 Transport{}, 73 TransportOnly, 73 TruncateSKRange, 24, 104 Truncation, 37 TunnelingAndDos, 88 TypeNames, 10 TypesAndCoordinates, 10 UnpairedElectrons, 32 UseFromStart, 31 UseOmpThreads, 62 v{}, 19 Velocities, 16 Verbose, 96 Verbosity, 23, 83 WeightFactor, 30 WideBand, 76 WriteAutotestTag, 54 WriteBandOut, 56 INDEX WriteChargesAsText, 54 WriteCoefficients, 60 WriteDetailedOut, 54 WriteDetailedXML, 54 WriteEigenvectors, 56, 60 WriteHS, 54 WriteHSDInput, 61, 93 WriteLDOS, 88 WriteMulliken, 60 WriteRealHS, 54 WriteResultsTag, 54 WriteSPTransitions, 60 WriteStatusArnoldi, 60 WriteTransitionDipole, 60 WriteTransitions, 60 WriteTunn, 88 WriteXMLInput, 61, 93 WriteXplusY, 60 WScaling, 50 x{}, 19 XCH.DAT, 70 Xlbomd, 16 XlbomdFast, 16 XMakeMol, 94 XplusY.DAT, 70 XREST.DAT, 71 Yes, 29 141