DCDFTBMD USER MANUAL
Version 1.0

Nakai Group
Department of Chemistry and Biochemistry, Waseda University

November 1, 2018
Abstract
DCDFTBMD allows the users to perform quantum mechanical molecular
dynamics (MD) simulation for large systems by using divide-and-conquer(DC) density functional tight-binding (DFTB) method. The scaling of DCDFTB with respect to system size is close to linear. A high parallel efficiency
has been achieved in the implementation with the hybrid of MPI and
OpenMP. The program has functionalities not only for MD but also forgeometry optimization, vibrational frequency calculation, and property
evaluation.

1

Contents
1 Input for DCDFTBMD……………………………………….3
1.1 Keyword section………………………………………………….4
1.1.1 SCC…………………………………………………………………….5
1.1.2 DC…………………………………………………………………….11
1.1.3 OPT/OPTIMIZE……………………………………………………...14
1.1.4
1.1.5
1.1.6
1.1.7
1.1.8
1.1.9

MD……………………………………………………………………16NCC…………………………………………………………………..23
FREQ/FREQUENCY………………………………………………...24
DISP/DISPERSION………………………………………………….26
PBC…………………………………………………………………...32
MISC………………………………………………………………….33

1.2 Title section……………………………………………………...36
1.3 Parameter section………………………………………………..37
1.3.1 General overview……………………………………………………..37
1.3.2 Specification of additional constants…………………………………37

1.4 Geometry section………………………………………………..40
1.4.1 General overview……………………………………………………..40
1.4.2 Specification of additional information………………………………41
1.4.3 Manual construction of subsystems in DC-DFTB…………………...42

2 Running DCDFTBMD……………………………………...44
3 Output of DCDFTBMD…………………………………….45
3.1 Standard output file……………………………………………..45
3.2 Detailed data file………………………………………………..48
3.3 Other output files……………………………………………….53

References……………………………………………………...55
2

1 Input for DCDFTBMD
To run the program DCDFTBMD, you need to prepare the input file and parameters(Slater-Koster files). The input file, which is named “dftb.inp” by default, mainly
consists of 4 parts below:


Information about the calculation (Keyword section)



Brief description of the calculation (Title section)

Specification of Slater-Koster parameters (Parameter section)

Specification of molecular coordinates (Geometry section)

Each section is terminated with blank line (starting new line is indicated as  in thisdocument). For example, the input file looks like as follows in the case of SCC-DFTB
single point energy calculation of water molecule:
SCC=TRUE˽DC=FALSEKeyword section
TITLE
2
O˽2
oo.spl˽oh.spl
H˽1
ho.spl˽hh.spl


Title section
Parameter section

3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000


Geometry section

Additional input may be required after geometry section if subsystems are specified
manually in DC approach (see section 1.4.3). Note that blank area (represented as ˽ inthis document) should be written with spaces (DO NOT USE TABS!).

3

1.1 Keyword section
SCC=TRUE˽DC=FALSESCC and DC keywords areturned on and off

The present example specifies two keywords in a single line. The same specification ispossible by giving two separated lines as follows:
SCC keyword is turned on
DC keyword is turned off

SCC=TRUE
DC=FALSE

The following subsections account for the available keywords at present. Two keywords
SCC and DC are activated if no keywords are given in the input. Other keywords can beturned on using the following forms:
keyword=TRUE
The above specification employs the default value of options in curly bracket. Several
keywords have options for controlling the calculation. The italic type in parenthesis
stands for the data type of option. The options can be specified by the following forms:


keyword=(option1=value1˽option2=value2˽…)



keyword=(option1=value1)
keyword=(option2=value2˽…)

If option is written in red or blue character, that option is available in only serial and
OpenMP or MPI/OpenMP version, respectively.

4

1.1.1 SCC
The keyword SCC sets whether a self consistent charge (SCC) calculation is made or not.
As already mentioned, the default is SCC=TRUE, namely the SCC-DFTB (DFTB2) [1]
calculation is performed.


MAXITER (integer) {200}Maximum number of SCC cycles to reach convergence.



BROYITER/MAXBROYDEN (integer) {70}Maximum number of charge mixing using Broyden method.



ENERGYCONV/ECONV (real) {1.0e−09}
Convergence criteria for the energy change between two SCC cycles in atomic
unit.



DENSITYCONV/DCONV (real) {1.0e−06}
Convergence criteria for the maximum difference in any charge density between
two SCC cycles in atomic unit.



READCHARGE (method) {NONE}
Enables to read initial charge information from file.
 NONE constructs initial Mulliken population based on the number ofvalence electrons.
 BINARY reads Mulliken population from the binary file “chrgfile”. If thisoption is set but “chrgfile” is not found in the same directory as the one
including the input file, this option will be turned off.
 ASCII reads net atomic charges from the ASCII file “charge.dat”. If thisoption is set but “charge.dat” is not found in the same directory as the one
including the input file, this option will be turned off.
Examples for water molecule where oxygen corresponds to atom label 1 aregiven below.
For serial and OpenMP versions, the first line of “charge.dat” is either“ALLATOM” or “ATOMSPECIES”.
 ALLATOM lists the net charge for every atom from the second line.
ALLATOM
5

−0.6

Mulliken population of O
(atom label 1) is 6.6
Mulliken population of H
(atom label 2) is 0.7
Mulliken population of H
(atom label 3) is 0.7

−0.3
−0.3

 ATOMSPECIES lists the chemical symbol and the net charge in thefirst and second columns from the second line until detecting blank line
if any.
ATOMSPECIES
O˽−0.6

Mulliken population of O
(atom label 1) is 6.6
Mulliken population of H
(atom labels 2 and 3) is 0.7

H˽−0.3

For MPI/OpenMP version, the first line of “charge.dat” is either“ALLSHELL” or “SHELLSPECIES”.
 ALLSHELL lists the net charge for every shell from the second line.
ALLSHELL
−0.2

Mulliken population of sshell of O (atom label 1) is
1.8
Mulliken population of p

−0.8

shell of O (atom label 1) is
4.8
Mulliken population of sshell of H (atom label 2) is
0.7
Mulliken population of sshell of H (atom label 3) is
0.7

−0.3

−0.3

 SHELLSPECIES lists the chemical symbol, subshell label (s/p/d), and6

the net charge in the first, second, and third columns from the secondline until detecting blank line if any. If system is open-shell, alpha/beta
spin (a/b) needs to be specified in the fourth column.
SHELLSPECIESO˽s˽−0.2

Mulliken population of sshell of O (atom label 1) is
1.8
Mulliken population of p
shell of O (atom label 1) is
4.8
Mulliken population of sshell of H (atom labels 2
and 3) is 0.7O˽p˽−0.8H˽s˽−0.3



DAMPXH (logical) {FALSE}
Enables the use of damped γ-function for X–H pairs. [2–4]

DAMPXHZETA (real) {4.0e+00}
Empirical parameter for the damped γ-function for X–H pairs.



THIRDDIAG (logical) {FALSE}
Enables the use of third order corrections only for diagonal terms. [3]Hubbard derivatives need to be specified in Parameter section (see section 1.3.2).



THIRDFULL (logical) {FALSE}
Enables the use of full third order corrections. [4]Hubbard derivatives need to be specified in Parameter section (see section 1.3.2).
The so-called DFTB3 [4] requires DAMPXH=TRUE.

SPIN (logical) {FALSE}
Enables to perform collinear spin polarized calculation. [5, 6]
Spin constants have to be stored in different file for each atom (see section 1.3.2).



ALMIX (real) {0.2e+00}
Charge mixing parameter used in charge mixing method (0 ≤ ALMIX ≤ 1).

MIXER (integer) {1}
Charge mixing method to accelerate SCC convergence.
7

Note that modified Broyden mixing is used to solve coupled perturbed DFTB
equation in analytical Hessian calculation.
 1: Modified Broyden [7] 2: Simple mixing 3: Modified Anderson [8] 4: DIIS [9]

NGEN (integer) {4}
Number of iteration to start Anderson and DIIS charge mixing.



SOLVER (integer) {1}
Eigensolver for diagonalization of Hamiltonian matrix.
 1: DSYGVD
 2: DSYGV
 3: DM_VGEVPH
This may be available if the platform is the K computer, Fujitsu FX10, andFujitsu FX100.

EFIELD (logical) {FALSE}
Enables to include external static electric field.
This is available for only cluster models.
Static electric field vector is set after atomic coordinates starting from thecharacter “EF” (see section 1.4.1).



POLAR (logical) {FALSE}
Enables to calculate dipole moment and static polarizability using external static
electric field.
This is available for only cluster models.



POLEFSTR (real) {2.0e−03}
Strength of external static electric field used in dipole moment and static
polarizability calculation.



POINTCHARGE (logical) {FALSE}
Enables to include point charges.
Coordinates and charge value are set after atomic coordinates starting from thecharacter “PC” (see section 1.4.1).



NUMPC (integer) {0}Number of point charges in system.
8



MAYER (logical) {FALSE}
Enables to perform Mayer bond order analysis. [10]

CM3 (logical) {FALSE}
Enables to perform CM3 charge analysis. [11, 12]

ORBITALPOP/ORBPOP (logical) {FALSE}
Enables to print Mulliken population of each orbital.



ZEROCHARGE (logical) {FALSE}
Enables to permanently set values of neutral atoms as initial charge guess in SCC
calculation during geometry optimization and molecular dynamics.
Deactivation uses the converged charges of previous geometry step as initial
guess (better than zero charge guess).



