DCDFTBMDDivide-and-Conquer Density Functional Tight-Binding Molecular Dynamics

User Manual
Version 2.0
December, 2019

Nakai Group
Department of Chemistry and Biochemistry,
School of Advanced Science and Engineering, Waseda University

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 DC-DFTB with respect tosystem size is close to linear. A high parallel efficiency has been achieved in theimplementation with the hybrid of MPI and OpenMP. The program has
functionalities not only for MD but also for geometry optimization, vibrationalfrequency calculation, and property evaluation.

DCDFTBMD website
General information about DCDFTBMD, including development background,program download, documentation, and publication can be found on thefollowing website:
http://www.chem.waseda.ac.jp/dcdftbmd/en/Improvements from old version
New features added to version 2.0






Metadynamics simulation
Geometry optimization with Quick-min or FIRE method
LSQ method for accelerating SCC convergence
Interruption and restarting numerical Hessian calculation
Numerical gradient calculation






Treating extended Slater-Koster parameters and f-electrons
DFTB+U approach for SCC calculation
Onsite correction for SCC calculation
D3H5 correction for F element

Below is available for only MPI/OpenMP version:
 MPI-based Broyden charge mixing for reducing memory usage
 Better efficiency of orbitally resolved SCC calculation
 Spline interpolation algorithm for orbitally resolved SCC calculation
Default settings changed from version 1.0Threshold for scaling velocities in constant temperature MD simulation
Old: 100 K → New: 10000 K
Step size for numerical differentiation of DFTB parameterized integrals to
obtain first and second derivatives
Old: Fixed (0.01 atomic unit) → New: Adjustable parameter
Step size for numerical differentiation of D3 gradient to obtain second
derivatives
Old: Fixed (0.0001 atomic unit) → New: Adjustable parameter

Reference Hubbard value to identify hydrogen in damped γ-function
Old: Fixed (0.4195 atomic unit) → New: Adjustable parameter
Below is changed for only serial and OpenMP versions: No use of disk I/O for storing some intermediate dataBugs fixed from version 1.0
 DFTB3 atom resolved energy Zero charge guess for open-shell system
Below is fixed for only MPI/OpenMP version:

Hubbard derivative parameter check in orbitally resolved SCC calculation

Contents
1 Input for DCDFTBMD

1

1.1 Keyword section2

1.1.1 General overview
1.1.2 SCC

2
3Basic settings
Charge mixing
Initial charge guess
SolverThird-order correction

3
3
4
4
5

Modified γ-function
Open-shell and spin-polarization correction
External field effect
DFTB+U correction

5
6
6
7

Onsite correction
Electronic property evaluation
Additional options for only MPI/OpenMP version
1.1.3 NCC
SolverOpen-shell
Electronic property evaluation
1.1.4 DC
Basic settings

7
7
8
10
10
10
10
11
11

Advanced options
Additional options for only MPI/OpenMP version
1.1.5 DISP/DISPERSION
Basic settings
Advanced options for Grimme’s DFT-D2
Advanced options for Slater-Kirkwood
Advanced options for Lennard-Jones
Advanced options for Grimme’s DFT-D3
Advanced options for Řezáč’s D3H4
Advanced options for DFT-ulg

11
12
14
14
15
15
15
15
16
17

i

Advanced options for dDMC
Advanced options for Řezáč’s D3H5
Advanced options for Kubillus’s D3X
Advanced options for frequency calculation
Additional options for only MPI/OpenMP version
1.1.6 PBC
Basic settings
Advanced options for Ewald summation

17
17
18
19
19
21
21
21

Advanced options for multipole expansion
Additional options for only MPI/OpenMP version
1.1.7 OPT/OPTIMIZE
Basic settings
Constraint optimization

22
22
23
23
24

Lattice optimization
Acceleration of SCC convergence
Advanced options for BFGS
Advanced options for steepest descent
Advanced options for conjugate gradient
Advanced options for Quick-min
Advanced options for FIRE
1.1.8 FREQ/FREQUENCY
Basic settings
Advanced options
Additional options for only MPI/OpenMP version
1.1.9 MD
Basic settings
Constant energy ensemble
Constant temperature ensemble
Advanced options for thermostat
Constant pressure ensemble
Advanced options for barostat
Random number generation
Soft boundary potential
Constraint MD with RATTLE
Acceleration of SCC convergence
Annealing simulation
ii

24
25
25
25
25
26
26
27
27
28
29
30
30
30
31
32
32
33
33
34
36
36
37

Metadynamics
Advanced options for metadynamics
1.1.10 MISC
Additional analysis and property evaluation
Global options
Additional options for only serial and OpenMP versions

37
38
40
40
41
42

1.2 Title section

43

1.3 Parameter section

44

1.3.1 General overview
1.3.2 Specification of additional constants
1.3.3 Format of ASCII files for additional constants

44
45
47

Spin constants fileU–J parameters file
Onsite parameters file
Slater-Kirkwood dispersion constants file

47
47
48
48

Lennard-Jones dispersion parameters file48

1.4 Geometry section

49

1.4.1 General overview
1.4.2 Specification of additional information
1.4.3 Manual construction of subsystems in DC-DFTB

49
51
52

1.5 Auxiliary input files

53

1.5.1 Binary format files

53

1.5.2 ASCII format files
Reading initial charge in SCC calculation
Reading Hessian matrix in frequency calculation
Reading initial velocity in MD simulation
Manual construction of bond constraints in RATTLE
Setting temperature profile in annealing simulation
Setting CVs in metadynamics simulation
Reading bias potential in metadynamics simulation
Setting wall potential in metadynamics simulation

54
54
55
56
56
57
58
69
70

iii

2 Running DCDFTBMD

71

3 Output of DCDFTBMD

72

3.1 Standard output file72

3.2 Detailed data file75

3.3 Other output files

80

3.3.1 Output files in DFTB calculation
3.3.2 Output files in frequency calculation
3.3.3 Output files in MD simulation

80
81
82

3.3.4 Output files in metadynamics

83

References

84

iv

Contents

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:Keyword section: Information about the calculation (section 1.1)

Title section: Brief description of the calculation (section 1.2)

Parameter section: Specification of Slater-Koster parameters (section 1.3)



Geometry section: Specification of molecular coordinates (section 1.4)

Each section is terminated with blank line (starting new line is indicated as  inthis document). For example, the input file looks like as follows in the case of SCCDFTB single point energy calculation of water molecule:
Keyword section

SCC=TRUE˽DC=FALSE

Title section

TITLE
Parameter section2
O˽2
oo.spl˽oh.spl
H˽1
ho.spl˽hh.spl


Geometry 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


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 ˽ in this document) should be written with spaces (DO NOT USE
TABS!).

1

Contents

1.1 Keyword section
1.1.1 General overview
SCC and DC keywords are

SCC=TRUE˽DC=FALSE

turned on and off

The present example specifies two keywords in a single line. The same
specification is possible by giving two separated lines as follows:
SCC=TRUESCC keyword is turned on

DC=FALSEDC keyword is turned off

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 be turned 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 inparenthesis stands for the data type of option. The options can be specified by
the following forms:
keyword=(option1=value1˽option2=value2˽…)

Setting option1 and option2
in a single line

keyword=(option1=value1)

Setting option1 and option 2

keyword=(option2=value2˽…)

in two separated lines

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

2

Contents

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

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



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



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

Charge mixing

MIXER (integer) {1}
Charge mixing method to accelerate SCC convergence.
Note that modified Broyden mixing is used to solve coupled perturbed
DFTB equation in analytical Hessian calculation.
 1: Modified Broyden [2] 2: Simple mixing 3: Modified Anderson [3] 4: DIIS [4]

ALMIX (real) {2.0e−01}
Charge mixing parameter used in charge mixing method.
The value should be larger than zero and less than one.

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



NGEN (integer) {4}
Number of iteration to start Anderson and DIIS charge mixing.
Simple mixing is used until reaching this value.

3

Contents

Initial charge guess


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
this option 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” (seesection 1.5.2 for the format). If this option 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.



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

Solver

SOLVER (integer) {1}
Eigensolver for diagonalization of Hamiltonian matrix.
 1: DSYGVD
 2: DSYGVPSEUDODIAG/PDIAG (logical) {FALSE}
Enables to perform psuedodiagonalization using Jacobi rotations [5].
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.
4

Contents

Third-order correction

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



THIRDFULL (logical) {FALSE}
Enables the use of full third order corrections [7].
Hubbard derivatives need to be specified in Parameter section (seesection 1.3.2).
The so-called DFTB3 model [7] requires damped γ-function for X–H pairs(DAMPXH=TRUE).

Modified γ-function

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



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



DAMPXHREFU (real) {4.195e−01}
Reference Hubbard value to identify hydrogen used in DAMPXH=TRUE.

UDERIVGAUSS (logical) {FALSE}
Enables to add charge dependent Gaussian function to Hubbard
derivative in on-site third-order correction (THIRDDIAG=TRUE) [9].



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.

5(positionofcenter)usedinContents



GAMMASPLINE (logical) {FALSE}
Enables to use cubic spline interpolation for exponentially decaying
short-range part of γ-function.
See reference [10] for details.
For MPI/OpenMP version, this is incompatible with one of shell resolved
calculation (RUNTYPE=SHELL0).

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

Open-shell and spin-polarization correctionSPIN (logical) {FALSE}
Enables to perform collinear spin polarized calculation [11, 12].
Spin constants have to be stored in different file for each atom (seesections 1.3.2 and 1.3.3).



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



NELECTHRESH (real) {1.0e−06}
Threshold for conservation
SPINRELAX=TRUE.

of

total

numberofelectronsin

UNPAIREDELEC (real) {0.0e+00}
Number of unpaired electrons.
It is possible to set non-integer positive value.
Note that the specified spin multiplicity in Geometry section (see section
1.4.1) is discarded except for zero unpaired electrons.



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

External field effectPOINTCHARGE (logical) {FALSE}
Enables to include point charges.
Coordinates and charge value are set after atomic coordinates startingfrom the character “PC” (see section 1.4.1).
6

