Add unit.py to manage unit conversions

dpdata/abacus/scf.py

@@ -1,11 +1,10 @@
 import os,sys
 import numpy as np
+from ..unit import EnergyConversion, PressureConversion, LengthConversion
 
-bohr2ang = 0.5291770
-ry2ev = 13.605698
-kbar2evperang3 = 1e3 / 1.6021892e6
-# The consts are cited from $ABACUS_ROOT/source/src_global/constant.h
-
+bohr2ang = LengthConversion("bohr", "angstrom").value()
+ry2ev = EnergyConversion("rydberg", "eV").value()
+kbar2evperang3 = PressureConversion("kbar", "eV/angstrom^3").value()
 
 def get_block (lines, keyword, skip = 0, nlines = None):
     ret = []

dpdata/amber/md.py

@@ -3,9 +3,11 @@
 from scipy.io import netcdf
 import numpy as np
 from dpdata.amber.mask import pick_by_amber_mask
+from dpdata.unit import EnergyConversion
+from ..periodic_table import ELEMENTS
 
-kcalmol2eV= 0.04336410390059322
-symbols = ['X', 'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Nh', 'Fl', 'Mc', 'Lv', 'Ts', 'Og']
+kcalmol2eV = EnergyConversion("kcal_mol", "eV").value()
+symbols = ['X'] + ELEMENTS
 
 energy_convert = kcalmol2eV
 force_convert = energy_convert

dpdata/amber/sqm.py

@@ -1,6 +1,8 @@
 import numpy as np
-from ..periodic_table import ELEMENTS
+from dpdata.periodic_table import ELEMENTS
+from dpdata.unit import EnergyConversion
 
+kcal2ev = EnergyConversion("kcal_mol", "eV").value()
 
 START = 0
 READ_ENERGY = 1
@@ -47,7 +49,8 @@ def parse_sqm_out(fname):
                 if len(coords) == len(charges):
                     flag = START
             elif flag == READ_FORCES:
-                forces.append([float(x) for x in line.strip().split()[-3:]])
+                ll = line.strip()
+                forces.append([float(ll[-60:-40]), float(ll[-40:-20]), float(ll[-20:])])
                 if len(forces) == len(charges):
                     flag = START
 
@@ -62,7 +65,7 @@ def parse_sqm_out(fname):
     data['coords'] = np.array([coords])
 
     energies = np.array(energies)
-    forces = np.array([forces], dtype=np.float32)
+    forces = np.array([forces], dtype=np.float32) * kcal2ev
     if len(forces) > 0:
         data['energies'] = energies
         data['forces'] = forces
@@ -99,4 +102,4 @@ def make_sqm_in(data, fname=None, frame_idx=0, **kwargs):
     if fname is not None:
         with open(fname, 'w') as fp:
             fp.write(ret)
-    return ret
+    return ret

dpdata/cp2k/output.py

@@ -3,11 +3,12 @@
 import re
 from collections import OrderedDict
 from .cell import cell_to_low_triangle
+from ..unit import EnergyConversion, LengthConversion, ForceConversion, PressureConversion
 
 #%%
-AU_TO_ANG = 5.29177208590000E-01
-AU_TO_EV = 2.72113838565563E+01
-AU_TO_EV_EVERY_ANG = AU_TO_EV/AU_TO_ANG
+AU_TO_ANG = LengthConversion("bohr", "angstrom").value()
+AU_TO_EV = EnergyConversion("hartree", "eV").value()
+AU_TO_EV_EVERY_ANG = ForceConversion("hartree/bohr", "eV/angstrom").value()
 delimiter_patterns=[]
 delimiter_p1 = re.compile(r'^ \* GO CP2K GO! \*+')
 delimiter_p2 = re.compile(r'^ \*+')
@@ -221,9 +222,9 @@ def get_frames (fname) :
     coord_flag = False
     force_flag = False
     stress_flag = False
-    eV = 2.72113838565563E+01 # hatree to eV
-    angstrom = 5.29177208590000E-01 # Bohr to Angstrom
-    GPa = 160.21766208 # 1 eV/(Angstrom^3) = 160.21 GPa
+    eV = EnergyConversion("hartree", "eV").value()
+    angstrom = LengthConversion("bohr", "angstrom").value()
+    GPa = PressureConversion("eV/angstrom^3", "GPa").value()
     atom_symbol_list = []
     cell = []
     coord = []

dpdata/gaussian/log.py

@@ -1,14 +1,12 @@
 import numpy as np
+from ..unit import LengthConversion, EnergyConversion, ForceConversion
+from ..periodic_table import ELEMENTS
 
+length_convert = LengthConversion("bohr", "angstrom").value()
+energy_convert = EnergyConversion("hartree", "eV").value()
+force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
 
