diff --git a/README.md b/README.md index 828223e5d..1e86a7196 100644 --- a/README.md +++ b/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. diff --git a/dpdata/dftbplus/__init__.py b/dpdata/dftbplus/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/dpdata/dftbplus/output.py b/dpdata/dftbplus/output.py new file mode 100644 index 000000000..ba8f6c840 --- /dev/null +++ b/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 diff --git a/dpdata/plugins/dftbplus.py b/dpdata/plugins/dftbplus.py index e69de29bb..5c8b46828 100644 --- a/dpdata/plugins/dftbplus.py +++ b/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, + } diff --git a/tests/dftbplus/detailed.out b/tests/dftbplus/detailed.out new file mode 100644 index 000000000..21fe93fc5 --- /dev/null +++ b/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 + + diff --git a/tests/dftbplus/dftb_pin.hsd b/tests/dftbplus/dftb_pin.hsd new file mode 100644 index 000000000..02e0ed228 --- /dev/null +++ b/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 = {} diff --git a/tests/test_dftbplus.py b/tests/test_dftbplus.py new file mode 100644 index 000000000..2a2913a52 --- /dev/null +++ b/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()