Skip to content

Commit 452a0e5

Browse files
ruio248pre-commit-ci[bot]liaoruihao
authored
new feature DFTB+ (#511)
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: liaoruihao <liaoruihao@liaoruihaodeMacBook-Pro.local>
1 parent 5b57ac4 commit 452a0e5

7 files changed

Lines changed: 376 additions & 1 deletion

File tree

README.md

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,11 +88,12 @@ The `System` or `LabeledSystem` can be constructed from the following file forma
8888
| Amber | multi | True | True | LabeledSystem | 'amber/md' |
8989
| Amber/sqm | sqm.out | False | False | System | 'sqm/out' |
9090
| Gromacs | gro | True | False | System | 'gromacs/gro' |
91-
| ABACUS | STRU | False | False | System | 'abacus/stru' |
91+
| ABACUS | STRU | False | False | System | 'abacus/stru' |
9292
| ABACUS | STRU | False | True | LabeledSystem | 'abacus/scf' |
9393
| ABACUS | cif | True | True | LabeledSystem | 'abacus/md' |
9494
| ABACUS | STRU | True | True | LabeledSystem | 'abacus/relax' |
9595
| ase | structure | True | True | MultiSystems | 'ase/structure' |
96+
| DFTB+ | dftbplus | False | True | LabeledSystem | 'dftbplus' |
9697

9798

9899
The Class `dpdata.MultiSystems` can read data from a dir which may contains many files of different systems, or from single xyz file which contains different systems.

dpdata/dftbplus/__init__.py

Whitespace-only changes.

dpdata/dftbplus/output.py

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
from typing import Tuple
2+
3+
import numpy as np
4+
5+
6+
def read_dftb_plus(fn_1: str, fn_2: str) -> Tuple[str, np.ndarray, float, np.ndarray]:
7+
"""Read from DFTB+ input and output.
8+
9+
Parameters
10+
----------
11+
fn_1 : str
12+
DFTB+ input file name
13+
fn_2 : str
14+
DFTB+ output file name
15+
16+
Returns
17+
-------
18+
str
19+
atomic symbols
20+
np.ndarray
21+
atomic coordinates
22+
float
23+
total potential energy
24+
np.ndarray
25+
atomic forces
26+
27+
"""
28+
coord = None
29+
symbols = None
30+
forces = None
31+
energy = None
32+
with open(fn_1) as f:
33+
flag = 0
34+
for line in f:
35+
if flag == 1:
36+
flag += 1
37+
elif flag == 2:
38+
components = line.split()
39+
flag += 1
40+
elif line.startswith("Geometry"):
41+
flag = 1
42+
coord = []
43+
symbols = []
44+
elif flag in (3, 4, 5, 6):
45+
s = line.split()
46+
components_num = int(s[1])
47+
symbols.append(components[components_num - 1])
48+
coord.append([float(s[2]), float(s[3]), float(s[4])])
49+
flag += 1
50+
if flag == 7:
51+
flag = 0
52+
with open(fn_2) as f:
53+
flag = 0
54+
for line in f:
55+
if line.startswith("Total Forces"):
56+
flag = 8
57+
forces = []
58+
elif flag in (8, 9, 10, 11):
59+
s = line.split()
60+
forces.append([float(s[1]), float(s[2]), float(s[3])])
61+
flag += 1
62+
if flag == 12:
63+
flag = 0
64+
elif line.startswith("Total energy:"):
65+
s = line.split()
66+
energy = float(s[2])
67+
flag = 0
68+
69+
symbols = np.array(symbols)
70+
forces = np.array(forces)
71+
coord = np.array(coord)
72+
assert coord.shape == forces.shape
73+
74+
return symbols, coord, energy, forces

dpdata/plugins/dftbplus.py

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
import numpy as np
2+
3+
from dpdata.dftbplus.output import read_dftb_plus
4+
from dpdata.format import Format
5+
from dpdata.unit import EnergyConversion, ForceConversion
6+
7+
energy_convert = EnergyConversion("hartree", "eV").value()
8+
force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
9+
10+
11+
@Format.register("dftbplus")
12+
class DFTBplusFormat(Format):
13+
"""The DFTBplusFormat class handles files in the DFTB+ format.
14+
15+
This class provides a method to read DFTB+ files from a labeled system and
16+
returns a dictionary containing various properties of the system.For more
17+
information, please refer to the official documentation at the following URL:
18+
https://dftbplus.org/documentation
19+
20+
Attributes
21+
----------
22+
None
23+
24+
Methods
25+
-------
26+
from_labeled_system(file_paths, **kwargs): Reads system information from files.
27+
28+
"""
29+
30+
def from_labeled_system(self, file_paths, **kwargs):
31+
"""Reads system information from the given DFTB+ file paths.
32+
33+
Parameters
34+
----------
35+
file_paths : tuple
36+
A tuple containing the input and output file paths.
37+
- Input file (file_in): Contains information about symbols and coord.
38+
- Output file (file_out): Contains information about energy and force.
39+
**kwargs : dict
40+
other parameters
41+
42+
"""
43+
file_in, file_out = file_paths
44+
symbols, coord, energy, forces = read_dftb_plus(file_in, file_out)
45+
last_occurrence = {v: i for i, v in enumerate(symbols)}
46+
atom_names = np.array(sorted(np.unique(symbols), key=last_occurrence.get))
47+
atom_types = np.array([np.where(atom_names == s)[0][0] for s in symbols])
48+
atom_numbs = np.array([np.sum(atom_types == i) for i in range(len(atom_names))])
49+
natoms = coord.shape[0]
50+
51+
return {
52+
"atom_types": atom_types,
53+
"atom_names": list(atom_names),
54+
"atom_numbs": list(atom_numbs),
55+
"coords": coord.reshape((1, natoms, 3)),
56+
"energies": np.array([energy * energy_convert]),
57+
"forces": (forces * force_convert).reshape((1, natoms, 3)),
58+
"cells": np.zeros((1, 3, 3)),
59+
"orig": np.zeros(3),
60+
"nopbc": True,
61+
}

tests/dftbplus/detailed.out

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
Fermi distribution function
2+
3+
Calculation with static geometry
4+
5+
6+
********************************************************************************
7+
iSCC Total electronic Diff electronic SCC error
8+
88 -0.33554958E+01 -0.31884693E-07 0.11352275E-06
9+
********************************************************************************
10+
11+
Total charge: -0.00000000
12+
13+
Atomic gross charges (e)
14+
Atom Charge
15+
1 -0.48336965
16+
2 0.03575828
17+
3 0.22380553
18+
4 0.22380583
19+
20+
Nr. of electrons (up): 8.00000000
21+
Atom populations (up)
22+
Atom Population
23+
1 5.48336965
24+
2 0.96424172
25+
3 0.77619447
26+
4 0.77619417
27+
28+
l-shell populations (up)
29+
Atom Sh. l Population
30+
1 1 0 1.65733979
31+
1 2 1 3.82602986
32+
2 1 0 0.96424172
33+
2 2 1 0.00000000
34+
3 1 0 0.77619447
35+
3 2 1 0.00000000
36+
4 1 0 0.77619417
37+
4 2 1 0.00000000
38+
39+
Orbital populations (up)
40+
Atom Sh. l m Population Label
41+
1 1 0 0 1.65733979 s
42+
1 2 1 -1 1.21678996 p_y
43+
1 2 1 0 1.46556336 p_z
44+
1 2 1 1 1.14367654 p_x
45+
2 1 0 0 0.96424172 s
46+
2 2 1 -1 0.00000000 p_y
47+
2 2 1 0 0.00000000 p_z
48+
2 2 1 1 0.00000000 p_x
49+
3 1 0 0 0.77619447 s
50+
3 2 1 -1 0.00000000 p_y
51+
3 2 1 0 0.00000000 p_z
52+
3 2 1 1 0.00000000 p_x
53+
4 1 0 0 0.77619417 s
54+
4 2 1 -1 0.00000000 p_y
55+
4 2 1 0 0.00000000 p_z
56+
4 2 1 1 0.00000000 p_x
57+
58+
Fermi level: -0.2370222488 H -6.4497 eV
59+
Band energy: -3.2187242078 H -87.5859 eV
60+
TS: 0.0000000000 H 0.0000 eV
61+
Band free energy (E-TS): -3.2187242078 H -87.5859 eV
62+
Extrapolated E(0K): -3.2187242078 H -87.5859 eV
63+
Input / Output electrons (q): 8.0000000000 8.0000000000
64+
65+
Energy H0: -3.3599884076 H -91.4299 eV
66+
Energy SCC: 0.0056352830 H 0.1533 eV
67+
Energy 3rd: -0.0011426808 H -0.0311 eV
68+
Total Electronic energy: -3.3554958054 H -91.3077 eV
69+
Repulsive energy: 0.0590974170 H 1.6081 eV
70+
Total energy: -3.2963983884 H -89.6996 eV
71+
Extrapolated to 0: -3.2963983884 H -89.6996 eV
72+
Total Mermin free energy: -3.2963983884 H -89.6996 eV
73+
Force related energy: -3.2963983884 H -89.6996 eV
74+
75+
SCC converged
76+
77+
Total Forces
78+
1 0.016567056203 0.002817951422 0.005634574270
79+
2 -0.018803818530 -0.000002880649 -0.000006015442
80+
3 0.001118562874 -0.005291070259 -0.000870711110
81+
4 0.001118199454 0.002475999486 -0.004757847718
82+
83+
Maximal derivative component: 0.188038E-01 au
84+
85+
Dipole moment: -0.06792979 -0.20495079 -0.40960550 au
86+
Dipole moment: -0.17266032 -0.52093298 -1.04111341 Debye
87+
88+

tests/dftbplus/dftb_pin.hsd

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
Geometry = GenFormat {
2+
4 C
3+
N H
4+
1 1 1.014150 0.112320 0.047370
5+
2 2 3.909390 0.037985 -0.101159
6+
3 2 0.702550 -0.851820 -0.060860
7+
4 2 0.702550 0.603740 -0.789160
8+
}
9+
Driver = {}
10+
Hamiltonian = DFTB {
11+
SCC = Yes
12+
MaxAngularMomentum = {
13+
N = "p"
14+
H = "p"
15+
}
16+
SlaterKosterFiles = {
17+
N-N = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/N-N.skf"
18+
H-H = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/H-H.skf"
19+
H-N = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/H-N.skf"
20+
N-H = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/N-H.skf"
21+
}
22+
ThirdOrderFull = Yes
23+
HubbardDerivs = {
24+
H = -0.1857
25+
N = -0.1535
26+
}
27+
HCorrection = Damping {
28+
Exponent = 4.0
29+
}
30+
PolynomialRepulsive = {}
31+
ShellResolvedSCC = No
32+
OldSKInterpolation = No
33+
RangeSeparated = None {}
34+
ReadInitialCharges = No
35+
InitialCharges = {}
36+
SCCTolerance = 1.0000000000000001E-005
37+
ConvergentSCCOnly = Yes
38+
SpinPolarisation = {}
39+
ElectricField = {}
40+
Solver = RelativelyRobust {}
41+
Charge = 0.0000000000000000
42+
MaxSCCIterations = 100
43+
OnSiteCorrection = {}
44+
Dispersion = {}
45+
Solvation = {}
46+
Electrostatics = GammaFunctional {}
47+
ThirdOrder = No
48+
Differentiation = FiniteDiff {
49+
Delta = 1.2207031250000000E-004
50+
}
51+
ForceEvaluation = "Traditional"
52+
Mixer = Broyden {
53+
MixingParameter = 0.20000000000000001
54+
InverseJacobiWeight = 1.0000000000000000E-002
55+
MinimalWeight = 1.0000000000000000
56+
MaximalWeight = 100000.00000000000
57+
WeightFactor = 1.0000000000000000E-002
58+
}
59+
Filling = Fermi {
60+
Temperature = 0.0000000000000000
61+
}
62+
}
63+
Options = {
64+
WriteDetailedOut = Yes
65+
WriteAutotestTag = No
66+
WriteDetailedXML = No
67+
WriteResultsTag = No
68+
RestartFrequency = 20
69+
RandomSeed = 0
70+
WriteHS = No
71+
WriteRealHS = No
72+
MinimiseMemoryUsage = No
73+
ShowFoldedCoords = No
74+
TimingVerbosity = 1
75+
WriteChargesAsText = No
76+
}
77+
Analysis = {
78+
CalculateForces = Yes
79+
ProjectStates = {}
80+
WriteEigenvectors = No
81+
WriteBandOut = Yes
82+
MullikenAnalysis = Yes
83+
WriteNetCharges = No
84+
AtomResolvedEnergies = No
85+
}
86+
ParserOptions = {
87+
ParserVersion = 11
88+
WriteHSDInput = Yes
89+
StopAfterParsing = No
90+
IgnoreUnprocessedNodes = No
91+
}
92+
Reks = None {}
93+
ExcitedState = {}

tests/test_dftbplus.py

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
import unittest
2+
3+
import numpy as np
4+
from comp_sys import CompLabeledSys, IsNoPBC
5+
from context import dpdata
6+
7+
8+
class TestDeepmdLoadAmmonia(unittest.TestCase, CompLabeledSys, IsNoPBC):
9+
def setUp(self):
10+
energy_convert = dpdata.unit.EnergyConversion("hartree", "eV").value()
11+
force_convert = dpdata.unit.ForceConversion(
12+
"hartree/bohr", "eV/angstrom"
13+
).value()
14+
15+
self.system_1 = dpdata.LabeledSystem(
16+
("dftbplus/dftb_pin.hsd", "dftbplus/detailed.out"), fmt="dftbplus"
17+
)
18+
19+
self.system_2 = dpdata.LabeledSystem(
20+
data={
21+
"atom_types": np.array([0, 1, 1, 1]),
22+
"atom_names": ["N", "H"],
23+
"atom_numbs": [1, 3],
24+
"coords": np.array(
25+
[
26+
[
27+
[1.014150, 0.112320, 0.047370],
28+
[3.909390, 0.037985, -0.101159],
29+
[0.702550, -0.851820, -0.060860],
30+
[0.702550, 0.603740, -0.789160],
31+
]
32+
]
33+
),
34+
"energies": np.array([-3.2963983884]) * energy_convert,
35+
"forces": np.array(
36+
[
37+
[
38+
[0.016567056203, 0.002817951422, 0.005634574270],
39+
[-0.018803818530, -0.000002880649, -0.000006015442],
40+
[0.001118562874, -0.005291070259, -0.000870711110],
41+
[0.001118199454, 0.002475999486, -0.004757847718],
42+
]
43+
]
44+
)
45+
* force_convert,
46+
"cells": np.zeros((1, 3, 3)),
47+
"orig": np.zeros(3),
48+
"nopbc": True,
49+
}
50+
)
51+
self.places = 6
52+
self.e_places = 6
53+
self.f_places = 6
54+
self.v_places = 6
55+
56+
57+
if __name__ == "__main__":
58+
unittest.main()

0 commit comments

Comments
 (0)