-hartree2ev = 27.211386018
-bohr2ang = 0.52917721067
-
-length_convert = bohr2ang
-energy_convert = hartree2ev
-force_convert = energy_convert / length_convert
-
-symbols = ['X', 'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Nh', 'Fl', 'Mc', 'Lv', 'Ts', 'Og']
+symbols = ["X"] + ELEMENTS
 
 def to_system_data(file_name, md=False):
     data = {}

dpdata/gromacs/gro.py

@@ -1,11 +1,10 @@
 #!/usr/bin/env python3
 import re
 import numpy as np
-import json
-import os
+from ..unit import LengthConversion
 
-nm2ang = 10.
-ang2nm = 1. / nm2ang
+nm2ang = LengthConversion("nm", "angstrom").value()
+ang2nm = LengthConversion("angstrom", "nm").value()
 cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7]
 
 def _format_atom_name(atom_name):

dpdata/pwmat/atomconfig.py

@@ -1,5 +1,5 @@
 #!/usr/bin/python3 
-
+from ..periodic_table import ELEMENTS
 import numpy as np
 
 def _to_system_data_lower(lines) :
@@ -39,13 +39,6 @@ def _to_system_data_lower(lines) :
         for jj in range(ii) :
             atom_types.append(idx)
     system['atom_types'] = np.array(atom_types, dtype = int)
-    ELEMENTS=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', \
-            'Sc', 'Ti', 'V', 'Cr','Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', \
-            'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',\
-            'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', \
-            'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', \
-            'Md', 'No', 'Lr']
-
     system['atom_names'] = [ELEMENTS[ii-1] for ii in np.unique(sorted(atomic_number))]
     return system
 
@@ -80,12 +73,6 @@ def from_system_data(system, f_idx = 0, skip_zeros = True) :
     for ii, jj in zip(atom_numbs, atom_names):
         for kk in range(ii):
             symbal.append(jj) 
-    ELEMENTS=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', \
-            'Sc', 'Ti', 'V', 'Cr','Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', \
-            'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',\
-            'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', \
-            'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', \
-            'Md', 'No', 'Lr']
     atomic_numbers = []
     for ii in symbal:
         atomic_numbers.append(ELEMENTS.index(ii)+1)

dpdata/pwmat/movement.py

@@ -1,4 +1,5 @@
 import numpy as np
+from ..periodic_table import ELEMENTS
 
 def system_info (lines, type_idx_zero = False) :
     atom_names = []
@@ -32,12 +33,6 @@ def system_info (lines, type_idx_zero = False) :
                 atom_types.append(idx)
             else :
                 atom_types.append(idx+1)
-    ELEMENTS=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', \
-            'Sc', 'Ti', 'V', 'Cr','Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', \
-            'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',\
-            'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', \
-            'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', \
-            'Md', 'No', 'Lr']
     for ii in np.unique(sorted(atomic_number)):
         atom_names.append(ELEMENTS[ii-1])
     return atom_names, atom_numbs, np.array(atom_types, dtype = int), nelm

dpdata/qe/traj.py

@@ -1,16 +1,14 @@
 #!/usr/bin/python3 
-
 import numpy as np
-import dpdata,warnings
+import dpdata, warnings
+from ..unit import EnergyConversion, LengthConversion, ForceConversion, PressureConversion
 
-ry2ev = 13.605693009
-hartree2ev = 27.211386018
-bohr2ang = 0.52917721067
-kbar2evperang3 = 1e3 / 1.602176621e6
+ry2ev = EnergyConversion("rydberg", "eV").value()
+kbar2evperang3 = PressureConversion("kbar", "eV/angstrom^3").value()
 
-length_convert = bohr2ang
-energy_convert = hartree2ev
-force_convert = energy_convert / length_convert
+length_convert = LengthConversion("bohr", "angstrom").value()
+energy_convert = EnergyConversion("hartree", "eV").value()
+force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
 
 def load_key (lines, key) :
     for ii in lines :

dpdata/unit.py