Contents



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



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).

DFTB+U correction

PLUSU (logical) {FALSE}
Enables to add DFTB+U correction [13].
U–J parameters have to be stored in different file for each atom (seesections 1.3.2 and 1.3.3).



PLUSUTYPE (integer) {1}
Method of DFTB+U correction.
 1: Fully localized limit (FLL) 2: Pseudo self-interaction correction (pSIC)

Onsite correction

ONSITE (logical) {FALSE}
Enables to add onsite correction [14].
Onsite parameters have to be stored in different file for each atom (seesections 1.3.2 and 1.3.3).



STOREBLOCKPOP (logical) {TRUE}
Enables to store block population of last geometry step in memory forthe improved initial guess in simulation with PLUSU=TRUE andONSITE=TRUE.

Electronic property evaluation


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



CM3 (logical) {FALSE}
Enables to perform CM3 charge analysis [16, 17].

7

Contents



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



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.

Additional options for only MPI/OpenMP version


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 third-order correction (THIRDFULL=TRUE
and THIRDDIAG=TRUE).
 SHELL: Based on valence shells.
This is incompatible with on-site third-order correction(THIRDDIAG=TRUE).
 ATOM: Based on atoms.
This is incompatible with orbitally resolved SCC calculation(ORSCC=TRUE).

ORSCC (logical) {FALSE}
Enables the use of orbitally resolved SCC calculation.
Different Hubbard parameters for different atomic angular momentum
shells are used instead of the use of s-shell only.
The shell resolved calculation (RUNTYPE=SHELL) is mandatory.
The DFTB2 [12] and DFTB3 [18] methods can be extended with thistreatment.



BROYPARALLEL (logical) {FALSE}
Enables to reduce memory usage for Broyden method by using MPI
communications.

8Contents



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.

GSSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for computing exponentially decaying
short-range
partof
γ-functionofperiodic
systemin
GAMMASPLINE=TRUE.

GSDOMAIN (logical) {FALSE}
Enables the use of cell-linked list (CLL) algorithm [19, 20] for computing
exponentially decaying short-range part of γ-function 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.

9

Contents

1.1.3 NCC
The keyword NCC sets whether a non-SCC calculation (NCC-DFTB or DFTB1) [21,22] is performed or not. NCC=TRUE is identical to set SCC=FALSE.
Solver

SOLVER (integer) {1}
Eigensolver for diagonalization of Hamiltonian matrix.
 1: DSYGVD
 2: DSYGV

Open-shell

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

Electronic property evaluation


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 [15].
Mulliken population analysis (PRINTMULLIKEN=TRUE) is required.



CM3 (logical) {FALSE}
Enables to perform CM3 charge analysis [16, 17].
Mulliken population analysis (PRINTMULLIKEN=TRUE) is required.
Note that the fitted parameters are taken from DFTB2.

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

10

Contents

1.1.4 DCThe keyword DC sets whether a divide-and-conquer approach [23, 24] isemployed or not. This keyword is turned on without any specification as well as
SCC keyword. DC=FALSE is necessary if DC-DFTB calculation is deactivated.
Basic settings


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.
 ATOM defines one atom as each subsystem.
 SPECIFY defines subsystems one by one for each atom (see section






1.4.3).
SYSTEM defines all atoms as one subsystem.
AUTO constructs subsystems automatically by dividing total systemwith the specified cubic or rectangular grid.
Default is 5.0 Å cubic box.
AUTOXH constructs subsystems similarly to SUBTYPE=AUTO
without cutting covalent bonds involving hydrogen atom.
SEMIAUTO defines subsystems as the hybrid of AUTO and SPECIFY(see section 1.4.3).

 SEMIAUTOXH
constructs
subsystems
similarly
to
SUBTYPE=SEMIAUTO without cutting covalent bonds involvinghydrogen atom.
Advanced options


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.
11

Contents



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.



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



TRANSCOORD (logical) {TRUE}
Transformation of Cartesian
SUBTYPE=SEMIAUTO,

coordinate in SUBTYPE=AUTO,SUBTYPE=AUTOXH,
andSUBTYPE=SEMIAUTOXH.
 TRUE: Standard orientation is used.
 FALSE: Input geometry is used as it is.


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

Additional options for only MPI/OpenMP version


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



DOMAIN (logical) {FALSE}
Enables the use of CLL algorithm [19, 20] for accelerating assignment of
subsystems to localization 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.
12

Contents



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 [19, 20] for judging the presence ofcovalent 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}
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.

13

Contents

1.1.5 DISP/DISPERSIONThe keyword DISP/DISPERSION sets whether the noncovalent interaction(dispersion, hydrogen bond, and halogen bond) is empirically corrected or not.
Basic settings


DISPTYPE (integer) {1}
Method of noncovalent interaction correction.
 1: Grimme’s DFT-D2 [27]
 2: Slater-Kirkwood [28]The required dispersion constants have to be stored in different filefor each element (see sections 1.3.2 and 1.3.3).
 3: Lennard-Jones [29]The default values of Lennard-Jones parameters are taken fromvalues of UFF force field [30], while it is also possible to specify them
manually in different file for each element (READUFF=TRUE).
 4: Grimme’s DFT-D3 (zero-damping) [31]
 5: Grimme’s DFT-D3(BJ) (Becke and Johnson damping) [31, 32]
 6: Řezáč’s D3H4 [33]Since this corrects hydrogen bonding interaction empirically, damped
γ-function for X–H pairs may not be combined together(DAMPXH=FALSE in SCC keyword). See, for example, reference [34].
 7: DFT-ulg [35]
 8: dDMC [36]
 9: Řezáč’s D3H5 [37]Since H5 correction is designated to modify the original γ-function,
namely DAMPXH=FALSE needs to be used in SCC keyword.
This is incompatible with non-SCC calculations (NCC=TRUE or
SCC=FALSE).
 10: Kubillus’s D3X [38]14Contents

Advanced options for Grimme’s DFT-D2 (DISPTYPE=1)


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

Advanced options for Slater-Kirkwood (DISPTYPE=2)


UPDATESKPARAMS (logical) {FALSE}
Update Slater-Kirkwood dispersion parameters in each DFTB calculation.
In the case of UPDATESKPARAMS=FALSE, dispersion parameters are
kept fixed with those of initial geometry during the entire calculation.

Advanced options for Lennard-Jones (DISPTYPE=3)


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

Advanced options for Grimme’s DFT-D3 (DISPTYPE=4, 5)
ZEROS8 (real) {6.73e−01}
Global scaling factor for r−8 term 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 for r−8 term 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}
Enables to Include three-body term for Grimme’s DFT-D3 and DFT-D3(BJ).

15

Contents

Advanced options for Řezáč’s D3H4 (DISPTYPE=6)


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).



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).

16

Contents

Advanced options for DFT-ulg (DISPTYPE=7)


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.

Advanced options for dDMC (DISPTYPE=8)


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



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



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

Advanced options for Řezáč’s D3H5 (DISPTYPE=9)


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).



D3H5RS6 (real) {1.25e+00}
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).



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

Contents



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).



D3H5KFH (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).

Advanced options for Kubillus’s D3X (DISPTYPE=10)


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).
18

Contents



D3XDNI (real) {1.521e+00}
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).

Advanced options for frequency calculation


D3STEPSIZE (real) {1.0e−04}
Step size for numerical differentiation of analytically calculated Grimme’s
DFT-D3 first derivatives (DISPTYPE=4–6, 9, 10) in harmonic vibrationalfrequency analysis (FREQTYPE=1 in FREQ/FREQUENCY keyword) in
atomic unit.

Additional options for only MPI/OpenMP version


H4IDOMAIN (logical) {FALSE}
Enables the use of CLL algorithm [19, 20] for index of hydrogen atoms foraccelerating Ř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.

19

Contents



H4VDOMAIN (logical) {FALSE}
Enables the use of CLL algorithm [19, 20] 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.



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.

20

Contents

1.1.6 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 are detected in Geometry section even if PBC=TRUE is not explicitly
written (see section 1.4.1).
Basic settings
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 [39]See reference [10] for details.
There are several limitations at present.
 Stress tensor calculation is not supported (i.e. STRESS=FALSE).
 Point
chargesshouldnotbepresent(i.e.
POINTCHARGE=FALSE in SCC keyword).
 Atom resolved energies cannot be evaluated (i.e.
PRINTATOME=FALSE in MISC keyword).
 Shell resolved calculation in MPI/OpenMP version is notsupported (i.e. RUNTYPE=ATOM in SCC keyword).
 Only cubic unit cell has been tested.

Advanced options for Ewald summation (COULOMBTYPE=1)


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).
Otherwise, value is automatically determined (default and recommended).



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

21

Contents

Advanced options for multipole expansion (COULOMBTYPE=2)


FMMP (integer) {8}
Expansion order for multipole expansion.



FMMN (integer) {3}
Number of aggregated units at each level of telescoping procedure in
multipole expansion.



FMMK (integer) {5}
Maximum number of level for multipole expansion.

Additional options for only MPI/OpenMP version


EWALDSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for computation of reciprocal spacesummation in Ewald summation.
Note that the computational cost is increased.



EWALDMEMORYGB (real) {8.0e+00}
Maximum memory size for temporary arrays used in reciprocal spacesummation in Ewald summation in gigabyte (GB) unit.
If the estimated memory size for temporary arrays is larger than thespecified value, the reduced memory usage approach is used incomputations irrespective of EWALDSAVEMEMORY=FALSE.

FMMSAVEMEMORY (logical) {FALSE}
Enables to reduce memory usage for computation of solid harmonics inmultipole expansion.
Note that the computational cost is increased.

22

Contents

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

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



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



OPTTYPE (integer) {1}
Method for geometry optimization.
 1: Broyden-Fletcher-Goldfarb-Shanno (BFGS)
This uses covalent radius of atoms to construct initial Hessian matrix,
which are taken from reference [25].
 2: Steepest descent 3: Conjugate gradient 4: Quick-min [40]
 5: Fast inertial relaxation engine (FIRE) [41]MAXMOVE (real) {3.0e−01}