PSEUDODIAG/PDIAG (logical) {FALSE}
Enables to perform psuedodiagonalization using Jacobi rotations. [13]After reaching the SCC convergence with psuedodiagonalization, full
diagonalization is performed to obtain final eigenvalues and eigenvectors.



PDIAGTHRESH (real) {1.0e−02}
Threshold to activate psuedodiagonalization.
Full diagonalization is performed to update reference eigenvalues until charge
density difference becomes less than the specified value.



PDIAGKAPPA (real) {4.0e−02}
Accuracy control parameter for psuedodiagonalization.



GAMMASPLINE (logical) {FALSE}
Enables to use cubic spline interpolation for exponentially decaying short-range
part of γ-function.
For MPI/OpenMP version, this is available in only RUNTYPE=ATOM.
See reference [14] for theory, implementation, and numerical assessments of thepresent option.



GAMMASPLINEGRID (real) {4.0e−02}
Grid size of cubic spline interpolation for exponentially decaying short-range
part of γ-function in atomic unit.



UNPAIREDELEC (real) {0.0e+00}
Number of unpaired electrons.
9

It is possible to set non-integer positive value.
Note that the specified spin multiplicity is discarded except for zero unpaired
electrons.


SPINRELAX (logical) {FALSE}
Enables to relax total spin in spin-polarized calculation.



NELECTHRESH (real) {1.0e−06}
Threshold for conservation of total number of electrons in SPINRELAX=TRUE.

RESTRICT (logical) {FALSE}
Enables to use restricted approach in open-shell doublet calculation.
Default is unrestricted approach.



UDERIVGAUSS (logical) {FALSE}
Enables to add charge dependent Gaussian function to Hubbard derivative in
THIRDDIAG=TRUE [15].



UDERIVGAUSSD0 (real) {−9.0e−02}
Parameter of Gaussian function (height) used in UDERIVGAUSS=TRUE.

UDERIVGAUSSG0 (real) {1.61e+01}
Parameter of Gaussian function (width) used in UDERIVGAUSS=TRUE.

UDERIVGAUSSQ0 (real) {7.5e−01}
Parameter of Gaussian function
UDERIVGAUSS=TRUE.

(positionofcenter)usedin

RUNTYPE (method) {ATOM}
Selects how to resolve DFTB calculations.
Note that serial and OpenMP versions employ atom resolved calculation.
 SHELL0: Based on valence shells. This is incompatible with
THIRDFULL=TRUE and THIRDDIAG=TRUE.
 SHELL: Based on valence shells. This is incompatible with
THIRDDIAG=TRUE.
 ATOM: Based on atoms. This is incompatible with ORSCC=TRUE.

ORSCC (logical) {FALSE}
Enables the use of orbitally resolved SCC calculation.
Different Hubbard parameters for different atomic angular momentum shells areused instead of the use of s-shell only.
10The DFTB2 [6] and DFTB3 [16] methods can be extended with this treatment.


GSDOMAIN (logical) {FALSE}
Enables the use of cell-linked list (CLL) algorithm for accelerating short-range
part of charge-charge interaction calculation in GAMMASPLINE=TRUE.

GSDOMAINBLOCKXYZ/GSDOMAINBLOCK (integer) {5}
Block size in each of x-, y-, and z-directions used in GSDOMAIN=TRUE.

GSDOMAINBLOCKX (integer) {5}
Block size in x-direction used in GSDOMAIN=TRUE.

GSDOMAINBLOCKY (integer) {5}
Block size in y-direction used in GSDOMAIN=TRUE.

GSDOMAINBLOCKZ (integer) {5}
Block size in z-direction used in GSDOMAIN=TRUE.

GSDOMAINNEIGHBOR (integer) {1}
Number of linked neighbors referred to as adjacent cells in GSDOMAIN=TRUE.

MIXERSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for charge mixing methods. This is
implemented for atom resolved calculations (RUNTYPE=ATOM) and for
closed-shell singlet calculations.

1.1.2 DCThe keyword DC sets whether a divide-and-conquer approach [17, 18] is employed ornot. This keyword is turned on without any specification as well as SCC keyword.
DC=FALSE is necessary if DC-DFTB calculation is deactivated.


BETA (real) {8.0e+02}
Inverse temperature parameter for the Fermi function in atomic unit.



BUFRAD/BUFFERRADIUS (real) {5.0e+00}
Buffer radius to construct buffer region in Angstrom (Å) unit.
Atoms included in spheres centered at atoms in the central region are selected.



SUBTYPE (method) {AUTO}
Method to construct central region.
 AUTO constructs subsystems automatically by dividing total system with
11







the specified cubic or rectangular grid.
Default is 5.0 Å cubic box.
SPECIFY defines subsystems one by one for each atom (see section 1.4.3).
ATOM defines one atom as each subsystem.
SYSTEM defines all atoms as one subsystem.
SEMIAUTO defines subsystems as the hybrid of AUTO and SPECIFY (seesection 1.4.3).
AUTOXH constructs subsystems similarly to AUTO without cutting

covalent bonds involving hydrogen atom.
 SEMIAUTOXH constructs subsystems similarly to SEMIAUTO without
cutting covalent bonds involving hydrogen atom.


DELTARXYZ/DELTAR (real) {5.0e+00}
Grid length of cubic grid used in SUBTYPE=AUTO in Å unit.



DELTARX (real) {5.0e+00}
Grid length of x-component used in SUBTYPE=AUTO in Å unit.



DELTARY (real) {5.0e+00}
Grid length of y-component used in SUBTYPE=AUTO in Å unit.



DELTARZ (real) {5.0e+00}
Grid length of z-component used in SUBTYPE=AUTO in Å unit.

NOFERMI (logical) {FALSE}
Common Fermi level is not determined during DC calculation.



TRANSCOORD (logical) {TRUE}
Transformation of Cartesian coordinate in SUBTYPE=AUTO.
 TRUE: Standard orientation is used.
 FALSE: Initial geometry is used as it is.



SUBXHSCALE (real) {1.1e+00}
Scaling factor for the sum of covalent radii [19] to judge the presence of covalentbonds involving hydrogen atom in SUBTYPE=AUTOXH and
SUBTYPE=SEMIAUTOXH.

FERMITYPE (method) {PARALLEL}
Algorithm to determine Fermi level. See reference [20] for details.
 PARALLEL: CRI (compute-reduce-interpolate) algorithm is used.
12 SERIAL: GSS (gather-sort-solve) algorithm is used.


DOMAIN (logical) {FALSE}
Enables the use of CLL algorithm for accelerating assignment of subsystems tolocalization regions.



DOMAINBLOCKXYZ/DOMAINBLOCK (integer) {5}
Block size in each of x-, y-, and z-directions used in DOMAIN=TRUE.

DOMAINBLOCKX (integer) {5}
Block size in x-direction used in DOMAIN=TRUE.

DOMAINBLOCKY (integer) {5}
Block size in y-direction used in DOMAIN=TRUE.

DOMAINBLOCKZ (integer) {5}
Block size in z-direction used in DOMAIN=TRUE.

DOMAININDEXTYPE (method) {ATOM}
Definition of particle for the index of linked list in DOMAIN=TRUE.
 SUBSYSTEM: The index of subsystem will be given if all constituent atoms
are located in the same domain.
 ATOM: The index will be assigned to all atoms.



DOMAINNEIGHBOR (integer) {1}
Number of linked neighbors referred to as adjacent cells in DOMAIN=TRUE.

SINGLEPREC (logical) {FALSE}
Enables the use of single precision floating point number for distancecalculations in assignment of subsystems to localization regions.



XHDOMAIN (logical) {FALSE}
Enables the use of CLL algorithm for judging the presence of covalent bonds
involving hydrogen atom in SUBTYPE=AUTOXH.

XHDOMAINBLOCKXYZ/XHDOMAINBLOCK (integer) {5}
Block size in each of x-, y-, and z-directions used in XHDOMAIN=TRUE.

XHDOMAINBLOCKX (integer) {5}
Block size in x-direction used in XHDOMAIN =TRUE.

XHDOMAINBLOCKY (integer) {5}13

Block size in y-direction used in XHDOMAIN =TRUE.
XHDOMAINBLOCKZ (integer) {5}
Block size in z-direction used in XHDOMAIN =TRUE.

XHDOMAINNEIGHBOR (integer) {1}
Number of linked neighbors referred to as adjacent cells in XHDOMAIN=TRUE.

1.1.3 OPT/OPTIMIZE
The keyword OPT sets whether the geometry optimization is performed or not.


MAXITER (integer) {50}Maximum number of optimization steps.



GRADCONV/GCONV (real) {1.0e−04}
Convergence criteria for the maximum absolute value of force component inatomic unit.



OPTTYPE (integer) {1}
Method for geometry optimization.
 1: BFGS
This uses covalent radius of atoms to construct initial Hessian matrix, which
are taken from reference [19].
 2: Steepest descent 3: Conjugate gradientSTEPSIZE (real) {1.0e+02}
Step size along the forces for steepest descent method in atomic unit.

CONSTRAINT (logical) {FALSE}
Enables to specify geometry constraints.
Constraint vector of x, y, z components is necessary for atoms to be constrained(see section 1.4.2).



LATTICEOPT (logical) {FALSE}
Enables to optimize translation vectors.
It is not allowed to impose partial geometry constraints and to set external point
charges.
14Note that conjugate gradient optimizer is used when OPTTYPE=1 (BFGS).
EXTPRESSURE (real) {0.0e+00}
External pressure during optimization of translation vectors in atomic unit.
If value is not zero, enthalpy (E + pV) of the system is additionally printed.