@@ -0,0 +1,169 @@
+from abc import ABC
+from scipy import constants
+
+AVOGADRO = constants.Avogadro                             # Avagadro constant
+ELE_CHG  = constants.elementary_charge                    # Elementary Charge, in C
+BOHR     = constants.value("atomic unit of length")       # Bohr, in m
+HARTREE  = constants.value("atomic unit of energy")       # Hartree, in Jole
+RYDBERG  = constants.Rydberg * constants.h * constants.c  # Rydberg, in Jole
+
+# energy conversions
+econvs = {
+    "eV": 1.0,
+    "hartree": HARTREE / ELE_CHG,
+    "kJ_mol": 1 / (ELE_CHG * AVOGADRO / 1000),
+    "kcal_mol": 1 / (ELE_CHG * AVOGADRO / 1000 / 4.184),
+    "rydberg": RYDBERG / ELE_CHG,
+    "J": 1 / ELE_CHG,
+    "kJ": 1000 / ELE_CHG
+}
+
+# length conversions
+lconvs = {
+    "angstrom": 1.0,
+    "bohr": BOHR * 1E10,
+    "nm": 10.0,
+    "m": 1E10,
+}
+
+def check_unit(unit):
+    if unit not in econvs.keys() and unit not in lconvs.keys():
+        try:
+            eunit = unit.split("/")[0]
+            lunit = unit.split("/")[1]
+            if eunit not in econvs.keys():
+                raise RuntimeError(f"Invaild unit: {unit}")
+            if lunit not in lconvs.keys():
+                raise RuntimeError(f"Invalid unit: {unit}")
+        except:
+            raise RuntimeError(f"Invalid unit: {unit}")
+
+class Conversion(ABC):
+    def __init__(self, unitA, unitB, check=True):
+        """
+        Parent class for unit conversion
+
+        Parameters
+        ----------
+            unitA : str, unit to be converted
+            unitB : str, unit which unitA is converted to, i.e. `1 unitA = self._value unitB`
+            check : bool, whether to check unit validity
+
+        Examples
+        --------
+        >>> conv = Conversion("foo", "bar", check=False)
+        >>> conv.set_value("10.0")
+        >>> print(conv)
+        1 foo = 10.0 bar
+        >>> conv.value()
+        10.0
+        """
+        if check:
+            check_unit(unitA)
+            check_unit(unitB)
+        self.unitA = unitA
+        self.unitB = unitB
+        self._value = 0.0
+
+    def value(self):
+        return self._value
+
+    def set_value(self, value):
+        self._value = value
+
+    def __repr__(self):
+        return f"1 {self.unitA} = {self._value} {self.unitB}"
+
+    def __str__(self):
+        return self.__repr__()
+
+class EnergyConversion(Conversion):
+    def __init__(self, unitA, unitB):
+        """
+        Class for energy conversion
+
+        Examples
+        --------
+        >>> conv = LengthConversion("eV", "kcal_mol")
+        >>> conv.value()
+        23.06054783061903
+        """
+        super().__init__(unitA, unitB)
+        self.set_value(econvs[unitA] / econvs[unitB])
+
+class LengthConversion(Conversion):
+    def __init__(self, unitA, unitB):
+        """
+        Class for length conversion
+
+        Examples
+        --------
+        >>> conv = LengthConversion("angstrom", "nm")
+        >>> conv.value()
+        0.1
+        """
+        super().__init__(unitA, unitB)
+        self.set_value(lconvs[unitA] / lconvs[unitB])
+
+class ForceConversion(Conversion):
+    def __init__(self, unitA, unitB):
+        """
+        Class for force conversion
+
+        Parameters
+        ----------
+            unitA, unitB : str, in format of "energy_unit/length_unit"
+
+        Examples
+        --------
+        >>> conv = ForceConversion("kJ_mol/nm", "eV/angstrom")
+        >>> conv.value()
+        0.0010364269656262175
+        """
+        super().__init__(unitA, unitB)
+        econv = EnergyConversion(unitA.split("/")[0], unitB.split("/")[0]).value()
+        lconv = LengthConversion(unitA.split("/")[1], unitB.split("/")[1]).value()
+        self.set_value(econv / lconv)
+
+class PressureConversion(Conversion):
+    def __init__(self, unitA, unitB):
+        """
+        Class for pressure conversion
+
+        Parameters
+        ----------
+            unitA, unitB : str, in format of "energy_unit/length_unit^3", or in `["Pa", "pa", "kPa", "kpa", "bar", "kbar"]`
+
+        Examples
+        --------
+        >>> conv = PressureConversion("kbar", "eV/angstrom^3")
+        >>> conv.value()
+        0.0006241509074460763
+        """
+        super().__init__(unitA, unitB, check=False)
+        unitA, factorA = self._convert_unit(unitA)
+        unitB, factorB = self._convert_unit(unitB)
+        eunitA, lunitA = self._split_unit(unitA)
+        eunitB, lunitB = self._split_unit(unitB)
+        econv = EnergyConversion(eunitA, eunitB).value() * factorA / factorB
+        lconv = LengthConversion(lunitA, lunitB).value()
+        self.set_value(econv / lconv**3)
+
+    def _convert_unit(self, unit):
+        if unit == "Pa" or unit == "pa":
+            return "J/m^3", 1
+        elif unit == "kPa" or unit == "kpa":
+            return "kJ/m^3", 1
+        elif unit == "GPa" or unit == "Gpa":
+            return "kJ/m^3", 1E6
+        elif unit == "bar":
+            return "J/m^3", 1E5
+        elif unit == "kbar":
+            return "kJ/m^3", 1E5
+        else:
+            return unit, 1
+
+    def _split_unit(self, unit):
+        eunit = unit.split("/")[0]
+        lunit = unit.split("/")[1][:-2]
+        return eunit, lunit

tests/cp2k/aimd/deepmd/energy.raw

@@ -1,3 +1,3 @@
 -7.597057812500000000e+04
--7.597049218750000000e+04
--7.597026562500000000e+04
+-7.597050000000000000e+04
+-7.597027343750000000e+04