Maximum movable size of atomic relaxation for OPTTYPE=2–5 in Å unit.



PRINTCOORD (logical) {TRUE}
Enables to print molecular coordinate at each geometry step in thestandard output file.
In the case of PRINTCOORD=FALSE, only the final molecular coordinate
is printed in the standard output file.



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

23

Contents

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

Lattice optimization


LATTICEOPT (logical) {FALSE}
Enables to optimize translation vectors.
It is not allowed to impose partial geometry constraints(CONSTRAINT=TRUE) and to set external point charges(POINTCHARGE=TRUE in SCC keyword).
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 oftranslation 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 oftranslation 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 oftranslation vectors with FIXANGLE=TRUE.
24

Contents



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
z-coordinate for all atoms (see section 1.4.2).

Acceleration of SCC convergenceLSMULLIKEN (logical) {FALSE}
Enables to predict better initial Mulliken charges using least squares
fitting of geometric information to accelerate SCC convergence [42, 43].
This is not supported for block population used in DFTB+U correction(PLUSU=TRUE in SCC keyword) and onsite correction (ONSITE=TRUE
in SCC keyword).



LSGEN (integer) {5}
Maximum number of geometric steps to consider for least squares fittingin LSMULLIKEN=TRUE.

LSFITGEOM (logical) {FALSE}
Enables to superimpose Cartesian coordinates of previous geometricsteps on those of geometric step to be predicted in LSMULLIKEN=TRUE.
This is available for only cluster models and unnecessary to be activated.

Advanced options for BFGS (OPTTYPE=1)

BFGSSOLVER (integer) {1}
Eigensolver for diagonalization of Hessian matrix in BFGS method.
 1: DSYEVD
 2: DSYEVAdvanced options for steepest descent (OPTTYPE=2)

STEPSIZE/STEPSTDC (real) {1.0e+02}
Step size along the forces for steepest descent method in atomic unit.

Advanced options for conjugate gradient (OPTTYPE=3)


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

Contents



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

Advanced options for Quick-min (OPTTYPE=4)

DELTATQMIN/DELTASTEPQMIN (real) {1.0e−15}
Time interval between two Quick-min steps in second (s) unit.

Advanced options for FIRE (OPTTYPE=5)

DELTATFIRE/DELTASTEPFIRE (real) {1.0e−15}
Time interval between two FIRE steps in s unit.



FIREMAXDELTAT/FIREMAXDELTASTEP (real) {1.0e−14}
The largest time interval used in FIRE method in s unit.



FIRENMIN (integer) {5}Parameter for FIRE method.

FIREFINC (real) {1.1e+00}
Parameter for FIRE method.

FIREFDEC (real) {5.0e−01}
Parameter for FIRE method.

FIREALPHASTART (real) {1.0e−01}
Parameter for FIRE method.

FIREFALPHA (real) {9.9e−01}
Parameter for FIRE method.

26

Contents

1.1.8 FREQ/FREQUENCY
The keyword FREQ/FREQUENCY sets whether the harmonic vibrational frequency
analysis is performed or not. If this keyword is turned on together with geometry
optimization (OPT=TRUE or OPTIMIZE=TRUE), frequency calculation is carriedout for the last optimization step of given structures. The frequency calculation iscurrently available for cluster models.
Basic settings


FREQTYPE (integer) {1}
Method of Hessian matrix construction.
 1: Analytical second derivatives
This is not implemented for the following cases.
 Open-shell system is calculated.
 Point charges are included (POINTCHARGE=TRUE in SCC
keyword).
 Spline interpolation for exponentially decaying short-range partof γ-function is used (GAMMASPLINE=TRUE in SCC keyword).
 DFTB+U correction is used (PLUSU=TRUE in SCC keyword).
 Onsite correction is used (ONSITE=TRUE in SCC keyword).
 Shell resolved calculation is chosen in MPI/OpenMP version(RUNTYPE=SHELL, SHELL0 and ORSCC=TRUE in SCC
keyword).
 2: Numerical second derivatives using first derivatives



STEPSIZE (real) {1.0e−04}
Step size for numerical differentiation of first derivatives in FREQTYPE=2
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.



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

27

Contents



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”(see section 1.5.2 for the format).

Advanced options
HESSSOLVER (integer) {1}
Eigensolver for diagonalization of Hessian matrix.
 1: DSYEVD
 2: DSYEV

CPCONV (real) {1.0e−06}
Convergence criteria for coupled-perturbed DFTB equations forcomputation of analytical Hessian in SCC calculation in atomic unit.



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



PROJECTION (logical) {TRUE}
Projection out translational and rotational modes from normal coordinate
analysis.
In the case of PROJECTION=FALSE, full Hessian matrix is diagonalized.



READISOTOPE (logical) {FALSE}
Enables to specify arbitrary atomic mass (see section 1.4.2 details).



THERMOTEMP (real) {2.9815e+02}
Target temperature for themochemical analysis in Kelvin (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.
28

Contents



INTERRUPTCP (logical) {FALSE}
Enables to interrupt solving coupled-perturbed DFTB equations incomputation of analytical Hessian in SCC calculation.
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 beconstructed in INTERRUPTCP=TRUE.

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

INTERRUPTSEMINUM (logical) {FALSE}
Enables to interrupt performing numerical second derivatives calculationin FREQTYPE=2.
The partially constructed Hessian matrix is stored in the ASCII file
“hess.dat” for restarting purpose.



SEMINUMSTART (integer) {1}
Atom index to specify the beginning of partial Hessian matrix to beconstructed in INTERRUPTSEMINUM=TRUE.

SEMINUMEND (integer) {1}
Atom index to specify the end of partial Hessian matrix to be constructedin INTERRUPTSEMINUM=TRUE.

Additional options for only MPI/OpenMP version


SAVEMEMORYANA (logical) {FALSE}
Enables to reduce memory usage for computation of analytical Hessianin SCC calculation.
Note that the computational cost is increased.

29

Contents

1.1.9 MDThe keyword MD sets whether the molecular dynamics (MD) using velocity Verlet
integration scheme [44] is performed or not.
Basic settings


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



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



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



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



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



READVELOCITY (logical) {FALSE}
Enables to read velocities for all atoms from the ASCII file.
 TRUE: Read the file “veloc.dat” (see section 1.5.2 for the format).
 FALSE: No file to be read.

Constant energy ensemble

NVE (logical) {TRUE}
Enables to run MD with NVE ensemble.
In the case of NVE=FALSE, NVT ensemble calculation is performed.



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

30

Contents

Constant temperature ensemble

NVT (logical) {FALSE}
Enables to run MD with NVT ensemble.
In the case of NVT=FALSE, NVE ensemble calculation is performed.



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



NVTTYPE (integer) {3}
Method to control temperature for achieving constant temperatureensemble.
 1: Velocity scaling [45]
 2: Nosé-Hoover chain [46]
 3: Nosé-Hoover chain [46]In comparison to NVTTYPE=2, NVTTYPE=3 allows to control somethermostat settings.
 4: Berendsen [47]
 5: Andersen [48]

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



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



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



CALCPRESSURE (logical) {FALSE}
Enables to calculate virial pressure in constant energy or constanttemperature ensemble.
This is available in only simulations of periodic system.
31

Contents

Advanced options for thermostat


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



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



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).



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 (fs is femtosecond).



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.

Constant pressure ensemble

NPH (logical) {FALSE}
Enables to run MD with NPH ensemble.
In the case of NPH=FALSE, NVE ensemble calculation is performed.



NPT (logical) {FALSE}
Enables to run MD with NPT ensemble.
In the case of NPT=FALSE, NVE ensemble calculation is performed.

32

Contents



EXTPRESSURE (real) {1.01325e+05}
Target pressure during constant pressure ensemble simulation in Pascal(Pa) unit.



NPHTYPE (integer) {1}
Method to control pressure for achieving constant pressure ensemble.
 1: Berendsen [47]

NPTTYPE (integer) {1}
Method to control pressure and temperature for achieving (isotropic)
NPT ensemble.
 1: Hoover-Andersen [46]ISOSCALING (logical) {FALSE}
Enables to scale translation vectors uniformly during constant pressure
ensemble simulation.
Default is anisotropic change of unit cell.

Advanced options for barostat


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



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

Random number generationRANDTYPE (integer) {1}
Method to generate random numbers.
 1: Fortran intrinsic random number generator
Different platform may generate different random numbers using thesame seed.
 2: Random number generator taken from reference [49].

33Contents



