new feature DFTB+

README.md

@@ -88,11 +88,12 @@ The `System` or `LabeledSystem` can be constructed from the following file forma
 | Amber   | multi       | True         | True    | LabeledSystem | 'amber/md'           |
 | Amber/sqm | sqm.out   | False        | False   | System        | 'sqm/out'            |
 | Gromacs | gro         | True         | False   | System        | 'gromacs/gro'        |
-| ABACUS  | STRU        | False        | False   | System        | 'abacus/stru'         |
+| ABACUS  | STRU        | False        | False   | System        | 'abacus/stru'        |
 | ABACUS  | STRU        | False        | True    | LabeledSystem | 'abacus/scf'         |
 | ABACUS  | cif         | True         | True    | LabeledSystem | 'abacus/md'          |
 | ABACUS  | STRU        | True         | True    | LabeledSystem | 'abacus/relax'       |
 | ase     | structure   | True         | True    | MultiSystems  | 'ase/structure'      |
+| DFTB+   | dftbplus    | False        | True    | LabeledSystem | 'dftbplus'           |
 
 
 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/output.py

@@ -0,0 +1,74 @@
+from typing import Tuple
+
+import numpy as np
+
+
+def read_dftb_plus(fn_1: str, fn_2: str) -> Tuple[str, np.ndarray, float, np.ndarray]:
+    """Read from DFTB+ input and output.
+
+    Parameters
+    ----------
+    fn_1 : str
+        DFTB+ input file name
+    fn_2 : str
+        DFTB+ output file name
+
+    Returns
+    -------
+    str
+        atomic symbols
+    np.ndarray
+        atomic coordinates
+    float
+        total potential energy
+    np.ndarray
+        atomic forces
+
+    """
+    coord = None
+    symbols = None
+    forces = None
+    energy = None
+    with open(fn_1) as f:
+        flag = 0
+        for line in f:
+            if flag == 1:
+                flag += 1
+            elif flag == 2:
+                components = line.split()
+                flag += 1
+            elif line.startswith("Geometry"):
+                flag = 1
+                coord = []
+                symbols = []
+            elif flag in (3, 4, 5, 6):
+                s = line.split()
+                components_num = int(s[1])
+                symbols.append(components[components_num - 1])
+                coord.append([float(s[2]), float(s[3]), float(s[4])])
+                flag += 1
+                if flag == 7:
+                    flag = 0
+    with open(fn_2) as f:
+        flag = 0
+        for line in f:
+            if line.startswith("Total Forces"):
+                flag = 8
+                forces = []
+            elif flag in (8, 9, 10, 11):
+                s = line.split()
+                forces.append([float(s[1]), float(s[2]), float(s[3])])
+                flag += 1
+                if flag == 12:
+                    flag = 0
+            elif line.startswith("Total energy:"):
+                s = line.split()
+                energy = float(s[2])
+                flag = 0
+
+    symbols = np.array(symbols)
+    forces = np.array(forces)
+    coord = np.array(coord)
+    assert coord.shape == forces.shape
+
+    return symbols, coord, energy, forces

dpdata/plugins/dftbplus.py

@@ -0,0 +1,61 @@
+import numpy as np
+
+from dpdata.dftbplus.output import read_dftb_plus
+from dpdata.format import Format
+from dpdata.unit import EnergyConversion, ForceConversion
+
+energy_convert = EnergyConversion("hartree", "eV").value()
+force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
+
+
+@Format.register("dftbplus")
+class DFTBplusFormat(Format):
+    """The DFTBplusFormat class handles files in the DFTB+ format.
+
+    This class provides a method to read DFTB+ files from a labeled system and
+    returns a dictionary containing various properties of the system.For more
+    information, please refer to the official documentation at the following URL:
+    https://dftbplus.org/documentation
+
+    Attributes
+    ----------
+        None
+
+    Methods
+    -------
+        from_labeled_system(file_paths, **kwargs): Reads system information from files.
+
+    """
+
+    def from_labeled_system(self, file_paths, **kwargs):
+        """Reads system information from the given DFTB+ file paths.
+
+        Parameters
+        ----------
+        file_paths : tuple
+            A tuple containing the input and output file paths.
+            - Input file (file_in): Contains information about symbols and coord.
+            - Output file (file_out): Contains information about energy and force.
+        **kwargs : dict
+            other parameters
+
+        """
+        file_in, file_out = file_paths
+        symbols, coord, energy, forces = read_dftb_plus(file_in, file_out)
+        last_occurrence = {v: i for i, v in enumerate(symbols)}
+        atom_names = np.array(sorted(np.unique(symbols), key=last_occurrence.get))
+        atom_types = np.array([np.where(atom_names == s)[0][0] for s in symbols])
+        atom_numbs = np.array([np.sum(atom_types == i) for i in range(len(atom_names))])
+        natoms = coord.shape[0]
+
+        return {
+            "atom_types": atom_types,
+            "atom_names": list(atom_names),
+            "atom_numbs": list(atom_numbs),
+            "coords": coord.reshape((1, natoms, 3)),
+            "energies": np.array([energy * energy_convert]),
+            "forces": (forces * force_convert).reshape((1, natoms, 3)),
+            "cells": np.zeros((1, 3, 3)),
+            "orig": np.zeros(3),
+            "nopbc": True,
+        }