ISOSCALING (logical) {FALSE}
Enables to scale translation vectors uniformly during optimization of translation
vectors.
This is incompatible with FIXANGLE=TRUE.

FIXANGLE (logical) {FALSE}
Enables to fix angles between translation vectors during optimization oftranslation vectors.
This is incompatible with ISOSCALING=TRUE.

FIXANGLE1 (logical) {FALSE}
Enables to fix length of first translation vector during optimization of translation
vectors with FIXANGLE=TRUE.

FIXANGLE2 (logical) {FALSE}
Enables to fix length of second translation vector during optimization oftranslation vectors with FIXANGLE=TRUE.

FIXANGLE3 (logical) {FALSE}
Enables to fix length of third translation vector during optimization of translation
vectors with FIXANGLE=TRUE.

LATTICEONLY (logical) {FALSE}
Enables to fix atomic coordinates during optimization of translation vectors.
It is possible to apply the same option by putting a single asterisk (*) after zcoordinate for all atoms (see section 1.4.2).



STEPCGRATOM (real) {1.0e−01}
Maximum step size of atom relaxation for conjugate gradient method in atomic
unit.



STEPCGRLATTICE (real) {1.0e−01}
Maximum step size of lattice relaxation for conjugate gradient method in atomic
unit.

15



BFGSSOLVER (integer) {1}
Eigensolver for diagonalization of Hessian matrix in BFGS method.
 1: DSYEVD
 2: DSYEVPRINTCOORD (logical) {TRUE}
Enables to print molecular coordinate at each geometry step in the standard
output file.
Deactivation (PRINTCOORD=FALSE) corresponds to print only the final
molecular coordinate in the standard output file.



PRINTGRAD (logical) {TRUE}
Enables to print DFTB energy gradient at each geometry step in the standard
output file.
In the case of PRINTGRAD=FALSE, maximum and root mean square forces
are written in the standard output file.

1.1.4 MDThe keyword MD sets whether the molecular dynamics (MD) using velocity Verlet
integration scheme is performed or not.


NSTEP/MAXSTEP (integer) {10000}
Number of MD steps.



DELTAT/DELTASTEP (real) {1.0e−15}
Time interval between two MD steps in second (s) unit.



BATHTEMP/TEMPERATURE (real) {2.9815e+02}
Target temperature of thermostat in Kelvin (K) unit.

NVE (logical) {TRUE}
Enables to run MD with NVE ensemble.
If deactivated, NVT ensemble will be used.

NVT (logical) {FALSE}
Enables to run MD with NVT ensemble.
If deactivated, NVE ensemble will be used.



NVTTYPE (integer) {3}
Method to control temperature for achieving NVT or NPT ensemble.
16 1: Velocity scaling [21]
 2: Nosé-Hoover chain [22]
 3: Nosé-Hoover chain [22]In comparison to NVTTYPE=2, NVTTYPE=3 allows to control the
coupling strength.
 4: Berendsen [23]
 5: Andersen [24]

RESTART (logical) {FALSE}
Enables to read information for restarting MD from the binary file “restart”.



SOFT/SOFTPOT (logical) {FALSE}
Enables to use soft potential.



SOFTRANGE (real) {1.0e+01}
Radius of soft potential in Å unit.



SOFTK (real) {1.5e+00}
Hardness of soft potential in kcal/(mol∙Å2) unit.



INITTEMP (real) {2.9815e+02}
Target temperature of initial velocities in K unit.
Initial velocities are generated following the Maxwell-Boltzmann distribution.



SCALESTEP (integer) {20}Interval to scale velocities for velocity scaling method (NVTTYPE=1).



ERRORTEMP (real) {1.0e+02}
Threshold of temperature for scaling velocities in K unit.
If the absolute value of difference between the target temperature of thermostat
and the MD temperature is greater than ERRORTEMP, the velocities are scaled
so that the temperature becomes BATHTEMP for NVT methods except forvelocity scaling (NVTTYPE=2–5).



READVELOCITY (logical) {FALSE}
Enables to read velocities for all atoms from the ASCII file.
 TRUE: Read the file “veloc.dat” in Å/fs unit (fs is femtosecond). The ASCII
file should have 3 columns (x, y, z velocity components) for each atom of
the given structure in Geometry section.
 FALSE: No file to be read.
17

NOSECLEAN (logical) {FALSE}
Enables to initialize internal parameters of Nosé-Hoover chain thermostat when
MD is restarted.



PRINT (integer) {1}Interval to write coordinate, velocity, and Mulliken atomic charges to the files“traject”, “velocity”, and “mulliken”, respectively.



NBATH (integer) {3}
Number of particles for Nosé-Hoover chain thermostat (NVTTYPE=3).



COUPLINGSTR (real) {1.5e+03}
Frequency of oscillation of the thermostat particles in cm−1 unit for Nosé-Hoover
chain thermostat (NVTTYPE=3).



SEEDTYPE (integer) {1}
Method of random number generation for initial velocities.
 1: No specification (the same initial velocities may be reproduced on thesame platform for the same given structure in Geometry section)
 2: Generation of random number from the current time.
 3: Specification of arbitrary random number using RANDOMSEED option.

RANDOMSEED (integer) {0}Seed to generate random number for initial velocities using SEEDTYPE=3.



TAUTEMP (real) {1.0e−13}
Time constant for Berendsen thermostat (NVTTYPE=4) in s unit.



COLLISIONFREQ (real) {1.0e−02}
Random collision frequency for Andersen thermostat (NVTTYPE=5) in fs−1 unit.



MASSIVE (logical) {FALSE}
Method of scaling velocities for Andersen thermostat (NVTTYPE=5).
 TRUE: All velocities are scaled simultaneously with the specified
probability.
 FALSE: Each atomic velocity is scaled individually with the specified
probability.



LIMULLIKEN (logical) {FALSE}
Enables to predict better initial Mulliken charges using Lagrange interpolationto accelerate SCC convergence.
18



LIDEGREE (integer) {0}Degree of Lagrange polynomial in LIMULLIKEN=TRUE.
The default setting 0 employs the previous converged charges for initial guess(conventional approach).



LIMINDEG (logical) {FALSE}
Enables to obtain the appropriate degree of Lagrange polynomial within a range
from 0 to the value set in LIDEGREE in LIMULLIKEN=TRUE.
Degree which gives the smallest absolute error of the predicted charges from theconverged values at N step is used at N + 1 step.



RMROTATION (logical) {FALSE}
Removal of center of mass rotation in addition to translation.
This option is available in only cluster model simulations.



RMSTEP (integer) {10000}
Interval to remove center of mass translation
RMROTATION=TRUE) for NVTTYPE=1–4.