SEEDTYPE (integer) {1}
Method of random number generation for velocities.
 1: No specification (the same initial velocities may be reproduced on
the same 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
SEEDTYPE=3.

for initial velocities using

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



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



SOFTCENTERTYPE (method) {ATOM}
Method to specify center of soft potential.
 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 thecenterofsoft

potential

used

in

SOFTCENTERTYPE=ATOM.
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-component of the center of soft
SOFTCENTERTYPE=SPECIFY in Å unit.

34

potential

used

in

Contents



SOFTPOINTY (real) {0.0e+00}
y-component of the center of soft
SOFTCENTERTYPE=SPECIFY in Å unit.

potential

used

in

SOFTPOINTZ (real) {0.0e+00}
z-component of the center of soft
SOFTCENTERTYPE=SPECIFY in Å unit.

potential

used

in

SOFTSHAPETYPE (method) {SPHERE}
Method to specify shape of soft potential.
 SPHERE is sphere.
 SQUARE is cubic or rectangular box.
 ELLIPSOID is ellipsoid.



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



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=ELLIPSOIDin Å unit.



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

35

Contents



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

Constraint MD with RATTLERATTLE (method) {NONE}
Constraint MD simulation using RATTLE algorithm [50] for bond lengths.
 NONE corresponds to no constraints of bond lengths.
 COVALENTXH fixes X–H covalent bonds as the sum of their covalentradii, which are taken from reference [25].
 INITGEOMXH fixes X–H covalent bond distance as the value of initialstructure.
 SPECIFY sets any bond constraints using the file “rattle.dat” (seesection 1.5.2 for the format).



RATTLEXH (real) {1.2e+00}
Threshold of X–H pair as the covalent bond in RATTLE=COVALENTXH
and RATTLE=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.

Acceleration of SCC convergence

LIMULLIKEN (logical) {FALSE}
Enables to predict better initial Mulliken charges using Lagrange
interpolation to accelerate SCC convergence [43, 51].
This is not supported for block population used in DFTB+U correction(PLUSU=TRUE in SCC keyword) and onsite correction (ONSITE=TRUE
in SCC keyword).



LIDEGREE (integer) {0}Degree of Lagrange polynomial in LIMULLIKEN=TRUE.
Default uses the previous converged charges for initial guess(conventional approach).
36

Contents



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



LSMULLIKEN (logical) {FALSE}
Enables to predict better initial Mulliken charges using least squares
fitting of geometric information to accelerate SCC convergence [42, 43].
This is not supported for block population used in DFTB+U correction(PLUSU=TRUE in SCC keyword) and onsite correction (ONSITE=TRUE
in SCC keyword).



LSGEN (integer) {5}
Maximum number of geometric steps to consider for least squares fittingin LSMULLIKEN=TRUE.

LSFITGEOM (logical) {FALSE}
Enables to superimpose Cartesian coordinates of previous geometricsteps on those of geometric step to be predicted in LSMULLIKEN=TRUE.
This is available for only cluster models and unnecessary to be activated.

Annealing simulation


ANNEAL (logical) {FALSE}
Enables to perform annealing simulation.
Temperature scheduling is specified using the file “anneal.dat” up to 99
profiles (see section 1.5.2 for the format).

Metadynamics


METAD/METADYNAMICS (logical) {FALSE}
Enables to perform metadynamics simulation [52].
Information of collective variables (CVs) is specified using the file“metacv.dat” (see section 1.5.2 for the format and the implemented CVs).

METAMAXCV (integer) {10}Maximum number of CVs used in metadynamics.



METAFREQ (integer) {100}Interval to deposit Gaussian bias potential during metadynamics.
37

Contents



METAMAXGAUSS (integer) {10000}
Maximum number of Gaussian bias potentials to be deposited during
metadynamics.



METAHEIGHT (real) {4.0e−04}
Height of Gaussian to be deposited in Hartree unit.



METARESTART (logical) {FALSE}
Enables to read history of Gaussian bias potential for restarting
metadynamics from the binary file “restart_meta”.



METAREADBIAS (logical) {FALSE}
Enables to read history of Gaussian bias potential from the ASCII file.
 TRUE: Read the file “bias.dat” (see section 1.5.2 for the format).
 FALSE: No file to be read.

Advanced options for metadynamics


WTMETAD/WTMETADYNAMICS (logical) {FALSE}
Enables to perform well-tempered metadynamics simulation [53].



WTMETABIASTEMP (real) {1.2e+03}
Bias temperature of well-tempered metadynamics in K unit.



METAPBCMIC (logical) {TRUE}
Enables to apply minimum image convention to atom-pair distances forcalculation of CVs in metadynamics of periodic system.
Note that minimum image convention is not used to calculate center ofcoordinate or center of mass and some CVs (see section 1.5.2 for details).



METAPRINTFES (logical) {TRUE}
Enables to write free energy surface to the file “fes.dat” every depositionof Gaussian bias potential.
Information of grid must be specified in “metacv.dat” (see section 1.5.2
for details).
Note that scaling factor for the free energy surface is included in the caseof well-tempered metadynamics [53].

38Contents



METAWALL (logical) {FALSE}
Enables to apply wall potential to constrain sampling range of CVs in
metadynamics.
Information of wall potential must be specified in “metacv.dat” (seesection 1.5.2 for details).

39

Contents

1.1.10 MISC
The keyword MISC sets other miscellaneous options.
Additional analysis and property evaluation


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



PRINTRESOLVEDFORCE (logical) {FALSE}
Enables to evaluate force contributions (electronic, repulsive, anddispersion term if included).
If stress tensor calculation is performed (STRESS=TRUE in PBC
keyword), corresponding contributions of stress tensor are additionallyevaluated.
Note that total forces are printed in the standard output while therespective 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.



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



PRINTATOME (logical) {FALSE}
Enables to evaluate atomic energy contributions (electronic, repulsive,and dispersion term if included).
Note that the respective contributions are printed in the detailed data file.

40

Contents

Global options


FORCETYPE (integer) {1}
Method of gradient calculation.
 1: Analytical first derivatives.
 2: Numerical first derivatives using DFTB total energy.



FORCESTEPSIZE (real) {1.0e−04}
Step size for numerical differentiation of DFTB total energy in
FORCETYPE=2 in atomic unit.



HSSTEPSIZE (real) {1.0e−02}
Step size for numerical differentiation of DFTB parametrized integrals inanalytical first and second derivatives calculations in atomic unit.



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 andwrites the 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 filecontains 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.



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.



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

Contents



PRINTHS (logical) {FALSE}
Enables to print the lower triangle of non-perturbed Hamiltonian andoverlap 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.

Additional options for only serial and OpenMP versions


USESCRATCH (logical) {FALSE}
Enables to use disk I/O for storing some intermediate data.
In the case of USESCRATCH=FALSE, all intermediate data are stored inmemory.

42

Contents

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. All comments will be printed in the standard output file as it is (seesection 3.1).

43

Contents

1.3 Parameter section
1.3.1 General overview
2

Number of unique atom types

O˽2

Atom symbol and its highest angular momentum for valence shells

oo.spl˽oh.spl

Slater-Koster parameters: atom1–atom1 and atom1–atom2

H˽1

Atom symbol and its highest angular momentum for valence shells

ho.spl˽hh.spl

Slater-Koster parameters: atom2–atom1 and atom2–atom2

This section starts from the line giving the number of atom types in the given
structure. The present example is 2 since water molecule contains oxygen andhydrogen. Next, atom symbol and its maximum angular momentum are specified
for the first atom type (atom1). The integer 1, 2, 3, 4 corresponds to the maximum
angular momentum s, p, d, f, respectively. Then, the paths to Slater-Koster
parameter files including the first atom types are specified in a single line. The
order has to be consistent in the order of atom types specified in this section,
namely atom1 with atom1, atom1 with atom2 (the second atom type), atom1 with
atom3 (the third atom type), … These two lines are repeated until all atom types
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 thesame directory as the one containing the input file. When the prefix of current
directory is explicitly included, this section looks like as follows.
2
O˽2
./oo.spl˽./oh.spl
H˽1
./ho.spl˽./hh.spl

For treating f-electrons, please ensure that relevant Slater-Koster files have theextended format.

44

Contents

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 after the highest angularmomentum for valence shells
This is necessary for Hubbard derivatives in third order corrections(THIRDFULL=TRUE or 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 MPI/OpenMP version, values of s-, p-, d-, and f-shells need
to be provided separately.

O˽2˽−0.1575

Atom resolved case (usual)Cu˽3˽−0.20˽−0.0575˽−0.0575

Orbitally resolved case (order is (f), d, p, s)

Preparation of additional ASCII files
This is necessary for the following five options: Spin constants for spin polarized calculation (SPIN=TRUE in SCC
keyword) U–J parameters for DFTB+U correction (PLUSU=TRUE in SCC
keyword) Onsite parameters for onsite correction (ONSITE=TRUE in SCC
keyword) Atomic polarizabilities, van der Waals radii, and effective charge forSlater-Kirkwood dispersion (DISPTYPE=2 in DISP/DISPERSION
keyword) Lennard-Jones
parametersforLennard-Jones
dispersion(DISPTYPE=3 and READUFF=TRUE in DISP/DISPERSION keyword)
The path of additional file should be specified after the maximum angularmomentum (or Hubbard derivatives if they exist). In the case of specifying
multiple files, the path of dispersion constants file should come first,
followed by the path of spin constants file, subsequently the path of U–J
parameter file, and finally the path of onsite parameter file.

45

Contents

Examples for oxygen atom are given below.
O˽2˽o.wll

Only spin constants fileO˽2˽o.dispsk

Only dispersion constants fileO˽2˽o.uj

Only U–J parameters fileO˽2˽o.oc

Only onsite parameters fileO˽2˽o.dispsk˽o.wll

Both spin and dispersion constants filesO˽2˽o.dispsk˽o.uj

Both dispersion constants and U–J parameters
filesO˽2˽o.dispsk˽o.oc

Both dispersion constants and onsite parameters
filesO˽2˽o.wll˽o.uj

Both spin constants and U–J parameters filesO˽2˽o.wll˽o.oc

Both spin constants and onsite parameters filesO˽2˽o.uj˽o.oc

Both U–J parameters and onsite parameters filesO˽2˽o.dispsk˽o.wll˽o.uj

Dispersion constants, spin constants and U–J
parameters filesO˽2˽o.dispsk˽o.wll˽o.oc

Dispersion constants, spin constants and onsite
parameters filesO˽2˽o.dispsk˽o.ud˽o.oc

Dispersion constants, U–J parameters and onsite
parameters filesO˽2˽o.wll˽o.ud˽o.oc

Spin constants, U–J parameters and onsite
parameters filesO˽2˽o.dispsk˽o.wll˽o.ud˽o.oc

All dispersion constants, Spin constants, U–J
parameters and onsite parameters files

46

Contents

1.3.3 Format of ASCII files for additional constants
The ASCII files for specifying additional constants described in section 1.3.2 havedifferent format in each case. Here examples for carbon atom are shown.
Spin constants file−0.031
−0.025
−0.025
−0.023

The number of lines of spin constants file depends on the maximum angularmomentum (1, 4, 9, 16 for s, p, d, f). Each line contains the atomic spin constantin atomic unit in the following order:





s: ss
p: ss, sp, ps, pp
d: ss, sp, sd, ps, pp, pd, ds, dp, dd
f: ss, sp, sd, sf, ps, pp, pd, pf, ds, dp, dd, df, fs, fp, fd, ff

For reference, PBE values for H, C, N, O, S, Fe, and Ni listed in DFTB+ manual [54]are stored into the file with suffix “.wll” in the slko directory.
U–J parameters file0.178371875˽1
0.178371875˽1

The number of lines of U–J parameters file depends on the maximum angularmomentum (1, 2, 3, 4 for s, p, d, f). There should exist two columns in each line.
The first column specifies the magnitude of effective on-site Coulomb interaction
in atomic unit and the second column specifies the index of block for valenceshells. For example, if s- and p-orbitals needs to be treated as separated block,the second column of second line will be changed from “1” to “2”. Note that thisfile should be prepared with zero value for the first column even if no DFTB+U
correction is applied for some atom types.

47

Contents

Onsite parameters file−0.00000˽0.00000
−0.04973˽0.10512
−0.04973˽0.10512
−0.01203˽0.02643

The number of lines of onsite parameters file depends on the maximum angularmomentum (1, 4, 9, 16 for s, p, d, f). There should exist two columns in each line.
The first and second columns specify the same and different spin parameters inatomic unit in the following order:



s: ss
p: ss, sp, ps, pp

d: ss, sp, sd, ps, pp, pd, ds, dp, dd
f: ss, sp, sd, sf, ps, pp, pd, pf, ds, dp, dd, df, fs, fp, fd, ff

For reference, PBE values for H, C, N, O, and S listed in reference [14] or in DFTB+
manual [54] are stored into the file with suffix “.oc” in the slko directory.
Slater-Kirkwood dispersion constants file
−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 threehydrogen neighbors. For reference, H, C, N, O, and S dispersion constants listed
in DFTB+ manual [54] are stored into the file with suffix “.dispsk” in the slko
directory.
Lennard-Jones dispersion parameters file
−3.851˽0.105

There should exist two columns in a single line. The first column specifies nonbonded distance in Å unit and the second column specifies non-bonded energy
in kcal/mol unit. For reference, H, C, N, O, S, Fe, and Ni Lennard-Jones parameterslisted in reference [30] are stored into the file with suffix “.displj” in the slko
directory.
48

Contents

1.4 Geometry section
1.4.1 General overview
Number of atoms,

3˽0˽1

charge,spin multiplicity
O0˽−0.06275391˽−0.06275391˽−0.00000000

Cartesian coordinatesof atom label 1

H0˽−1.01910555˽−0.08185946˽−0.00000000

Cartesian coordinatesof atom label 2H0˽−0.08185946˽−1.01910555˽−0.00000000

Cartesian coordinatesof atom 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 threetranslation vectors with label 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

49

Contents

If the external static electric field is included (EFIELD=TRUE in SCC keyword), theelectric field vector can be specified by adding one line with label EF after the lastCartesian coordinates (EF should be an upper-case). The example of watermolecule 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
H0˽−0.08185946˽−1.01910555˽−0.00000000
EF˽−0.00100000˽−0.00100000˽−0.00000000

If the point charges are included (POINTCHARGE=TRUE in SCC keyword), thecoordinates and charge value can be specified by adding lines with label PC after
the last Cartesian coordinates (PC should be an upper-case). The example of watermolecule 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

50

Contents

1.4.2 Specification of additional information
Additional columns after Cartesian coordinates can be read when you would liketo specify atoms to be moved or atomic masses to be changed to investigate theisotope 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/OPTIMIZE 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 vibrationalfrequency analysis (READISOTOPE=TRUE in FREQ/FREQUENCY
keyword)
The form “ISO=value” changes atomic mass from the default value in
Slater-Koster parameter file 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
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 withatom label 2 can be moved along only the y-axis and the coordinate of hydrogenwith atom label 3 is never updated. If FREQ=(READISOTOPE=TRUE) is
additionally specified, the atomic 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.

51

Contents

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. Afterthe blank line of Geometry section, the indices of subsystems for each atom arelisted.
Here is an example of specification of 5 subsystems for the given structure
containing 22 atoms in SUBTYPE=SPECIFY:
1˽1˽1˽1˽1

Index of subsystem for atom label 1, 2, 3, 4, 52˽2˽2˽2

Index of subsystem for atom label 6, 7, 8, 93˽3˽3˽3

Index of subsystem for atom label 10, 11, 12, 13

4˽4˽4˽4

Index of subsystem for atom label 14, 15, 16, 17

5˽5˽5˽5˽5

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
to subsystem 2, atoms with label 10−13 belong to subsystem 3, and so on. It is
not necessarily 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 numberof atoms in the given structure.
In the case of SUBTYPE=SEMIAUTO, the integer “0” is used as an index foratoms to be 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, 50˽0˽0˽0

Index of subsystem for atom label 6, 7, 8, 90˽0˽0˽0

Index of subsystem for atom label 10, 11, 12, 130˽0˽0˽0

Index of subsystem for atom label 14, 15, 16, 172˽2˽2˽2˽2

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 to subsystem 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.

52

Contents

1.5 Auxiliary input files
In addition to the main input file, Slater-Koster parameter files, and ASCII files for
additional constants if any (see section 1.3.3), some auxiliary input files arenecessary depending on computational setup. Each auxiliary input file has the
specific filename and should be located in the same directory as the one
containing the main input file. There are two types of auxiliary input files: (i) using
binary file dumped in previous run of program and (ii) preparing ASCII file withspecific format.
1.5.1 Binary format files
Binary format files serve to providing information for restarting purpose. See
section 3.3 or relevant options in the description for the contents of respective
file. Note that the binary file should not be used with different versions (DO NOTMIX: (i) OpenMP (serial) and MPI/OpenMP and (ii) version 1.0 and 2.0!).
“chrgfile”
This is required if SCC=(READCHARGE=BINARY) is specified.
Please make sure that the previous run of program specifies SCC=TRUE
and WRITECHRGFILE=TRUE in MISC keyword.

“restart_hess”
This is required if FREQ=(RESTARTHESS=BINARY) is specified.
Please make sure that the previous run of program specifies
PRINTHESS=TRUE in FREQ/FREQUENCY keyword.

“restart”
This is required if MD=(RESTART=TRUE) is specified.
Not only “restart” itself, “restart_chk” may be used after renaming.

“restart_meta”
This is required if MD=(METARESTART=TRUE) is specified.
Not only “restart_meta” itself, “restart_meta_chk” may be used after
renaming.
Please make sure that the previous and present run of program specifies
METAD=TRUE, METADYNAMICS=TRUE, WTMETAD=TRUE, or
WTMETADYNAMICS=TRUE in MD keyword.
53

Contents

1.5.2 ASCII format files
ASCII format files mainly serve to providing detailed information which depends
on the given structure.
Reading initial charge in SCC calculation
“charge.dat”
This is required if SCC=(READCHARGE=ASCII) is specified.
Please make sure that the previous run of program specifies SCC=TRUE.
Examples for water molecule are given below, where oxygen corresponds
to atom label 1 and hydrogen corresponds to atom labels 2 and 3.
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
−0.6

Mulliken population of O (atom label 1) is 6.6−0.3

Mulliken population of H (atom label 2) is 0.7−0.3

Mulliken population of H (atom label 3) is 0.7

 ATOMSPECIES lists the atom symbol and the net charge in the firstand 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

H˽−0.3

Mulliken population of H (atom labels 2 and 3) is 0.7

54Contents

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 s shell of O (atom label 1) is 1.8−0.8

Mulliken population of p shell of O (atom label 1) is 4.8−0.3

Mulliken population of s shell of H (atom label 2) is 0.7−0.3

Mulliken population of s shell of H (atom label 3) is 0.7

 SHELLSPECIES lists the atom symbol, subshell symbol (s/p/d/f), andthe 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 symbol (a/b) is specified in the fourth column.
SHELLSPECIESO˽s˽−0.2

Mulliken population of s shell of O (atom label 1) is 1.8O˽p˽−0.8

Mulliken population of p shell of O (atom label 1) is 4.8H˽s˽−0.3

Mulliken population of s shell of H (atom labels 2 and
3) is 0.7

Reading Hessian matrix in frequency calculation
“hess.dat”
This is required if FREQ=(RESTARTHESS=ASCII) is specified.
Please make sure that the previous run of program specifies
INTERRUPTCP=TRUE
or
INTERRUPTSEMINUM=TRUE
in
FREQTYPE=1 or FREQTYPE=2 in FREQ/FREQUENCY keyword.
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 elementof Hessian matrix. It is enough for off-diagonal elements to set once inthe file. Namely, the program assumes reading 3N(3N + 1)/2 lines for the
given structure containing N atoms.
1˽1˽−0.000171

Hessian matrix element of (H)11 is −0.000171

2˽1˽−0.000000

Hessian matrix element of (H)21 and (H)12 is 0.000000

3˽1˽−0.000083

Hessian matrix element of (H)31 and (H)13 is 0.000083…

55

Contents

Reading initial velocity in MD simulation“veloc.dat”
This is required if MD=(READVELOCITY=TRUE) is specified.
Each line consists of three columns: x, y, z velocity components are given
for each atom of the given structure in Geometry section in Å/fs unit. (DONOT INCLUDE ATOM SYMBOL!)−0.00154193˽−0.00142221˽−0.00058536

Velocity components ofatom label 1−0.02998618˽−0.02026292˽−0.00706185

Velocity components ofatom label 2−0.00549584˽−0.00232603˽−0.00223541

Velocity components ofatom label 3…

Manual construction of bond constraints in RATTLE
“rattle.dat”
This is required if MD=(RATTLE=SPECIFY) is specified.
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
bond distance is fixed. The program assumes no blank line until the endof file.
 C: Bond distance is fixed to the sum of covalent radii, which are takenfrom reference [25].
 I: Bond distance is fixed to the value of initial structure.
 S: Bond distance is fixed to the value specified in the fourth columnin Å unit.
Note that the sum of covalent radii is used if no fourth column is set.
1˽2˽C

Bond distance of atom label 1 and 2 is the sum of
covalent radii

2˽3˽I

Bond distance of atom label 2 and 3 is the value of
initial structure

1˽3˽S˽0.95

Bond distance of atom label 1 and 3 is 0.95 Å

56

Contents

Setting temperature profile in annealing simulation“anneal.dat”
This is required if MD=(ANNEAL=TRUE) is specified.
The first column is an upper-case letter to determine how thetemperature is varied and the second column specifies the number of MD
steps that follows the corresponding profile. The third column istemperature value in K unit. The program assumes no blank line until theend of file.
 C: Temperature is set to constant value specified in the third column.
 L: Temperature varies linearly from the temperature of previous
profile to the value specified in the third column.
 G: Temperature varies geometrically from the temperature ofprevious profile to the value specified in the third column.
 E: Temperature varies exponentially from the temperature of previous
profile to the value specified in the third column.
The shape of exponential function is controlled by the value specified
in the fourth column.
Note that the value of 5.0 is used if no fourth column is set.
 S: Temperature varies sigmoidally from the temperature of previous
profile to the value specified in the third column.
The shape of sigmoid function is controlled by the value specified in
the fourth column.
Note that the value of 3.0 is used if no fourth column is set.
C˽1001˽0300.0

The first step is set to 300 K
INITTEMP option is updated with the specified valueL˽1000˽1000.0

Next 1000 steps increase linearly from 300 to 1000 K

C˽1000˽1000.0

Next 1000 steps keep the temperature 1000 K(equilibration)

E˽1000˽0500.0˽3.0

Next 1000 steps decrease exponentially with decay
factor of 3.0 from 1000 to 500 K

C˽1000˽0500.0

Next 1000 steps keep the temperature 500 K(equilibration)

S˽1000˽0300.0˽2.0

Next 1000 steps decrease sigmoidally with gain factorof 2.0 from 500 to 300 K

57

Contents

Setting CVs in metadynamics simulation“metacv.dat”
This is required if MD=(METAD=TRUE), MD=(METADYNAMICS=TRUE),
MD=(WTMETAD=TRUE), or MD=(WTMETADYNAMICS=TRUE) is
specified.
The first column of first row is the type of CVs, followed by the width ofdeposited Gaussian and the associated options. If free energy surface iswritten during the dynamics (METAPRINTFES=TRUE in MD keyword),the information of grid is additionally required. This is specified by adding
three columns in the order of lower bound, upper bound, and grid
spacing (unit is identical to that of Gaussian width). For some CVs, theintegers to specify atom labels or species for a group of atoms are given
in second and third rows.
If one needs to specify multiple CVs, the information of second and thesubsequent CVs is added from the next line of the end of definition ofpreceding CVs.
The implemented CVs are listed with each example below, where total
number of atoms in the system is assumed to be 10.
 ABSOLUTEPOSITION
The absolute position of an atom or a group of atoms.
For periodic system, minimum image convention is not applied even
if METAPBCMIC=TRUE in MD keyword.
ABSOLUTEPOSITION˽0.1˽1˽Z˽LABELS˽COC
5

The z-axis coordinate of atom label 5 is set to CV

where the second to sixth columns of the first row specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species in a group 4th: Direction to be restricted (X/Y/Z)
 5th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 6th: Center of a group of atoms
COC ··· Center of coordinate
COM ··· Center of mass
58

Contents

 BONDDISTANCE
The distance between two atoms (atom1–atom2).
Bond distance of atom label 1 and 2 is set

BONDDISTANCE˽0.1˽1˽2

to CV

where the second to fourth columns specify
 2nd: Gaussian width in Å unit
 3rd: Label specifying atom1
 4th: Label specifying atom2
 BONDANGLE
The angle defined by three atoms (atom1–atom2–atom3).
Bond angle of atom label 1–2–3 is set to

BONDANGLE˽0.1˽1˽2˽3

CV

where the second to fifth columns specify
 2nd: Gaussian width in degree (deg) unit
 3rd: Label specifying atom1
 4th: Label specifying atom2
 5th: Label specifying atom3
 BONDDIHEDRAL
The dihedral angle defined by four atoms (atom1–atom2–atom3–
atom4).
BONDDIHEDRAL˽0.1˽1˽2˽3˽4

Bond dihedral of atom label 1–2–3–4 is
set to CV

where the second to sixth columns specify
 2nd: Gaussian width in deg unit
 3rd: Label specifying atom1
 4th: Label specifying atom2
 5th: Label specifying atom3
 6th: Label specifying atom459

Contents

 BONDDISTANCEDIFFERENCE
The difference between two bond distances in which each bonddistance is calculated as the distance between two atoms as follows:

BONDDISTANCEDIFFERENCE˽0.1˽1˽2˽3˽4

Distance r12 − r34 is set toCV

where the second to sixth columns specify

2nd: Gaussian width in Å unit
3rd: Label specifying atom1
4th: Label specifying atom2
5th: Label specifying atom3

 6th: Label specifying atom4
 BONDDISTANCEADDITION
The summation over two bond distances in which each bond distanceis calculated as the distance between two atoms as follows:

BONDDISTANCEADDITION˽0.1˽1˽2˽3˽4

Distance r12 + r34 is set toCV

where the second to sixth columns specify
 2nd: Gaussian width in Å unit
 3rd: Label specifying atom1
 4th: Label specifying atom2
 5th: Label specifying atom3
 6th: Label specifying atom460

Contents

 MINIMUMDISTANCE
The minimum value of atom pair distances

given as follows:

MINIMUMDISTANCE˽0.1˽3˽2˽500˽LABELS1˽2˽3

Group1 contains 3 atoms with labels 1, 2, and 3

4˽5

Group2 contains 2 atoms with labels 4 and 5

where the second to sixth columns of the first row specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species in group1
 4th: Number of atoms or species in group2
 5th: Positive real parameterto control smoothness
 6th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section MEANDISTANCE
The mean distance between two groups of atoms given as follows:

MEANDISTANCE˽0.1˽3˽2˽6˽LABELS1˽2˽3

Group1 contains 3 atoms with labels 1, 2, and 3

4˽5

Group2 contains 2 atoms with labels 4 and 5

where the second to sixth columns of the first row specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species in group1
 4th: Number of atoms or species in group2
 5th: Positive integer
to specify the power of inverse of distance 6th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
61

Contents

 ATOMPOINTPLANEDISTANCE
The distance between the point defined by an atom position (atom1)
and the plane defined by three atoms (atom2, atom3, and atom4).
For periodic system, minimum image convention is not applied even
if METAPBCMIC=TRUE in MD keyword.
ATOMPOINTPLANEDISTANCE˽0.1˽1˽2˽3˽4

Distance between atomlabel 1 and plane definedby atom labels 2, 3, and 4
is set to CV

where the second to sixth columns specify
 2nd: Gaussian width in Å unit
 3rd: Label specifying atom1 in the point 4th: Label specifying atom2 in the plane 5th: Label specifying atom3 in the plane 6th: Label specifying atom4 in the plane ATOMPLANEPLANEANGLE
The angle between two planes in which each plane is defined by threeatoms (atom1, atom2, atom3 and atom4, atom5, atom6).
For periodic system, minimum image convention is not applied even
if METAPBCMIC=TRUE in MD keyword.
ATOMPLANEPLANEANGLE˽0.1˽1˽2˽3˽4˽5˽6

Angle

betweenplanesdefined by atom labels 1,
2, 3 and 4, 5, 6 is set to CV

where the second to eighth columns specify








2nd: Gaussian width in deg unit
3rd: Label specifying atom1 in the first plane
4th: Label specifying atom2 in the first plane
5th: Label specifying atom3 in the first plane
6th: Label specifying atom4 in the second plane
7th: Label specifying atom5 in the second plane
8th: Label specifying atom6 in the second plane62

Contents MIDPOINTDISTANCE
The distance between two points in which each point is calculated aseither center of coordinate or center of mass of a group of atoms.
MIDPOINTDISTANCE˽0.1˽3˽3˽LABELS˽COC
1˽2˽3

The first point is calculated using 3 atoms with labels
1, 2, and 3

4˽5˽6

The second point is calculated using 3 atoms withlabels 4, 5, and 6

where the second to sixth columns of the first row specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species for the first point 4th: Number of atoms or species for the second point 5th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 6th: Center of a group of atoms
COC ··· Center of coordinate
COM ··· Center of mass

63

Contents

 GYRATIONRADIUS
The radius of gyration of a group of atoms.
GYRATIONRADIUS˽0.1˽10˽LABELS˽COC

All 10 atoms are included
to calculate non-weighted
radius of gyration

where the second to fifth columns specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species in a group
Note that no second line if total number of atoms is specified.
 4th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 5th: Center of a group of atoms
COC ··· Center of coordinate, which calculates the non-weighted
radius of gyration given by

wheredenotes a coordinate.
COM ··· Center of mass, which calculates the mass-weighted
radius of gyration given bywhere

is an atomic mass.

64

Contents

 INERTIAMOMENTThe moment of inertia of a group of atoms.
INERTIAMOMENT˽0.1˽10˽LABELS˽COC

All 10 atoms are included
to calculate non-weighted
moment of inertia

where the second to fifth columns specify
 2nd: Gaussian width 3rd: Number of atoms or species in a group
Note that no second line if total number of atoms is specified.
 4th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 5th: Center of a group of atoms
COC ··· Center of coordinate, which calculates the non-weighted
moment of inertia given by

wheredenotes a coordinate. Note that Gaussian width isspecified with Å2 unit.
COM ··· Center of mass, which calculates the mass-weighted
moment of inertia given bywhere
is an atomic mass. Note that Gaussian width is specified
with AMU∙Å2 unit (AMU is atomic mass unit).

65

Contents

 RATIONALCOORDINATIONNUMBER
The coordination number between two groups of atoms using arational function given as follows:

RATIONALCOORDINATIONNUMBER˽0.1˽3˽2˽6˽12˽2.0˽LABELS˽AVERAGE
1˽2˽3

Group1 contains 3 atoms with labels 1, 2, and 3

4˽5

Group2 contains 2 atoms with labels 4 and 5

where the second to ninth columns of the first row specify
 2nd: Gaussian width
Note that coordination number is unitless.
 3rd: Number of atoms or species in group1
 4th: Number of atoms or species in group2
 5th: Positive integer
to specify the power appearing in
numerator
 6th: Positive integer
denominator 7th: Reference distanceto specify the power appearing inin Å unit

 8th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 9th: Normalization factor
SUM ···AVERAGE ···66

Contents

 FERMICOORDINATIONNUMBER
The coordination number between two groups of atoms using a
Fermi function given as follows:

FERMICOORDINATIONNUMBER˽0.1˽3˽2˽8.0˽2.0˽LABELS˽AVERAGE
1˽2˽3

Group1 contains 3 atoms with labels 1, 2, and 3

4˽5

Group2 contains 2 atoms with labels 4 and 5

where the second to eighth columns of the first row specify
 2nd: Gaussian widthNote that coordination number is unitless.
3rd: Number of atoms or species in group1
4th: Number of atoms or species in group2
5th: Positive real parameterto control steepness of Fermi
function
6th: Reference distancein Å unit

 7th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 8th: Normalization factor
SUM ···AVERAGE ···67

Contents

 RMSD
The root mean square deviation (RMSD) of a group of atoms.
RMSD˽0.1˽10˽LABELS˽COC˽FIT˽ref.xyz

RMSD with respect to thestructure in ref.xyz is set toCV

where the second to seventh columns specify
 2nd: Gaussian width in Å unit
 3rd: Number of atoms or species in a group
Note that no second line if total number of atoms is specified.
 4th: Method to define atoms
LABELS ··· Atom labels
SPECIES ··· Order of the listed elements in Parameter section
 5th: Center of a group of atoms
COC ··· Center of coordinate, which calculates the non-weighted
RMSD given by

wheredenotes a coordinate.
COM ··· Center of mass, which calculates the mass-weighted
RMSD given bywhere
is an atomic mass.
 6th: Method to compute RMSD
AS-IS ··· Using Cartesian components without any transformation
FIT ··· Superimposing present and reference coordinates byquaternion
This is available for only cluster models.
 7th: Filename containing reference Cartesian coordinates
This is xyz format.

68

Contents

Reading bias potential in metadynamics simulation“bias.dat”
This is required if MD=(METAREADBIAS=TRUE) is specified.
Please make sure that the previous and present run of program specifies
METAD=TRUE, METADYNAMICS=TRUE, WTMETAD=TRUE, or
WTMETADYNAMICS=TRUE in MD keyword.
The format should be identical to the output file of metadynamics
simulation, “biaspot”, which contains the history of Gaussian height,
width, and value of CVs (see also section 3.3.4).
An example for CV of bond distance is given below.
GAUSSIAN˽BIAS˽POTENTIAL:˽˽˽˽˽˽˽˽˽1***˽AT˽T=˽˽˽˽˽˽˽˽˽10.00˽FSEC,˽THIS˽RUN’S˽STEP˽NO.=˽˽˽˽˽˽20
˽˽˽˽Gaussian˽height˽˽˽˽˽=˽˽˽˽˽˽˽˽0.0040000000˽a.u.
˽˽˽˽Collective˽variable˽=˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽1˽˽˽˽˽˽Coordinate˽˽˽˽˽˽˽˽=˽˽˽˽˽˽˽˽1.6146539213˽Angstrom˽˽˽˽˽˽Gaussian˽width˽˽˽˽=˽˽˽˽˽˽˽˽0.1000000000˽Angstrom
…

The following remarks should be noted:
 Entry of Gaussian bias potential (first line) should not exceed thenumber specified in METAMAXGAUSS option.
 Information of simulation time and MD step (second line) is ignored
when the program reads this file.
 Gaussian height (third line) is always constant (value specified in
METAHEIGHT option) in standard metadynamics while it gradually
decreases in well-tempered metadynamics.
 Collective variable, coordinate, and Gaussian width (fourth, fifth, sixthlines) are printed repeatedly if multiple CVs are specified in“metacv.dat”.

69Contents

Setting wall potential in metadynamics simulation
This is required if MD=(METAWALL=TRUE) is specified.
The following even polynomial function can be imposed as wall potential:

The parameters are specified in “metacv.dat” together with direction of wall
potential.
Here is an example of applying quadratic potential if the value of first CV is greaterthan 2.0 in two-dimensional metadynamics.
<Description of first CV>

See setting CVs in metadynsmics simulation for details<Description of second CV>

See setting CVs in metadynsmics simulation for details

One blank line to separate definitions of CVs and wall
potentials

U˽1˽1.0˽2.0˽2

Description of wall potential

where the first to fifth columns of the last line specify
 1st: Direction of wall potential to be imposed

L ··· Lower bound
U ··· Upper bound2nd: Index of CV
3rd: Coefficientin kcal/(mol∙Å2) unit
4th: Limit of CV value

Note that the unit is identical to that of Gaussian width.
5th: Exponent in positive even integer

If one needs to specify multiple wall potentials, the information of second and thesubsequent wall potentials is added from the next line of the end of definition ofpreceding wall potentials.

70

Contents

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

OpenMP version<installation path>/dftb_openmp.00.x

MPI/OpenMP version

mpiexec˽–np˽(numberofprocesses)˽<installation

path>/dftb_mpiomp.00.x

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

OpenMP version<installation path>/dftb_openmp.00.x file

MPI/OpenMP version

mpiexec˽–np˽(numberofprocesses)˽<installation

path>/dftb_mpiomp.00.x fileFor OpenMP parallelization, number of threads must be set before starting the
calculation.
bash

export˽OMP_NUM_THREADS=(number of threads)

csh, tcsh

setenv˽OMP_NUM_THREADS˽(number of threads)

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

ulimit˽–s˽unlimited
export˽OMP_STACKSIZE=(stack size value)

csh, tcsh

unlimit
setenv˽OMP_STACKSIZE˽(stack size value)

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

71

Contents

3 Output of DCDFTBMD
If the calculation terminates normally, the standard output file and the detailed
data file will be ordinarily appeared. By default (using “dftb.inp” as the input file),
the corresponding filenames are “dftb.out” and “dftb.dat”, respectively. Whenthe header of input file is specified arbitrarily (e.g. “file.inp” as described in
section 2), the header of output and data 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 of water 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 is slightly different for
the information of subsystem and Mulliken population. In the following, symbols
for representing new line () and blank area (˽) are omitted.

3.1 Standard output file
The standard output file begins with the header of program. The version number
is varied depending on the released executables. After punching of start time, thecomment line is printed as described in Title section (see section 1.2).
*********************************************************************************************************************
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 2.0 (Dec 2019)*********************************************************************************************************************
Execution of DCDFTBMD begun %a %b %d %H:%M:%S %Y-------TITLE--------

72

Contents

Next, the path of Slater-Koster parameter files is printed, followed by thesummary of calculation type.
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= True------------------------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…
DC= False------------------------DISP= False------------------------PBC= False------------------------OPT= False------------------------FREQ= False------------------------MD= False------------------------MISC= False-------------------------

73

Contents

Then, the system dependent quantities are printed.
Total number of basis set shells = 4
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

Total number of atoms, charge, and spin multiplicity correspond to values
specified in Geometry section (see section 1.4.1). The numbers of basis set shells,
basis functions, and electrons are calculated for water molecule as follows: Number of basis set shells: O(s) + O(p) + H(s) + H(s) = 4
 Number of basis functions: O(s) + O(px) + O(py) + O(pz) + H(s) + H(s) = 6
 Number of electrons: O(s2p4) + H(s1) + H(s1) = 8
Note that these quantities also depend on the used Slater-Koster parameter files
and may be different from the chemical intuition considering minimal basis set
and valence electrons.
Finally, the SCC-DFTB single point calculation starts in the present example.
*** Start single point energy calculation ************************************** 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.

74

Contents

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 MPI/OpenMP version is used, this part looks concise excluding the case of
PRINTLOCALREGION=TRUE in MISC keyword.
SUBSYSTEM--------------------1

1

1

Next, the eigenvalues of the molecular orbitals in Hartree unit (first column) andthe 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

75

Contents

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-----------------------------------------------------------------------------------------1

OCC.
1

O

-0.8489865909

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

S

0 -0.0000000000

P-10.4626963766

P

00.0000000000

P

1 -0.4626963766

2

H

S

0 -0.3908041098

3

H

S

0

2.0000000000

2.0000000000

0.3908041098…

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

6.5926170634-0.5926170634

20.7036914683-0.2963085317

30.7036914683-0.2963085317

76

Contents

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

2

O

6.5926170634

S

0 1.7342141199

P-1 1.4292014718

P

0 2.0000000000

P

1 1.4292014718

H

0.7036914683
S

3-0.5926170634-0.2963085317

0 0.7036914683

H

0.7036914683
S-0.2963085317

0 0.7036914683

If MPI/OpenMP version is used, 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 are shown 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

77Contents

If MAYER=TRUE in SCC or NCC keyword, atomic valences (diagonal elements)and Mayer 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.912201

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 :

1.1372233679

z0.0000000000

1.6082767104

78

Contents

Finally, the different energy contributions to the total energy are shown in Hartree
unit. If MPI/OpenMP version is used, PRINTATOME=TRUE in MISC keyword isnecessary 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 thetotal, electronic, and repulsive energies are additionally printed in Hartree unit.
ATOM RESOLVED ENERGY CONTRIBUTION [Eh]-----------------------------------------------------------------------TOTALELECTRONICREPULSIVE

1-3.3763986686-3.4127218718-0.0363232032

2-0.3507697335-0.3689313351-0.0181616016

3-0.3507697335-0.3689313351-0.0181616016

79

Contents

3.3 Other output files
Additional output files will be generated depending on the type of calculations.
3.3.1 Output files in DFTB calculation“chrgfile”
This is generated if SCC=TRUE and MISC=(WRITECHRGFILE=TRUE).
The result of Mulliken population analysis is contained in binary format.
The file may be used for setting initial charge in READCHARGE=BINARY
in SCC keyword (see section 1.5.1).

Below may be found in serial and OpenMP versions.
“zacaz”, “zbcbz”, “zzszz”
These are generated if MISC=(USESCRATCH=TRUE).
The calculated orbital coefficients (zacaz, zbcbz for alpha and beta spins)
and overlap matrix (zzszz) of localization regions are contained in binary
format.
The files are temporary and thus not important to examine the results.

80

Contents

3.3.2 Output files in frequency calculation“restart_hess”
This is generated if FREQ=(PRINTHESS=TRUE).
The calculated Hessian matrix elements are contained in binary format.
The file may be used for setting Hessian matrix in
RESTARTHESS=BINARY in FREQ/FREQUENCY keyword (see section
1.5.1).

“hess.dat”
Thisisgenerated
if
FREQ=(INTERRUPTCP=TRUE)
orFREQ=(INTERRUPTSEMINUM=TRUE)
inFREQTYPE=1
or
FREQTYPE=2 in FREQ/FREQUENCY keyword.
The indices and calculated Hessian matrix elements are printed (seesection 1.5.2 for the format).
The file may be used for setting Hessian matrix in RESTARTHESS=ASCIIin FREQ/FREQUENCY keyword (see section 1.5.2).

“dftb_vib.out”
This is generated if FREQ=TRUE.
The information of harmonic vibrational frequency analysis is containedwith Gaussian like format.

81

Contents

3.3.3 Output files in MD simulation“restart”
This is generated if MD=TRUE.
The information of last MD step is contained in binary format.
The file may be used for restarting purpose in RESTART=TRUE in MD
keyword (see section 1.5.1).

“restart_chk”
This is generated if MD=TRUE.
The information of last checkpoint step is contained in binary format.
The file may be used for restarting purpose in RESTART=TRUE in MD
keyword after renaming to “restart” (see section 1.5.1).

“traject”
This is generated if MD=TRUE.
Cartesian coordinates during MD simulation are printed in xyz format.

“velocity”
This is generated if MD=TRUE.
Velocity components during MD simulation are printed in Å/fs unit in xyz
format.

“mulliken”
This is generated if MD=TRUE.
Mulliken net atomic charges during MD simulation are printed with thefollowing format.
 Serial and OpenMP versions
 Line 1: Number of atoms (Natom)
 Line 2: Comments (simulation time and MD step) Line 3–(Natom + 2): Atom label (1st column), atom symbol (2ndcolumn), and Mulliken net atomic charges (3rd column)
 MPI/OpenMP version
 Line 1: Number of shells (Nshell, 1st column) and Natom (2ndcolumn) Line 2: Comments (simulation time and MD step) Line 3–(Nshell + 2): Atom label (1st column), atom symbol (2ndcolumn), subshell symbol (3rd column), and Mulliken net charges
for the shell (4th column)
82

Contents

3.3.4 Output files in metadynamics
Additional output files from MD simulation (see section 3.3.3) are listed.
“restart_meta”
Thisisgenerated
if
MD=(METAD=TRUE),
MD=(METADYNAMICS=TRUE),
MD=(WTMETAD=TRUE),
or
MD=(WTMETADYNAMICS=TRUE).
The history of Gaussian bias potentials deposited during metadynamics
is contained in binary format.
The file may be used for restarting purpose in METARESTART=TRUE in
MD keyword (see section 1.5.1).

“restart_meta_chk”
Thisisgenerated
if
MD=(METAD=TRUE),
MD=(METADYNAMICS=TRUE),
MD=(WTMETAD=TRUE),
or
MD=(WTMETADYNAMICS=TRUE).
The history of Gaussian bias potentials deposited during metadynamics
at the last checkpoint step is contained in binary format.
The file may be used for restarting purpose in METARESTART=TRUE in
MD keyword after renaming to “restart_meta” (see section 1.5.1).

“biaspot”
Thisisgenerated
if
MD=(METAD=TRUE),
MD=(METADYNAMICS=TRUE),
MD=(WTMETAD=TRUE),
or
MD=(WTMETADYNAMICS=TRUE).
The history of Gaussian bias potentials deposited during metadynamics
are printed (see section 1.5.2 for the format).
The file may be used for restarting purpose in METAREADBIAS=TRUE
in MD keyword after renaming to “bias.dat” (see section 1.5.2).

“fes.dat”
This is generated if MD=(METAPRINTFES=TRUE).
The reconstructed free energy surface data are contained with thefollowing format.
 1st–(Number of CVs) columns: Grid point of CVs (unit is identical tothat of Gaussian width) Last column: Free energy in Hartree unit
83

Contents

References[1] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G.
Seifert, Phys. Rev. B 58, 7260 (1998).
[2] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[3] V. Eyert, J. Comput. Phys. 124, 271 (1996).
[4] P. Pulay, J. Comput. Chem. 3, 556 (1982).
[5] J. J. P. Stewart, P. Császár, P. Pulay, J. Comput. Chem. 3, 227 (1982).
[6] Y. Yang, H. Yu, D. York, Q. Cui, M. Elstner, J. Phys. Chem. A 111, 10861 (2007).
[7] M. Gaus, Q. Cui, M. Elstner, J. Chem. Theory Comput. 7, 931 (2011).
[8] M. Elstner, J. Phys. Chem. A 111, 5614 (2007).
[9] Y. Yang, H. Yu, D. York, M. Elstner, Q. Cui, J. Chem. Theory Comput. 4, 2067 (2008).
[10] Y. Nishimura, H. Nakai, J. Comput. Chem. 39, 105 (2018).
[11] C. Köhler, G. Seifert, U. Gerstmann, M. Elstner, H. Overhof, T. Frauenheim, Phys. Chem.
Chem. Phys. 3, 5109 (2001).
[12] C. Köhler, G. Seifert, T. Frauenheim, Chem. Phys. 309, 23 (2005).
[13] B. Hourahine, S. Sanna, B. Aradi, C. Köhler, T. Niehaus, T. Frauenheim, J. Phys. Chem. A
111, 5671 (2007).
[14] A. Domínguez, T. A. Niehaus, T. Frauenheim, J. Phys. Chem. A 119, 3535 (2015).
[15] I. Mayer, Chem. Phys. Lett. 97, 270 (1983).
[16] J. A. Kalinowski, B. Lesyng, J. D. Thompson, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. A
108, 2545 (2004).
[17] S. Kaminski, M. Gaus, M. Elstner, J. Phys. Chem. A 116, 11927 (2012).
[18] M. Gaus, H. Jin, D. Demapan, A. S. Christensen, P. Goyal, M. Elstner, Q. Cui, J. Chem. Theory
Comput. 11, 4205 (2015).
[19] B. Quentrec, C. Brot, J. Comput. Phys. 13, 430 (1973).
[20] W. Mattson, B. M. Rice, Comput. Phys. Commun. 119, 135 (1999).
[21] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B 51, 12947 (1995).
[22] G. Seifert, D. Porezag, T. Frauenheim, Int. J. Quantum Chem. 58, 185 (1996).
[23] W. Yang, T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).
[24] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry andPhysics, R. Zalesny, M. G. Papadopoulos, P. G. Mezey, J. Leszczynski, Eds.; Springer:
Dordrecht, pp. 97–127 (2011).
[25] 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).