tests/dftbplus/detailed.out

@@ -0,0 +1,88 @@
+ Fermi distribution function
+
+Calculation with static geometry
+
+
+********************************************************************************
+ iSCC Total electronic   Diff electronic      SCC error    
+   88   -0.33554958E+01   -0.31884693E-07    0.11352275E-06
+********************************************************************************
+
+ Total charge:    -0.00000000
+
+ Atomic gross charges (e)
+ Atom           Charge
+    1      -0.48336965
+    2       0.03575828
+    3       0.22380553
+    4       0.22380583
+
+Nr. of electrons (up):      8.00000000
+Atom populations (up)
+ Atom       Population
+    1       5.48336965
+    2       0.96424172
+    3       0.77619447
+    4       0.77619417
+
+l-shell populations (up)
+ Atom Sh.   l       Population
+    1   1   0       1.65733979
+    1   2   1       3.82602986
+    2   1   0       0.96424172
+    2   2   1       0.00000000
+    3   1   0       0.77619447
+    3   2   1       0.00000000
+    4   1   0       0.77619417
+    4   2   1       0.00000000
+
+Orbital populations (up)
+ Atom Sh.   l   m       Population  Label
+    1   1   0   0       1.65733979  s
+    1   2   1  -1       1.21678996  p_y
+    1   2   1   0       1.46556336  p_z
+    1   2   1   1       1.14367654  p_x
+    2   1   0   0       0.96424172  s
+    2   2   1  -1       0.00000000  p_y
+    2   2   1   0       0.00000000  p_z
+    2   2   1   1       0.00000000  p_x
+    3   1   0   0       0.77619447  s
+    3   2   1  -1       0.00000000  p_y
+    3   2   1   0       0.00000000  p_z
+    3   2   1   1       0.00000000  p_x
+    4   1   0   0       0.77619417  s
+    4   2   1  -1       0.00000000  p_y
+    4   2   1   0       0.00000000  p_z
+    4   2   1   1       0.00000000  p_x
+
+Fermi level:                        -0.2370222488 H           -6.4497 eV
+Band energy:                        -3.2187242078 H          -87.5859 eV
+TS:                                  0.0000000000 H            0.0000 eV
+Band free energy (E-TS):            -3.2187242078 H          -87.5859 eV
+Extrapolated E(0K):                 -3.2187242078 H          -87.5859 eV
+Input / Output electrons (q):      8.0000000000      8.0000000000
+
+Energy H0:                          -3.3599884076 H          -91.4299 eV
+Energy SCC:                          0.0056352830 H            0.1533 eV
+Energy 3rd:                         -0.0011426808 H           -0.0311 eV
+Total Electronic energy:            -3.3554958054 H          -91.3077 eV
+Repulsive energy:                    0.0590974170 H            1.6081 eV
+Total energy:                       -3.2963983884 H          -89.6996 eV
+Extrapolated to 0:                  -3.2963983884 H          -89.6996 eV
+Total Mermin free energy:           -3.2963983884 H          -89.6996 eV
+Force related energy:               -3.2963983884 H          -89.6996 eV
+
+SCC converged
+
+Total Forces
+    1      0.016567056203      0.002817951422      0.005634574270
+    2     -0.018803818530     -0.000002880649     -0.000006015442
+    3      0.001118562874     -0.005291070259     -0.000870711110
+    4      0.001118199454      0.002475999486     -0.004757847718
+
+Maximal derivative component:        0.188038E-01 au
+
+Dipole moment:   -0.06792979   -0.20495079   -0.40960550 au
+Dipole moment:   -0.17266032   -0.52093298   -1.04111341 Debye
+
+

tests/dftbplus/dftb_pin.hsd