(and

rotation

if

RATTLE (method) {NONE}
Constraint MD simulation using RATTLE algorithm [25] for bond lengths.
 NONE corresponds to no constraints of bond lengths.
 COVALENTXH fixes X–H covalent bonds as the sum of their covalent radii,
which are taken from reference [19].
 INITGEOMXH fixes X–H covalent bond distance as the value of initialstructure.
 SPECIFY sets any bond constraints using the file “rattle.dat”.
The first and second columns specify atom labels whose distance is fixed.
The third column is an upper-case letter to determine how the specified bonddistance is fixed.
“C” fixes to the sum of
covalent radii“I” fixes to the value of
initial structure“S” fixes to the specified
value (0.95 Å in thisexample)1˽2˽C
2˽3˽I
1˽3˽S˽0.95

19Note that the sum of
covalent radii is used if nofourth column


RATTLEXH (real) {1.2e+00}
Threshold of X–H pair as the covalent bond in RATTLE=COVALENTXH andRATTLE=INITGEOMXH in Å unit.



RATTLEITER (integer) {100}Maximum number of cycles to satisfy RATTLE position and velocity constraints.



RATTLECONV (real) {1.0e−06}
Convergence criteria for the RATTLE constraints in atomic unit.



ANNEAL (logical) {FALSE}
Enables to perform simulated annealing simulation.
Temperature scheduling is specified using the file “anneal.dat” up to 99 profiles.
The first step is set toC˽1001˽0300.0

constant value of 300 K (C)INITTEMP is updated withthe specified value
Next 1000 steps increase
from 300 K to 1000 K
linearly (L)Next 1000 steps keep thetemperature
1000
K(equilibration)

L˽1000˽1000.0

C˽1000˽1000.0

Next 1000 steps decrease
from 1000 K to 500 K
exponentially (E)The fourth column controls
the shape of exponentialfunction
Note that the value of 5.0 isused if no fourth column
Next 1000 steps keep theE˽1000˽0500.0˽3.0

C˽1000˽0500.0

temperature 500 K
20

Next 1000 steps decrease
from 500 K to 300 K
sigmoidally (S)The fourth column controls
the shape of sigmoid
function
Note that the value of 3.0 isused if no fourth column

S˽1000˽0300.0˽2.0



SOFTCENTERTYPE (method) {ATOM} ATOM defines one of atoms as center of soft potential.
 SPECIFY manually sets the position of xyz coordinate.
 COC is the abbreviation of “center of coordinate”.
 COM is the abbreviation of “center of mass”.



SOFTATOMNO (integer) {1}
Atom label to be the center
SOFTCENTERTYPE=ATOM.

of

soft

potential

used

in



SOFTPOINTXYZ/SOFTPOINT (real) {0.0e+00}
Cartesian coordinate of the center of soft potential satisfying x = y = z used in
SOFTCENTERTYPE=SPECIFY in Å unit.



SOFTPOINTX (real) {0.0e+00}
x-componentof
thecenterofsoft
SOFTCENTERTYPE=SPECIFY in Å unit.





used

in

potential

used

in

potential

used

in

SOFTPOINTY (real) {0.0e+00}
y-componentof
thecenterofsoft
SOFTCENTERTYPE=SPECIFY in Å unit.



potential

SOFTPOINTZ (real) {0.0e+00}
z-componentof
thecenterofsoft
SOFTCENTERTYPE=SPECIFY in Å unit.
SOFTSHAPETYPE (method) {SPHERE}
Method to specify shape of soft potential.
 SPHERE is sphere.
 SQUARE is cubic or rectangular box.
21 ELLIPSOID is ellipsoid.


SOFTSIDEXYZ/SOFTSIDE (real) {1.0e+01}
Grid length of cubic grid used in SOFTSHAPETYPE=SQUARE in Å unit.



SOFTSIDEX (real) {1.0e+01}
Grid length of x-component used in SOFTSHAPETYPE=SQUARE in Å unit.



SOFTSIDEY (real) {1.0e+01}
Grid length of y-component used in SOFTSHAPETYPE=SQUARE in Å unit.



SOFTSIDEZ (real) {1.0e+01}
Grid length of z-component used in SOFTSHAPETYPE=SQUARE in Å unit.



SOFTSEMIAXISX (real) {1.0e+01}
Length of x-semi-principal axis used in SOFTSHAPETYPE=ELLIPSOID inÅ unit.



SOFTSEMIAXISY (real) {1.0e+01}
Length of y-semi-principal axis used in SOFTSHAPETYPE=ELLIPSOID inÅ unit.



SOFTSEMIAXISZ (real) {1.0e+01}
Length of z-semi-principal axis used in SOFTSHAPETYPE=ELLIPSOID inÅ unit.



NPH (logical) {FALSE}
Enables to run MD with NPH ensemble.
If deactivated, NVE ensemble will be used.



NPT (logical) {FALSE}
Enables to run MD with NPT ensemble.
If deactivated, NVE ensemble will be used.



NPHTYPE (integer) {1}
Method to control pressure for achieving NPH or NPT ensemble.
 1: Berendsen [23]

NPTTYPE (integer) {1}
Method to control pressure and temperature for achieving (isotropic) NPT
ensemble.
 1: Hoover-Andersen [22]22EXTPRESSURE (real) {1.01325e+05}
Target pressure during NPH or NPT ensemble simulation in Pascal (Pa) unit.



TAUPRESS (real) {2.0e−12}
Time constant for Berendsen barostat (NPHTYPE=1) in s unit.



ISOSCALING (logical) {FALSE}
Enables to scale translation vectors uniformly during NPH or NPT ensemble
simulation.
Default is anisotropic change of unit cell.



CALCPRESSURE (logical) {FALSE}
Enables to calculate virial pressure in NVE or NVT ensemble.
This option is supported for only periodic boundary condition.



COUPLINGSTRLAT (real) {1.5e+01}
Frequency of oscillation of the thermostat particles in cm−1 unit for HooverAndersen NPT ensemble (NPTTYPE=1).



CHECKPOINT (integer) {100}Interval to write the second restart file “restart_chk”.



RANDTYPE (integer) {1}
Method to generate random numbers.
 1: Fortran intrinsic random number generator
Different platform may generate different random numbers using the sameseed.
 2: Random number generator taken from reference [26].

1.1.5 NCC
The keyword NCC sets whether a non-SCC calculation (NCC-DFTB or DFTB1) [27, 28]
is performed or not. NCC=TRUE is identical to set SCC=FALSE.
PRINTMULLIKEN (logical) {FALSE}
Enables to perform Mulliken population analysis.
This option is required for further analysis of NCC Mulliken population.



MAYER (logical) {FALSE}
Enables to perform Mayer bond order analysis. [10]23Activation of PRINTMULLIKEN is required.


CM3 (logical) {FALSE}
Enables to perform CM3 charge analysis. [11, 12]
Activation of PRINTMULLIKEN is required.
Note that the fitted parameters are taken from DFTB2.

ORBITALPOP/ORBPOP (logical) {FALSE}
Enables to print Mulliken population of each orbital.
Activation of PRINTMULLIKEN is required.



SOLVER (integer) {1}
Eigensolver for diagonalization of Hamiltonian matrix.
 1: DSYGVD
 2: DSYGV
 3: DM_VGEVPH
This may be available if the platform is the K computer, Fujitsu FX10, andFujitsu FX100.

RESTRICT (logical) {FALSE}
Enables to use restricted approach in open-shell doublet calculation.
Default is unrestricted approach.

1.1.6 FREQ/FREQUENCY
The keyword FREQ sets whether the harmonic vibrational frequency analysis isperformed or not. If this keyword is turned on together with geometry optimization(OPT=TRUE), frequency calculation is carried out for the last optimization step of given
structures. The frequency calculation is currently available for cluster models.


FREQTYPE (integer) {1}
Method of Hessian matrix construction.
 1: Analytical second derivatives.
This
is
not
implementedfor
includingpointcharges(POINTCHARGE=TRUE in SCC keyword), for spline interpolation of γfunction (GAMMASPLINE=TRUE in SCC keyword), for open shell
systems, for including chemical-potential equalization correction(CPE=TRUE in SCC keyword), and for any shell resolved calculations in
MPI/OpenMP version (RUNTYPE=SHELL/SHELL0 in SCC keyword
24

and ORSCC=TRUE in SCC keyword).
 2: Numerical second derivatives using analytically calculated first
derivatives.


STEPSIZE (real) {1.0e−04}
Step size for numerical differentiation of analytically calculated first derivatives
in atomic unit.



PRINTHESS (logical) {FALSE}
Enables to print the lower triangle of Hessian matrix in the standard output file.
Additionally, Hessian matrix information is stored in the binary file“restart_hess” for restarting purpose.



THERMOTEMP (real) {2.9815e+02}
Target temperature for themochemical analysis in K unit.



THERMOPRES (real) {1.0e+00}
Target pressure for thermochemical analysis in standard atmosphere (atm) unit.



THERMOSCALE (real) {1.0e+00}
Frequency scaling factor for themochemical analysis.



PRECFORMAT (logical) {FALSE}
Enables to print high precision format for frequency eigenvectors.

READHESS (method) {NONE}
Enables to read Hessian matrix information from file.
 NONE calculates the Hessian matrix either analytically or numerically.
 BINARY reads Hessian matrix information from the binary file“restart_hess”.
 ASCII reads Hessian matrix information from the ASCII file “hess.dat”.
Each line consists of three columns: the first and second integers are used to
specify the index of Hessian matrix and the third value is an element ofHessian matrix.



READISOTOPE (logical) {FALSE}
Enables to specify arbitrary atomic mass.
The specification form (ISO=value) is necessary for atoms to be modified fromthe default mass in Slater-Koster parameter file (see section 1.4.2).



PROJECTION (logical) {TRUE}
25

Projection out translational and rotational modes from normal coordinate
analysis.
Deactivation (PROJECTION=FALSE) corresponds to the diagonalization offull Hessian matrix.


CPCONV (real) {1.0e−06}
Convergence criteria for coupled-perturbed DFTB equations in atomic unit.



EFSTR (real) {2.0e−03}
Strength of external static electric field used in infrared (IR) and Raman
intensities calculation.



INTERRUPTCP (logical) {FALSE}
Enables to interrupt solving coupled-perturbed DFTB equations.
The partially constructed Hessian matrix is stored in the ASCII file “hess.dat”
for restarting purpose.



CPSTART (integer) {1}
Atom index to specify the beginning of partial Hessian matrix to be constructedin INTERRUPTCP=TRUE.
CPEND (integer) {1}
Atom index to specify the end of partial Hessian matrix to be constructed in
INTERRUPTCP=TRUE.

HESSSOLVER (integer) {1}
Eigensolver for diagonalization of Hessian matrix.
 1: DSYEVD
 2: DSYEVSAVEMEMORYANA (logical) {FALSE}
Enables to reduce memory usage for computation of analytical Hessian in SCC
calculation in compensation for the increased computational cost.

1.1.7 DISP/DISPERSIONThe keyword DISP sets whether the dispersion interaction is empirically corrected or not.


DISPTYPE (integer) {1}
Method of dispersion correction.
 1: Grimme’s DFT-D2 [29]26 2: Slater-Kirkwood [30]
The required dispersion constants have to be stored in different file for eachelement (see section 1.3.2).
 3: Lennard-Jones [31]The default values of Lennard-Jones parameters are taken from values of
UFF force field [32].
It is also possible to specify Lennard-Jones parameters manually in different
file for each element (see section 1.3.2).
 4: Grimme’s DFT-D3 [33]This is so-called zero-damping.
 5: Grimme’s DFT-D3(BJ) [34]This is so-called Becke and Johnson damping.
 6: Řezáč’s D3H4 [35]Since this corrects hydrogen bonding interaction empirically, the damped γfunction may not be combined togerther (DAMPXH=FALSE in SCC
keyword). See, for example, reference [36].
 7: DFT-ulg [37]
 8: dDMC [38]
 9: Řezáč’s D3H5 [39]Since H5 correction is designated to modify the original γ-function,
DAMPXH=FALSE needs to be used in SCC keyword.
The correction is not available for non-SCC calculations (NCC=TRUE or
SCC=FALSE).
 10: Kubillus’s D3X [40]

UPDATESKPARAMS (logical) {FALSE}
Update Slater-Kirkwood dispersion parameters in each DFTB calculation.
Deactivation (UPDATESKPARAMS=FALSE) corresponds to keep dispersion
parameters fixed with those of initial geometry during the entire calculation.



READUFF (logical) {FALSE}
Enables to specify Lennard-Jones parameters manually for Lennard-Jones
dispersion (see section 1.3.2).



S6 (real) {1.0e+00}
Global scaling factor for Grimme’s DFT-D2.

ZEROS8 (real) {6.73e−01}27

Global scaling factor forterm of Grimme’s DFT-D3.



ZERORS6 (real) {1.235e+00}
Parameter for damping function of Grimme’s DFT-D3.



BJS8 (real) {5.883e−01}
Global scaling factor forterm of Grimme’s DFT-D3(BJ).



BJA1 (real) {5.719e−01}
Parameter for Becke and Johnson damping of Grimme’s DFT-D3(BJ).



BJA2 (real) {3.6017e+00}
Parameter for Becke and Johnson damping of Grimme’s DFT-D3(BJ).



THREEBODY (logical) {FALSE}
Inclusion of three-body term for Grimme’s DFT-D3 and DFT-D3(BJ).



D3H4S6 (real) {1.0e+00}
Parameter for Řezáč’s D3H4 correction (D3 dispersion).



D3H4RS6 (real) {1.215e+00}
Parameter for Řezáč’s D3H4 correction (D3 dispersion).



D3H4ALPHA (real) {3.0e+01}
Parameter for Řezáč’s D3H4 correction (D3 dispersion).



D3H4SHH (real) {3.0e−01}
Parameter for Řezáč’s D3H4 correction (H–H repulsion).



D3H4EHH (real) {1.431e+01}
Parameter for Řezáč’s D3H4 correction (H–H repulsion).



D3H4R0HH (real) {2.35e+00}
Parameter for Řezáč’s D3H4 correction (H–H repulsion).



D3H4CNN (real) {2.01e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).



D3H4CNO (real) {8.0e−01}
Parameter for Řezáč’s D3H4 correction (H4 correction).



D3H4CON (real) {2.58e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).
28



D3H4COO (real) {1.11e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).



D3H4CWAT (real) {1.32e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).



D3H4CSNHR3 (real) {2.33e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).



D3H4CSCOO (real) {1.22e+00}
Parameter for Řezáč’s D3H4 correction (H4 correction).



ULGS (real) {7.012e−01}
Parameter for DFT-ulg.
Note that default value is adjusted to PBE functional.



ULGB (real) {6.966e−01}
Parameter for DFT-ulg.
Note that default value is adjusted to PBE functional.



DDMCA (real) {1.018e+00}
Parameter for dDMC.



DDMCB0 (real) {1.857e+00}
Parameter for dDMC.



DDMCS (real) {4.6e+01}
Parameter for dDMC.



D3H5S6 (real) {1.0e+00}
Parameter for Řezáč’s D3H5 correction (D3 dispersion).



D3H5S8 (real) {4.9e−01}
Parameter for Řezáč’s D3H5 correction (D3 dispersion).



D3H5ALPHA (real) {2.961e+01}
Parameter for Řezáč’s D3H5 correction (D3 dispersion).



D3H5SHH (real) {3.0e−01}
Parameter for Řezáč’s D3H5 correction (H–H repulsion).



D3H5EHH (real) {1.431e+01}
Parameter for Řezáč’s D3H5 correction (H–H repulsion).
29

D3H5R0HH (real) {2.35e+01}
Parameter for Řezáč’s D3H5 correction (H–H repulsion).



D3H5SR (real) {7.14e−01}
Parameter for Řezáč’s D3H5 correction (H5 correction).



D3H5SW (real) {2.5e−01}
Parameter for Řezáč’s D3H5 correction (H5 correction).



D3H5KNH (real) {1.8e−01}
Parameter for Řezáč’s D3H5 correction (H5 correction).



D3H5KOH (real) {6.0e−02}
Parameter for Řezáč’s D3H5 correction (H5 correction).



D3H5KSH (real) {2.1e−01}
Parameter for Řezáč’s D3H5 correction (H5 correction).



D3XS8 (real) {3.209e+00}
Parameter for Kubillus’s D3X correction (D3 dispersion).



D3XA1 (real) {7.46e−01}
Parameter for Kubillus’s D3X correction (D3 dispersion).



D3XA2 (real) {4.191e+00}
Parameter for Kubillus’s D3X correction (D3 dispersion).



D3XC1 (real) {7.761e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XC2 (real) {5.0e−02}
Parameter for Kubillus’s D3X correction (X correction).



D3XC3 (real) {4.518e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XDNCL (real) {1.526e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XDNBR (real) {1.349e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XDNI (real) {1.521e+00}30

Parameter for Kubillus’s D3X correction (X correction).


D3XDOCL (real) {1.237e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XDOBR (real) {1.099e+00}
Parameter for Kubillus’s D3X correction (X correction).



D3XDOI (real) {1.313e+00}
Parameter for Kubillus’s D3X correction (X correction).



H4IDOMAIN (logical) {FALSE}
Enables the use of CLL algorithm for index of hydrogen atoms for acceleratingŘezáč’s D3H4 calculation (H4 correction) in DISPTYPE=6.



H4IDOMAINBLOCKXYZ/H4IDOMAINBLOCK (integer) {5}
Block size in each of x-, y-, and z-directions used in H4IDOMAIN=TRUE.



H4IDOMAINBLOCKX (integer) {5}
Block size in x-direction used in H4IDOMAIN =TRUE.



H4IDOMAINBLOCKY (integer) {5}
Block size in y-direction used in H4IDOMAIN =TRUE.



H4IDOMAINBLOCKZ (integer) {5}
Block size in z-direction used in H4IDOMAIN =TRUE.



H4IDOMAINNEIGHBOR (integer) {1}
Number of linked neighbors referred to as adjacent cells in H4IDOMAIN=TRUE.



H4VDOMAIN (logical) {FALSE}
Enables the use of CLL algorithm for evaluation of atom valence for acceleratingŘezáč’s D3H4 calculation (H4 correction) in DISPTYPE=6.



H4VDOMAINBLOCKXYZ/H4VDOMAINBLOCK (integer) {5}
Block size in each of x-, y-, and z-directions used in H4VDOMAIN=TRUE.



H4VDOMAINBLOCKX (integer) {5}
Block size in x-direction used in H4VDOMAIN =TRUE.



H4VDOMAINBLOCKY (integer) {5}
Block size in y-direction used in H4VDOMAIN =TRUE.
31

H4VDOMAINBLOCKZ (integer) {5}
Block size in z-direction used in H4VDOMAIN =TRUE.



H4VDOMAINNEIGHBOR (integer) {1}
Number of linked neighbors referred to as adjacent cells in H4VDOMAIN=TRUE.

1.1.8 PBC
The keyword PBC sets some parameters to control the calculation in periodic boundary
condition. Note that the system is treated as periodic if three translation vectors aredetected in Geometry section even if PBC=TRUE is not explicitly written (see section
1.4.1).
STRESS (logical) {FALSE}
Enables to calculate stress tensors and lattice derivatives.



COULOMBTYPE (integer) {1}
Method of Coulomb interaction calculation.
 1: Ewald summation
 2: Multipole expansion [41]See reference [14] for theory, implementation, and numerical assessmentsof the present option.
There are several limitations at present.
 The stress tensor calculation is not supported.
 Point charges should not be present (i.e. POINTCHARGE=FALSE in
SCC keyword).
 Atomresolvedenergies
cannotbeevaluated(i.e.
PRINTATOME=FALSE in MISC keyword).
 Shell
resolved
calculation
in
MPI/OpenMP
version(RUNTYPE=SHELL, RUNTYPE=SHELL0, or ORSCC=TRUE in
SCC keyword) is not available.
 Only cubic unit cell has been tested.



EWALDALPHA (real) {0.0e+00}
Fraction of Ewald summation between reciprocal and real space.
If value is larger than zero and smaller than one, manually specified value is used(in only single point calculation).
32

Otherwise, value is automatically determined (default and recommended).


MAXGAMMA (integer) {7}
Maximum number of unit cell repetition for Ewald summation.



FMMP (integer) {8}
Expansion order for multipole expansion in COULOMBTYPE=2.



FMMN (integer) {3}
Number of aggregated units at each level of telescoping procedure in
COULOMBTYPE=2.



FMMK (integer) {5}
Maximum number of level for multipole expansion in COULOMBTYPE=2.



EWALDSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for computation of reciprocal space summationin Ewald method in compensation for the increased computational cost.



EWALDMEMORYGB (real) {8.0e+00}
Maximum memory size for temporary arrays used in reciprocal space summationin Ewald method in gigabyte (GB) unit.
If the estimated memory size for temporary arrays is larger than the specified
value, the reduced memory usage approach is used in computations irrespectiveof EWALDSAVEMEMORY=FALSE.

FMMSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for computation of solid harmonics inmultipole-based algorithm (COULOMBTYPE=2) in compensation for theincreased computational cost.

1.1.9 MISC
The keyword MISC sets other miscellaneous options.


PRINTFORCE/FORCE (logical) {FALSE}
Enables to perform the force calculation after single point energy calculation.
If system includes external point charges, forces of point charges are additionallyprinted.
If STRESS=TRUE in PBC keyword, stress tensor, lattice derivatives, volume
and pressure are additionally printed.
33



PRINTRESOLVEDFORCE (logical) {FALSE}
Enables to evaluate force contributions (electronic, repulsive, and dispersion
term if included).
If STRESS=TRUE in PBC keyword, corresponding contributions of stresstensor are additionally evaluated.
Note that total forces are printed in the standard output while the respective
contributions are printed in the detailed data file.



PRINTMO (logical) {FALSE}
Enables to print eigenvectors.
Note that this option is available when the number of subsystem is equal to 1.



PRINTATOME (logical) {FALSE}
Enables to evaluate atomic energy contributions (electronic, repulsive, anddispersion term if included).



PRINTHS (logical) {FALSE}
Enables to print the lower triangle of non-perturbed Hamiltonian and overlap
matrices in the standard output.
This option is available when the number of subsystem is equal to 1.
Note that the program terminates immediately without proceeding further DFTB
computational steps.



PRINTDENSITY (logical) {FALSE}
Enables to print the lower triangle of density matrix in the standard output.



PRINTLOCALREGION (logical) {FALSE}
Enables to print number of subsystems, atoms, and basis functions in localization
regions in the standard output.
Note that the program terminates immediately without proceeding further DFTB
computational steps.



PRINTLAPTIME (logical) {FALSE}
Enables to print wall clock time for DFTB computational steps.
Total wall clock time is additionally printed at the end of program run.



WRITEDAT (method) {APPEND}
Enables to control dumping of the detailed data file.
 NONE never generates the data file.
 APPEND opens the data file in the beginning of program run and writes the34

information of each DFTB calculation. The size of detailed data file may
become huge if a large number of geometry steps is specified in geometry
optimization or MD.
 REPLACE opens/closes the data file at the beginning/end of DFTB
calculation. Namely, in geometry optimization or MD, the data file contains
the information of the last geometry step.


WRITECHRGFILE (logical) {TRUE}
Enables to control dumping of the binary file “chrgfile” containing Mulliken
charge information which is used in READCHARGE=BINARY in SCC
keyword.



USESCRATCH (logical) {TRUE}
Enables to use disk I/O for storing some intermediate data.
If USESCRATCH=FALSE, all intermediate data are stored in memory.

35

1.2 Title section
Comment line

TITLE

This is typically used for the comment of compound name, symmetry, electronic state,other relevant information, and so on. Each line should be less than 128 characters. Allcomments will be written in the standard output file as it is.

36

1.3 Parameter section
1.3.1 General overview
2
O˽2
oo.spl˽oh.spl
H˽1
ho.spl˽hh.spl

Number of unique elements
Atom symbol and its highest angular momentum for valenceshellsSlater-Koster parameters: atom1–atom1 and atom1–atom2
Atom symbol and its highest angular momentum for valenceshellsSlater-Koster parameters: atom2–atom1 and atom2–atom2

This section starts from the line giving the number of elements in the given structure. The
present example is 2 since water molecule contains oxygen and hydrogen. Next, atomic
symbol and its maximum angular momentum are specified for the first element (atom1).
The integer 1, 2, 3 corresponds to the maximum angular momentum s, p, d, respectively.
Then, the paths to Slater-Koster parameter files including the first element are specified
in a single line. The order has to be consistent in the order of elements specified in thissection, namely atom1 with atom1, atom1 with atom2 (the second element), atom1 with
atom3 (the third element), … These two lines are repeated until all elements in the given
structure are specified. Note that the present example assumes 4 Slater-Koster parameter
files (oo.spl, oh.spl, ho.spl, hh.spl) are located in the same directory as the one containing
the input file.
1.3.2 Specification of additional constants
Several options request atom dependent constants which are not hard-coded or written in
Slater-Koster files. Such constants are specified in two ways:

Direct specification of employed values
This is necessary for Hubbard derivatives in third order corrections(THIRDFULL=TRUE and THIRDDIAG=TRUE in SCC keyword).
The Hubbard derivative should be specified after the maximum angularmomentum.
In the case of orbitally resolved calculation (ORSCC=TRUE in SCC keyword)with Multi version, values of s-, p-, and d-shells need to be provided separately.

37

O˽2˽−0.1575
Cu˽3˽−0.20˽−0.0575˽−0.0575



Atom resolved case (usual)
Orbitally resolved case(order is d, p, s)

Preparation of additional ASCII files
This is necessary for the following three options: Spin constants for spin polarized calculation (SPIN=TRUE in SCC
keyword) Atomic polarizabilities, van der Waals radii, and effective charge for SlaterKirkwood dispersion (DISPTYPE=2 in DISP keyword) Lennard-Jones parameters for Lennard-Jones dispersion (DISPTYPE=3
and READUFF=TRUE in DISP keyword)
The path of additional file should be specified after the maximum angularmomentum (or Hubbard derivatives if they exist). In the case of specifying both
files, the path of dispersion constants file should come first, followed by the pathof spin constants file:
Only spin constants file
Only dispersion constants
file
Both spin and dispersion
constants filesO˽2˽o.wll
O˽2˽o.dispsk
O˽2˽o.dispsk˽o.wll

Here is an example of the spin constants file for carbon:−0.031
−0.025
−0.025
−0.023
The number of lines of spin constants file depends on the maximum angularmomentum (1, 4, 9 for s, p, d). Each line contains the atomic spin constant inatomic unit in the following order: s: ss p: ss, sp, ps, pp38 d: ss, sp, sd, ps, pp, pd, ds, dp, dd
For reference, PBE values for H, C, N, O, S, Fe, and Ni listed in DFTB+ manual[42] are stored into the file with suffix “.wll” in the slko directory.
Here is an example of the Slater-Kirkwood dispersion constants file for carbon:
1.382˽1.382˽1.382˽1.064˽3.8˽3.8˽3.8˽3.8˽2.5
There should exist nine columns in a single line. The first to fourth columns
specify atomic polarizabilities in Å3 unit, the fifth to eighth columns specify van
der Waals radii in Å unit, and the last column specifies effective charge. The
values for atomic polarizabilities and van der Waals radii are given for zero, one,
two, and three hydrogen neighbors. For reference, H, C, N, O, and S dispersion
constants listed in DFTB+ manual [42] are stored into the file with suffix“.dispsk” in the slko directory.
To manually specify the Lennard-Jones parameters, values of non-bonded
distance (in Å unit) and energy (in kcal/mol unit) are given in a single line.
The order of specification
is distance then energy

3.851˽0.105

For reference, H, C, N, O, S, Fe, and Ni Lennard-Jones parameters listed inreference [32] are stored into the file with suffix “.displj” in the slko directory.

39

1.4 Geometry section
1.4.1 General overview
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000

Number of atoms, charge,
spin multiplicity
Cartesian coordinates ofatom label 1
Cartesian coordinates ofatom label 2
Cartesian coordinates ofatom label 3This section begins with the line giving three integers; number of atoms, charge, and spin
multiplicity for the given structure. Then the position of atom is specified with Cartesian
coordinates in Å unit with floating point number (DO NOT USE INTEGERS!).
If the given system is periodic, the unit cell can be specified by adding three translation
vectors with symbol TV after the last Cartesian coordinates (TV should be an upper-case).
The following example is for diamond crystal structure (primitive unit cell).
2˽0˽1
C0˽−0.00000000˽−0.00000000˽−0.00000000
C0˽−0.89169750˽−0.89169750˽−0.89169750
TV˽−0.00000000˽−1.78339500˽−1.78339500
TV˽−1.78339500˽−0.00000000˽−1.78339500
TV˽−1.78339500˽−1.78339500˽−0.00000000
If the external static electric field is included, the electric field vector can be specified byadding one line with symbol EF after the last Cartesian coordinates (EF should be an
upper-case). The example of water molecule with electric field components of x = 0.001,
y = −0.001, and z = 0.000 is given below.
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
40H0˽−0.08185946˽−1.01910555˽−0.00000000
EF˽−0.00100000˽−0.00100000˽−0.00000000
If the point charges are included, the coordinates and charge value can be specified byadding lines with symbol PC after the last Cartesian coordinates (PC should be an uppercase). The example of water molecule with point charge of 1.0 positioned at x = 5.0 Å, y
= 5.0 Å, and z = 5.0 Å is given below.
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000
PC˽−5.00000000˽−5.00000000˽−5.00000000˽−1.00000000
1.4.2 Specification of additional information
Additional columns after Cartesian coordinates can be read when you would like tospecify atoms to be moved or atomic masses to be changed to investigate the isotope
effects of vibrational properties.


Freeze the coordinate in geometry optimization or MD
If a single asterisk (*) is put, the position of atom will not be moved.

Constraint geometry optimization (CONSTRAINT=TRUE in OPT keyword)
Constraint vector of x, y, z components is either 0.0 (fixed) or 1.0 (moved). For
example, the vector “0.0˽0.0˽1.0” will update only z component. The vector“0.0˽0.0˽0.0” is identical to put a single asterisk.

Specification of different atomic masses in harmonic vibrational frequency
analysis (READISOTOPE=TRUE in FREQ keyword)
The form ISO=value changes the default atomic mass to value. No space
between ISO= and value. This specification should come later than the single
asterisk or constraint vector when both options are activated.

Here is an example of the mixed specifications for water molecule.
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
41H0˽−1.01910555˽−0.08185946˽−0.00000000˽0.0˽1.0˽0.0
H0˽−0.08185946˽−1.01910555˽−0.00000000˽*˽ISO=2.014102
When the constraint geometry optimization is performed, the hydrogen with atom label 2
can be moved along only the y-axis and the coordinate of hydrogen with atom label 3 isnever updated. If FREQ=(READISOTOPE=TRUE) is additionally specified, theatomic mass of hydrogen with atom label 3 is changed to 2.014102, namely the harmonic
vibrational frequencies and corresponding thermochemical properties of HOD molecule
will be obtained.
1.4.3 Manual construction of subsystems in DC-DFTB
This section supplements how to specify the subsystems using SUBTYPE=SPECIFY
and SUBTYPE=SEMIAUTO options in DC keyword. After the blank line of Geometry
section, the indices of subsystems for each atom are listed.
Here is an example of specification of 5 subsystems for the given structure containing 22
atoms:
1˽1˽1˽1˽1
2˽2˽2˽2
3˽3˽3˽3
4˽4˽4˽4
5˽5˽5˽5˽5

Index of subsystem for atom label 1, 2, 3, 4, 5
Index of subsystem for atom label 6, 7, 8, 9
Index of subsystem for atom label 10, 11, 12, 13
Index of subsystem for atom label 14, 15, 16, 17
Index of subsystem for atom label 18, 19, 20, 21, 22

where atoms with label 1–5 belong to subsystem 1, atoms with label 6−9 belong tosubsystem 2, atoms with label 10−13 belong to subsystem 3, and so on. It is notnecessarily to use the same number of lines as the number of subsystems. The new line
will be read until the total number of indices is equal to the number of atoms in the given
structure.
In the case of SUBTYPE=SEMIAUTO, the integer “0” is used as an index for atoms tobe assigned to subsystems automatically by using the grid box. The following is an
example of system containing 22 atoms (using analogous format from the previous one
for comparison):1˽1˽1˽1˽1

Index of subsystem for atom label 1, 2, 3, 4, 542

0˽0˽0˽0
0˽0˽0˽00˽0˽0˽02˽2˽2˽2˽2

Index of subsystem for atom label 6, 7, 8, 9
Index of subsystem for atom label 10, 11, 12, 13
Index of subsystem for atom label 14, 15, 16, 17
Index of subsystem for atom label 18, 19, 20, 21, 22

where atoms with label 1–5 belong to subsystem 1, atoms with label 18–22 belong tosubsystem 2, and atoms with label 6–17 are assigned to subsystems 3–X automatically(X depends on the given structure and size of grid box). If all indices are zero and positive
integer, the option is identical to SUBTYPE=AUTO and SUBTYPE=SPECIFY,
respectively.

43

2 Running DCDFTBMD
Once the input file “dftb.inp” is created, other requisites such as Slater-Koster parameters
and restart files are needed to be placed in the appropriate directory. For example, in thecase of SCC-DFTB single point energy calculation of water molecule shown in section 1,
4 Slater-Koster parameter files (oo.spl, oh.spl, ho.spl, and hh.spl) must be placed in thesame directory as the input file. The executable is used without any arguments.
Serial version
OpenMP version
MPI/OpenMP version<installation path>/dftb_serial.00.x<installation path>/dftb_openmp.00.x
mpiexec –np (number of processes) <installation
path>/dftb_mpiomp.00.x

The header of input file may be specified arbitrarily by setting it as the first argument.
The example of using “file.inp” is as follows.
Serial version<installation path>/dftb_serial.00.x file

OpenMP version
MPI/OpenMP version<installation path>/dftb_openmp.00.x filempiexec –np (number of processes) <installation
path>/dftb_mpiomp.00.x file

For OpenMP parallelization, number of threads must be set before starting the calculation.
bash
csh, tcsh

export OMP_NUM_THREADS=(number of threads)
setenv OMP_NUM_THREADS (number of threads)

If a large-scale calculation stops, one may need to increase the stack size.
bash
csh, tcsh

ulimit –s unlimited; export OMP_STACKSIZE=(stack size value)
unlimit; setenv OMP_STACKSIZE (stack size value)

Note that the maximum number of basis sets in serial and OpenMP versions is limited to
16384.

44

3 Output of DCDFTBMD
If the calculation terminates normally, the standard output file and the detailed data file
will be appeared. By default (using “dftb.inp” as the input file), the corresponding
filenames are “dftb.out” and “dftb.dat”, respectively. When the header of input file is
specified arbitrarily (e.g. “file.inp” as described in section 2), the header of output anddata files will be changed appropriately (e.g. “file.out” and “file.dat”). Here, the output
files of serial and OpenMP versions for SCC-DFTB single point energy calculation ofwater molecule shown in section 1 are briefly outlined. Those of MPI/OpenMP version
are similar to serial and OpenMP versions although the format of detailed data file isslightly different for the information of subsystem and Mulliken population.

3.1 Standard output file
The standard output file begins with the summary of calculation type.
*********************************************************************************************************************
DDDDDD

CCCCCCDDDDDD

FFFFFFFF TTTTTTTTT BBBBBBBDDDDDDDD

CCCCCCCC DDDDDDDD FFFFFFFF TTTTTTTTT BBB

MMMM DDDDDBBB MMMMMM DDDDDDD

DDD

DDD CCCDDD

DDD FFTTTBBBBBB MM M M MM DDD

DDDDDD CCC

DDD

DDD

FFFFFFTTTBBBBBBBB MMDDDDDD CCC

DDD

DDD

FFFFFFTTTBBBMM DDDDDD

BBB MM

MM DDDDDD

BB MM

MM DDDDDDD

DDDDDDDD

CCCCCCCC DDDDDDDD

FF

TTT

BBB

DDDDDD

CCCCCCFFTTTBBBBBBBBDDDDDDMM

M

DDD

MM DDDDD

VERSION 1.0 (Nov 2018)*********************************************************************************************************************
Execution of DCDFTBMD begun %a %b %d %H:%M:%S %Y-------TITLE-------2 kinds of elements will be used
Pair of element 1 and element 1 : {path}/oo.spl
Pair of element 1 and element 2 : {path}/oh.spl
Pair of element 2 and element 1 : {path}/ho.spl
Pair of element 2 and element 2 : {path}/hh.spl
SCC= True45

------------------------Max cycle= 200

Energy convergence= 1.0E-09

Density convergence= 1.0E-06

Charge mixing method= Broyden

Broyden iteration= 70

Mixing parameter= 2.0E-01

Eigensolver= DSYGVD

Damping gamma X-H= FalseOrbital resolved SCC= FalseThird order diagonal= FalseThird order full

= FalseSpin contribution= FalseMayer bond order= FalsePolarizability= FalseExternal field= FalsePoint charge= FalsePseudodiagonalization = False
Cubic splined gamma = False
Direct SCC= FalseRead charge= NONEDC= False------------------------Optimize = False------------------------MD= False------------------------Frequency = False------------------------Dispersion = False------------------------PBC= False------------------------MISC= False------------------------Total number of basis set shells = 4

46

Number of basis functions=6

Number of electrons=8

Charge of system=0

Spin multiplicity=1

Number of occupied orbitals=4

Total number of atoms=3*** Start single point energy calculation ***Then, the SCC-DFTB single point calculation starts in the present example.
*********************************** Start SCC-DFTB calculation ***********************************
Iter. Total energy

Energy diff.

Density diff.

Fermi level------------------------------------------------------------------------------------------------------1-0.0726464067-0.0726464067-2.0000000000-0.0025000000

2-4.0712418245-4.1438882313-0.0200187809-0.0025000000

3-4.0741988428-0.0029570183-0.0590590622-0.3564583333

4-4.0779379555-0.0037391127-0.0004103232-0.3525000000

5-4.0779381355-0.0000001800-0.0000024570-0.3525000000

6-4.0779381355-0.0000000000-0.0000000001-0.3562500000

Final SCC-DFTB Energy = -4.0779381355 Eh after 6 iterations
Execution of DCDFTBMD terminated normally %a %b %d %H:%M:%S %Y

Here, the SCC converged total energy is −4.0779381355 Hartree.

47

3.2 Detailed data file
The detailed data file contains the relevant properties of the given structure. The
constructed subsystems and localization regions are initially printed. Since thiscalculation is DC=FALSE, all 3 atoms belong to subsystem 1.
SUBSYSTEM : 1------------------------1

2

3

LOCALIZATION REGION : 1----------------------------------------1

If you use MPI/OpenMP version, this part looks concise excluding the case of
PRINTLOCALREGION=TRUE in MISC keyword.
SUBSYSTEM--------------------1

1

1

Then, the eigenvalues of the molecular orbitals (first column) and the electron occupation
number of that level (second column) are shown.
EIGENVALUES [Eh] AND OCCUPANCIES------------------------------------------------------------0.8489865909

2.0000000000-0.4143390336

2.0000000000-0.3137537414

2.0000000000-0.2591755518

2.0000000000-0.3992760730

0.0000000000-0.5583952222

0.0000000000

If PRINTMO=TRUE in MISC keyword, the eigenvectors are additionally printed. The
occupied and virtual orbitals are distinguished with OCC. and VIRT., respectively.
EIGENVALUES [Eh], OCCUPANCIES, AND MO COEFFICIENTS

48

-----------------------------------------------------------------------------------------1

OCC.
1

O

-0.8489865909

2.0000000000

S

0 -0.8522740053

P-10.0017546133

P00.0000000000

P

10.0017546133

2

H

S

0 -0.1530727458

3

H

S

0 -0.1530727458

2

OCC.
1

O

-0.4143390336

2.0000000000

S

0 -0.0000000000

P-10.4626963766

P

00.0000000000

P

1 -0.4626963766

2

H

S

0 -0.3908041098

3

H

S

0

0.3908041098…

Then, the total Mulliken charge on the atom (second column) is shown together with thecorresponding net charges (third column).
MULLIKEN POPULATIONS AND NET ATOMIC CHARGES---------------------------------------------------------------------------------1

6.5926170634-0.5926170634

20.7036914683-0.2963085317

30.7036914683-0.2963085317

If ORBITALPOP=TRUE in SCC or NCC keyword, the Mulliken charges on each atomic
orbital are additionally printed.
MULLIKEN POPULATIONS AND NET ATOMIC CHARGES---------------------------------------------------------------------------------1

O

6.5926170634

S

0 1.7342141199

P-1 1.4292014718

P

0 2.0000000000

P

1 1.4292014718-0.5926170634

49

2

H

0.7036914683
S

3-0.2963085317

0 0.7036914683

H

0.7036914683
S-0.2963085317

0 0.7036914683

If you use MPI/OpenMP version, the shell resolved Mulliken population is printed.
MULLIKEN POPULATIONS--------------------------------------1.73421412

4.85840294

0.70369147

0.70369147

Then, the x, y, z dipole moment components and the norm using the Mulliken charges areshown in both atomic unit and Debye unit. This is not printed if periodic boundary
condition is used.
DIPOLE MOMENT [a.u.]----------------------------------xy0.4545263289

0.4545263289

NORM :

z0.0000000000

0.6427972988

DIPOLE MOMENT [Debye]--------------------------------------xyz

1.1552906835

1.1552906835

0.0000000000

NORM :

1.6338277531

If MAYER=TRUE in SCC or NCC keyword, atomic valences (diagonal elements) andMayer bond orders (off-diagonal elements) are shown with matrix representation.
ATOMIC VALENCIES AND MAYER BOND ORDERS-------------------------------------------------------------------------1

2

3

1

1.796085

0.898042

0.898042

2

0.898042

0.912201

0.014159

3

0.898042

0.014159

0.91220150

If CM3=TRUE in SCC or NCC keyword, the CM3 net atomic charges and dipole
moments are additionally printed.
CM3 NET ATOMIC CHARGES-------------------------------------------1-0.5833492665

2-0.2916746332

3-0.2916746332

CM3 DIPOLE MOMENT [a.u.]-----------------------------------------xyz

0.4474181000

0.4474181000

0.0000000000

NORM :

0.6327447451

CM3 DIPOLE MOMENT [Debye]---------------------------------------------xy1.1372233679
NORM :

z1.1372233679

0.0000000000

1.6082767104

Then, the different energy contributions to the total energy are shown in Hartree unit. If
you use MPI/OpenMP version, PRINTATOME=TRUE in MISC keyword is necessary
to obtain the corresponding information.
ENERGY CONTRIBUTION [Eh]--------------------------------------------ENERGY H0= -4.1689461692

ENERGY SCC= -0.0183616273

TOTAL ELECTRONIC ENERGY = -4.1505845420
REPULSIVE ENERGY= -0.0726464065

TOTAL ENERGY= -4.0779381355

If PRINTATOME=TRUE in MISC keyword, the atomic energy contributions to the total,electronic, and repulsive energies are additionally printed in Hartree unit.
ATOM RESOLVED ENERGY CONTRIBUTION [Eh]------------------------------------------------------------------------

51

TOTAL

ELECTRONICREPULSIVE

1-3.3763986686-3.4127218718-0.0363232032

2-0.3507697335-0.3689313351-0.0181616016

3-0.3507697335-0.3689313351-0.0181616016

52

3.3 Other output files
Additional output files will be generated depending on the type of calculations:“chrgfile” (if SCC=TRUE)
contains the information of Mulliken population for restarting purpose in binary
format. The file may be used in READCHARGE=BINARY in SCC keyword.

“restart” (if MD=TRUE)
contains the information of last MD step for restarting purpose in binary format.
The file may be used in RESTART=TRUE in MD keyword.

“restart_chk” (if MD=TRUE)
contains the information of last checkpoint step for restarting purpose in binary
format. The file may be used in RESTART=TRUE in MD keyword.

“traject” (if MD=TRUE)
contains the Cartesian coordinates during MD simulation in xyz format.

“velocity” (if MD=TRUE)
contains the velocity components during MD simulation in xyz format in Å/fs
unit.



“mulliken” (if MD=TRUE)
contains the Mulliken net atomic charges during MD simulation. For serial and
OpenMP versions, sequence of the following information is printed:
 Line 1: Number of atoms
 Line 2: Comments (simulation time and MD step) Line 3–(Number of atoms + 2): Atom index (1st column), atom label (2ndcolumn), and Mulliken net atomic charges (3rd column)
For MPI/OpenMP version, sequence of the following information is printed:
 Line 1: Nunber of shells (1st column) and number of atoms (2nd column)
 Line 2: Comments (simulation time and MD step) Line 3–(Number of shells + 2): Atom index (1st column), atom label (2ndcolumn), angular momentum label (3rd column), and Mulliken net charges
for the shell (4th column)

“dftb_vib.out” (if FREQ=TRUE)
53

contains the information of harmonic vibrational frequency analysis to visualize
the results using GaussView.
“restart_hess” (if PRINTHESS=TRUE in FREQ keyword)
contains the information of Hessian matrix elements for restarting purpose. The
file may be used in RESTARTHESS=BINARY in FREQ keyword.

“hess.dat” (if INTERRUPTCP=TRUE in FREQ keyword)
contains the information of partial Hessian matrix elements for restarting
purpose. The file may be used in RESTARTHESS=ASCII in FREQ keyword.

You may find other files (zacaz, zbcbz, zzszz) in serial and OpenMP versions if
USESCRATCH=TRUE in MISC keyword but these are intermediate binary files andnot important to examine the results.

54

References[1]

M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S.

[2]
[3]
[4]
[5]

Suhai, G. Seifert, Phys. Rev. B 58, 7260 (1998).
M. Elstner, J. Phys. Chem. A 111, 5614 (2007).
Y. Yang, H. Yu, D. York, Q. Cui, M. Elstner, J. Phys. Chem. A 111, 10861 (2007).
M. Gaus, Q. Cui, M. Elstner, J. Chem. Theory Comput. 7, 931 (2011).
C. Köhler, G. Seifert, U. Gerstmann, M. Elstner, H. Overhof, T. Frauenheim, Phys.

Chem. Chem. Phys. 3, 5109 (2001).
[6] C. Köhler, G. Seifert, T. Frauenheim, Chem. Phys. 309, 23 (2005).
[7] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[8] V. Eyert, J. Comput. Phys. 124, 271 (1996).
[9] P. Pulay, J. Comput. Chem. 3, 556 (1982).
[10] I. Mayer, Chem. Phys. Lett. 97, 270 (1983).
[11] J. A. Kalinowski, B. Lesyng, J. D. Thompson, C. J. Cramer, D. G. Truhlar, J. Phys.
Chem. A 108, 2545 (2004).
[12] S. Kaminski, M. Gaus, M. Elstner, J. Phys. Chem. A 116, 11927 (2012).
[13] J. J. P. Stewart, P. Császár, P. Pulay, J. Comput. Chem. 3, 227 (1982).
[14] Y. Nishimura, H. Nakai, J. Comput. Chem. 39, 105 (2018).
[15] Y. Yang, H. Yu, D. York, M. Elstner, Q. Cui, J. Chem. Theory Comput. 4, 2067
(2008).
[16] M. Gaus, H. Jin, D. Demapan, A. S. Christensen, P. Goyal, M. Elstner, Q. Cui, J.
Chem. Theory Comput. 11, 4205 (2015).
[17] W. Yang, T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).
[18] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational
Chemistry and Physics, R. Zalesny, M. G. Papadopoulos, P. G. Mezey, J.
Leszczynski, Eds.; Springer: Dordrecht, pp. 97–127 (2011).
[19] B. Cordero, V. Gómez, A. E. Platero-Prats, M. Revés, J. Echeverría, E. Cremades,
F. Barragán, S. Alvarez, Dalton Trans. 21, 2832 (2008).
[20] H. Nishizawa, Y. Nishimura, M. Kobayashi, S. Irle, H. Nakai, J. Comput. Chem.
37, 1983 (2016).
[21] L. V. Woodcock, Chem. Phys. Lett. 10, 257 (1971).
[22] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, M. L. Klein, Mol. Phys. 87, 1117
(1996).
[23] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak,
J. Chem. Phys. 81, 3684 (1984).
55

[24] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).
[25] H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
[26] G. Marsaglia, A. Zaman, W. W. Tsang, Stat. Probab. Lett. 8, 35 (1990).
[27] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B 51,
12947 (1995).
[28] G. Seifert, D. Porezag, T. Frauenheim, Int. J. Quantum Chem. 58, 185 (1996).
[29] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
[30] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 114,
5149 (2001).
[31] L. Zhechkov, T. Heine, S. Patchkovskii, G. Seifert, H. A. Duarte, J. Chem. Theory
Comput. 1, 841 (2005).
[32] A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, W. M. Skiff, J. Am.
Chem. Soc. 114, 10024 (1992).
[33] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 132, 154104 (2010).
[34] S. Grimme, S. Ehrlich, L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
[35] J. Řezáč, P. Hobza, J. Chem. Theory Comput. 8, 141 (2012).
[36] V. M. Miriyala, J. Řezáč, J. Comput. Chem. 38, 688 (2017).
[37] H. Kim, J.-M. Choi, W. A. Goddard III, J. Phys. Chem. Lett. 3, 360 (2012).
[38] R. Petraglia, S. N. Steinmann, C. Corminboeuf, Int. J. Quantum Chem. 115, 1265
(2015).
[39] J. Řezáč, J. Chem. Theory Comput. 13, 4804 (2017).
[40] M. Kubillus, T. Kubař, M. Gaus, J. Řezáč, M. Elstner, J. Chem. Theory Comput.
11, 332 (2015).
[41] C. G. Lambert, T. A. Darden, J. A. Board Jr., J. Comput. Phys. 126, 274 (1996).
[42] http://www.dftbplus.org/fileadmin/DFTBPLUS/public/dftbplus/latest/manual.pdf

56