84

Contents[26] H. Nishizawa, Y. Nishimura, M. Kobayashi, S. Irle, H. Nakai, J. Comput. Chem. 37, 1983(2016).
[27] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
[28] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 114, 5149 (2001).
[29] L. Zhechkov, T. Heine, S. Patchkovskii, G. Seifert, H. A. Duarte, J. Chem. Theory Comput. 1,
841 (2005).
[30] A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, W. M. Skiff, J. Am. Chem. Soc.
114, 10024 (1992).
[31] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 132, 154104 (2010).
[32] S. Grimme, S. Ehrlich, L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
[33] J. Řezáč, P. Hobza, J. Chem. Theory Comput. 8, 141 (2012).
[34] V. M. Miriyala, J. Řezáč, J. Comput. Chem. 38, 688 (2017).
[35] H. Kim, J.-M. Choi, W. A. Goddard III, J. Phys. Chem. Lett. 3, 360 (2012).
[36] R. Petraglia, S. N. Steinmann, C. Corminboeuf, Int. J. Quantum Chem. 115, 1265 (2015).
[37] J. Řezáč, J. Chem. Theory Comput. 13, 4804 (2017).
[38] M. Kubillus, T. Kubař, M. Gaus, J. Řezáč, M. Elstner, J. Chem. Theory Comput. 11, 332 (2015).
[39] C. G. Lambert, T. A. Darden, J. A. Board Jr., J. Comput. Phys. 126, 274 (1996).
[40] H. Jónsson, G. Mills, K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed
Phase Simulations, B. J. Berne, G. Ciccotti, D. F. Coker, Eds.; World Scientific: Singapore,pp. 385–404 (1998).
[41] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Phys. Rev. Lett. 97, 170201(2006).
[42] T. Atsumi, H. Nakai, Chem. Phys. Lett. 490, 102 (2010).
[43] Y. Nishimura, H. Nakai, J. Comput. Chem. 40, 1538 (2019).
[44] W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, J. Chem. Phys. 76, 637 (1982).
[45] L. V. Woodcock, Chem. Phys. Lett. 10, 257 (1971).
[46] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, M. L. Klein, Mol. Phys. 87, 1117 (1996).
[47] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem.
Phys. 81, 3684 (1984).
[48] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).
[49] G. Marsaglia, A. Zaman, W. W. Tsang, Stat. Probab. Lett. 8, 35 (1990).
[50] H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
[51] T. Atsumi, H. Nakai, J. Chem. Phys. 128, 094101 (2008).
[52] A. Laio, M. Parrinello, Proc. Natl. Acad. Sci. USA 99, 12562 (2002).
[53] A. Barducci, G. Bussi, M. Parrinello, Phys. Rev. Lett. 100, 020603 (2008).
[54] http://www.dftbplus.org/fileadmin/DFTBPLUS/public/dftbplus/latest/manual.pdf

85