@@ -0,0 +1,93 @@
+Geometry = GenFormat {
+4 C
+N H
+1     1        1.014150        0.112320        0.047370
+2     2        3.909390        0.037985       -0.101159
+3     2        0.702550       -0.851820       -0.060860
+4     2        0.702550        0.603740       -0.789160
+}
+Driver = {}
+Hamiltonian = DFTB {
+  SCC = Yes
+  MaxAngularMomentum = {
+    N = "p"
+    H = "p"
+  }
+  SlaterKosterFiles = {
+    N-N = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/N-N.skf"
+    H-H = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/H-H.skf"
+    H-N = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/H-N.skf"
+    N-H = "/home/jz748/devel/git/Programs/Amber18/dat/slko/3ob-3-1/N-H.skf"
+  }
+  ThirdOrderFull = Yes
+  HubbardDerivs = {
+    H = -0.1857
+    N = -0.1535
+  }
+  HCorrection = Damping {
+    Exponent = 4.0
+  }
+  PolynomialRepulsive = {}
+  ShellResolvedSCC = No
+  OldSKInterpolation = No
+  RangeSeparated = None {}
+  ReadInitialCharges = No
+  InitialCharges = {}
+  SCCTolerance = 1.0000000000000001E-005
+  ConvergentSCCOnly = Yes
+  SpinPolarisation = {}
+  ElectricField = {}
+  Solver = RelativelyRobust {}
+  Charge = 0.0000000000000000
+  MaxSCCIterations = 100
+  OnSiteCorrection = {}
+  Dispersion = {}
+  Solvation = {}
+  Electrostatics = GammaFunctional {}
+  ThirdOrder = No
+  Differentiation = FiniteDiff {
+    Delta = 1.2207031250000000E-004
+  }
+  ForceEvaluation = "Traditional"
+  Mixer = Broyden {
+    MixingParameter = 0.20000000000000001
+    InverseJacobiWeight = 1.0000000000000000E-002
+    MinimalWeight = 1.0000000000000000
+    MaximalWeight = 100000.00000000000
+    WeightFactor = 1.0000000000000000E-002
+  }
+  Filling = Fermi {
+    Temperature = 0.0000000000000000
+  }
+}
+Options = {
+  WriteDetailedOut = Yes
+  WriteAutotestTag = No
+  WriteDetailedXML = No
+  WriteResultsTag = No
+  RestartFrequency = 20
+  RandomSeed = 0
+  WriteHS = No
+  WriteRealHS = No
+  MinimiseMemoryUsage = No
+  ShowFoldedCoords = No
+  TimingVerbosity = 1
+  WriteChargesAsText = No
+}
+Analysis = {
+  CalculateForces = Yes
+  ProjectStates = {}
+  WriteEigenvectors = No
+  WriteBandOut = Yes
+  MullikenAnalysis = Yes
+  WriteNetCharges = No
+  AtomResolvedEnergies = No
+}
+ParserOptions = {
+  ParserVersion = 11
+  WriteHSDInput = Yes
+  StopAfterParsing = No
+  IgnoreUnprocessedNodes = No
+}
+Reks = None {}
+ExcitedState = {}

tests/test_dftbplus.py

@@ -0,0 +1,58 @@
+import unittest
+
+import numpy as np
+from comp_sys import CompLabeledSys, IsNoPBC
+from context import dpdata
+
+
+class TestDeepmdLoadAmmonia(unittest.TestCase, CompLabeledSys, IsNoPBC):
+    def setUp(self):
+        energy_convert = dpdata.unit.EnergyConversion("hartree", "eV").value()
+        force_convert = dpdata.unit.ForceConversion(
+            "hartree/bohr", "eV/angstrom"
+        ).value()
+
+        self.system_1 = dpdata.LabeledSystem(
+            ("dftbplus/dftb_pin.hsd", "dftbplus/detailed.out"), fmt="dftbplus"
+        )
+
+        self.system_2 = dpdata.LabeledSystem(
+            data={
+                "atom_types": np.array([0, 1, 1, 1]),
+                "atom_names": ["N", "H"],
+                "atom_numbs": [1, 3],
+                "coords": np.array(
+                    [
+                        [
+                            [1.014150, 0.112320, 0.047370],
+                            [3.909390, 0.037985, -0.101159],
+                            [0.702550, -0.851820, -0.060860],
+                            [0.702550, 0.603740, -0.789160],
+                        ]
+                    ]
+                ),
+                "energies": np.array([-3.2963983884]) * energy_convert,
+                "forces": np.array(
+                    [
+                        [
+                            [0.016567056203, 0.002817951422, 0.005634574270],
+                            [-0.018803818530, -0.000002880649, -0.000006015442],
+                            [0.001118562874, -0.005291070259, -0.000870711110],
+                            [0.001118199454, 0.002475999486, -0.004757847718],
+                        ]
+                    ]
+                )
+                * force_convert,
+                "cells": np.zeros((1, 3, 3)),
+                "orig": np.zeros(3),
+                "nopbc": True,
+            }
+        )
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 6
+
+
+if __name__ == "__main__":
+    unittest.main()