Skip to content

Commit 23e9f3b

Browse files
Ericwang6Ericwang6
authored
Add unit.py to manage unit conversions (#202)
1 parent 3c595ca commit 23e9f3b

File tree

18 files changed

+241
-84
lines changed

18 files changed

+241
-84
lines changed

‎dpdata/abacus/scf.py

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,10 @@
11
import os,sys
22
import numpy as np
3+
from ..unit import EnergyConversion, PressureConversion, LengthConversion
34

4-
bohr2ang = 0.5291770
5-
ry2ev = 13.605698
6-
kbar2evperang3 = 1e3 / 1.6021892e6
7-
# The consts are cited from $ABACUS_ROOT/source/src_global/constant.h
8-
5+
bohr2ang = LengthConversion("bohr", "angstrom").value()
6+
ry2ev = EnergyConversion("rydberg", "eV").value()
7+
kbar2evperang3 = PressureConversion("kbar", "eV/angstrom^3").value()
98

109
def get_block (lines, keyword, skip = 0, nlines = None):
1110
ret = []

‎dpdata/amber/md.py

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,11 @@
33
from scipy.io import netcdf
44
import numpy as np
55
from dpdata.amber.mask import pick_by_amber_mask
6+
from dpdata.unit import EnergyConversion
7+
from ..periodic_table import ELEMENTS
68

7-
kcalmol2eV= 0.04336410390059322
8-
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']
9+
kcalmol2eV = EnergyConversion("kcal_mol", "eV").value()
10+
symbols = ['X'] + ELEMENTS
911

1012
energy_convert = kcalmol2eV
1113
force_convert = energy_convert

‎dpdata/amber/sqm.py

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
import numpy as np
2-
from ..periodic_table import ELEMENTS
2+
from dpdata.periodic_table import ELEMENTS
3+
from dpdata.unit import EnergyConversion
34

5+
kcal2ev = EnergyConversion("kcal_mol", "eV").value()
46

57
START = 0
68
READ_ENERGY = 1
@@ -47,7 +49,8 @@ def parse_sqm_out(fname):
4749
if len(coords) == len(charges):
4850
flag = START
4951
elif flag == READ_FORCES:
50-
forces.append([float(x) for x in line.strip().split()[-3:]])
52+
ll = line.strip()
53+
forces.append([float(ll[-60:-40]), float(ll[-40:-20]), float(ll[-20:])])
5154
if len(forces) == len(charges):
5255
flag = START
5356

@@ -62,7 +65,7 @@ def parse_sqm_out(fname):
6265
data['coords'] = np.array([coords])
6366

6467
energies = np.array(energies)
65-
forces = np.array([forces], dtype=np.float32)
68+
forces = np.array([forces], dtype=np.float32) * kcal2ev
6669
if len(forces) > 0:
6770
data['energies'] = energies
6871
data['forces'] = forces
@@ -99,4 +102,4 @@ def make_sqm_in(data, fname=None, frame_idx=0, **kwargs):
99102
if fname is not None:
100103
with open(fname, 'w') as fp:
101104
fp.write(ret)
102-
return ret
105+
return ret

‎dpdata/cp2k/output.py

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,12 @@
33
import re
44
from collections import OrderedDict
55
from .cell import cell_to_low_triangle
6+
from ..unit import EnergyConversion, LengthConversion, ForceConversion, PressureConversion
67

78
#%%
8-
AU_TO_ANG = 5.29177208590000E-01
9-
AU_TO_EV = 2.72113838565563E+01
10-
AU_TO_EV_EVERY_ANG = AU_TO_EV/AU_TO_ANG
9+
AU_TO_ANG = LengthConversion("bohr", "angstrom").value()
10+
AU_TO_EV = EnergyConversion("hartree", "eV").value()
11+
AU_TO_EV_EVERY_ANG = ForceConversion("hartree/bohr", "eV/angstrom").value()
1112
delimiter_patterns=[]
1213
delimiter_p1 = re.compile(r'^ \* GO CP2K GO! \*+')
1314
delimiter_p2 = re.compile(r'^ \*+')
@@ -221,9 +222,9 @@ def get_frames (fname) :
221222
coord_flag = False
222223
force_flag = False
223224
stress_flag = False
224-
eV = 2.72113838565563E+01 # hatree to eV
225-
angstrom = 5.29177208590000E-01 # Bohr to Angstrom
226-
GPa = 160.21766208 # 1 eV/(Angstrom^3) = 160.21 GPa
225+
eV = EnergyConversion("hartree", "eV").value()
226+
angstrom = LengthConversion("bohr", "angstrom").value()
227+
GPa = PressureConversion("eV/angstrom^3", "GPa").value()
227228
atom_symbol_list = []
228229
cell = []
229230
coord = []

‎dpdata/gaussian/log.py

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,12 @@
11
import numpy as np
2+
from ..unit import LengthConversion, EnergyConversion, ForceConversion
3+
from ..periodic_table import ELEMENTS
24

5+
length_convert = LengthConversion("bohr", "angstrom").value()
6+
energy_convert = EnergyConversion("hartree", "eV").value()
7+
force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
38

4-
hartree2ev = 27.211386018
5-
bohr2ang = 0.52917721067
6-
7-
length_convert = bohr2ang
8-
energy_convert = hartree2ev
9-
force_convert = energy_convert / length_convert
10-
11-
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']
9+
symbols = ["X"] + ELEMENTS
1210

1311
def to_system_data(file_name, md=False):
1412
data = {}

‎dpdata/gromacs/gro.py

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,10 @@
11
#!/usr/bin/env python3
22
import re
33
import numpy as np
4-
import json
5-
import os
4+
from ..unit import LengthConversion
65

7-
nm2ang = 10.
8-
ang2nm = 1. / nm2ang
6+
nm2ang = LengthConversion("nm", "angstrom").value()
7+
ang2nm = LengthConversion("angstrom", "nm").value()
98
cell_idx_gmx2dp = [0, 4, 8, 1, 2, 3, 5, 6, 7]
109

1110
def _format_atom_name(atom_name):

‎dpdata/pwmat/atomconfig.py

Lines changed: 1 addition & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#!/usr/bin/python3
2-
2+
from ..periodic_table import ELEMENTS
33
import numpy as np
44

55
def _to_system_data_lower(lines) :
@@ -39,13 +39,6 @@ def _to_system_data_lower(lines) :
3939
for jj in range(ii) :
4040
atom_types.append(idx)
4141
system['atom_types'] = np.array(atom_types, dtype = int)
42-
ELEMENTS=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', \
43-
'Sc', 'Ti', 'V', 'Cr','Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', \
44-
'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',\
45-
'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', \
46-
'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', \
47-
'Md', 'No', 'Lr']
48-
4942
system['atom_names'] = [ELEMENTS[ii-1] for ii in np.unique(sorted(atomic_number))]
5043
return system
5144

@@ -80,12 +73,6 @@ def from_system_data(system, f_idx = 0, skip_zeros = True) :
8073
for ii, jj in zip(atom_numbs, atom_names):
8174
for kk in range(ii):
8275
symbal.append(jj)
83-
ELEMENTS=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', \
84-
'Sc', 'Ti', 'V', 'Cr','Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', \
85-
'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd',\
86-
'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', \
87-
'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', \
88-
'Md', 'No', 'Lr']
8976
atomic_numbers = []
9077
for ii in symbal:
9178
atomic_numbers.append(ELEMENTS.index(ii)+1)

‎dpdata/pwmat/movement.py

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
import numpy as np
2+
from ..periodic_table import ELEMENTS
23

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

‎dpdata/qe/traj.py

Lines changed: 7 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,14 @@
11
#!/usr/bin/python3
2-
32
import numpy as np
4-
import dpdata,warnings
3+
import dpdata, warnings
4+
from ..unit import EnergyConversion, LengthConversion, ForceConversion, PressureConversion
55

6-
ry2ev = 13.605693009
7-
hartree2ev = 27.211386018
8-
bohr2ang = 0.52917721067
9-
kbar2evperang3 = 1e3 / 1.602176621e6
6+
ry2ev = EnergyConversion("rydberg", "eV").value()
7+
kbar2evperang3 = PressureConversion("kbar", "eV/angstrom^3").value()
108

11-
length_convert = bohr2ang
12-
energy_convert = hartree2ev
13-
force_convert = energy_convert / length_convert
9+
length_convert = LengthConversion("bohr", "angstrom").value()
10+
energy_convert = EnergyConversion("hartree", "eV").value()
11+
force_convert = ForceConversion("hartree/bohr", "eV/angstrom").value()
1412

1513
def load_key (lines, key) :
1614
for ii in lines :

‎dpdata/unit.py

Lines changed: 169 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,169 @@
1+
from abc import ABC
2+
from scipy import constants
3+
4+
AVOGADRO = constants.Avogadro # Avagadro constant
5+
ELE_CHG = constants.elementary_charge # Elementary Charge, in C
6+
BOHR = constants.value("atomic unit of length") # Bohr, in m
7+
HARTREE = constants.value("atomic unit of energy") # Hartree, in Jole
8+
RYDBERG = constants.Rydberg * constants.h * constants.c # Rydberg, in Jole
9+
10+
# energy conversions
11+
econvs = {
12+
"eV": 1.0,
13+
"hartree": HARTREE / ELE_CHG,
14+
"kJ_mol": 1 / (ELE_CHG * AVOGADRO / 1000),
15+
"kcal_mol": 1 / (ELE_CHG * AVOGADRO / 1000 / 4.184),
16+
"rydberg": RYDBERG / ELE_CHG,
17+
"J": 1 / ELE_CHG,
18+
"kJ": 1000 / ELE_CHG
19+
}
20+
21+
# length conversions
22+
lconvs = {
23+
"angstrom": 1.0,
24+
"bohr": BOHR * 1E10,
25+
"nm": 10.0,
26+
"m": 1E10,
27+
}
28+
29+
def check_unit(unit):
30+
if unit not in econvs.keys() and unit not in lconvs.keys():
31+
try:
32+
eunit = unit.split("/")[0]
33+
lunit = unit.split("/")[1]
34+
if eunit not in econvs.keys():
35+
raise RuntimeError(f"Invaild unit: {unit}")
36+
if lunit not in lconvs.keys():
37+
raise RuntimeError(f"Invalid unit: {unit}")
38+
except:
39+
raise RuntimeError(f"Invalid unit: {unit}")
40+
41+
class Conversion(ABC):
42+
def __init__(self, unitA, unitB, check=True):
43+
"""
44+
Parent class for unit conversion
45+
46+
Parameters
47+
----------
48+
unitA : str, unit to be converted
49+
unitB : str, unit which unitA is converted to, i.e. `1 unitA = self._value unitB`
50+
check : bool, whether to check unit validity
51+
52+
Examples
53+
--------
54+
>>> conv = Conversion("foo", "bar", check=False)
55+
>>> conv.set_value("10.0")
56+
>>> print(conv)
57+
1 foo = 10.0 bar
58+
>>> conv.value()
59+
10.0
60+
"""
61+
if check:
62+
check_unit(unitA)
63+
check_unit(unitB)
64+
self.unitA = unitA
65+
self.unitB = unitB
66+
self._value = 0.0
67+
68+
def value(self):
69+
return self._value
70+
71+
def set_value(self, value):
72+
self._value = value
73+
74+
def __repr__(self):
75+
return f"1 {self.unitA} = {self._value} {self.unitB}"
76+
77+
def __str__(self):
78+
return self.__repr__()
79+
80+
class EnergyConversion(Conversion):
81+
def __init__(self, unitA, unitB):
82+
"""
83+
Class for energy conversion
84+
85+
Examples
86+
--------
87+
>>> conv = LengthConversion("eV", "kcal_mol")
88+
>>> conv.value()
89+
23.06054783061903
90+
"""
91+
super().__init__(unitA, unitB)
92+
self.set_value(econvs[unitA] / econvs[unitB])
93+
94+
class LengthConversion(Conversion):
95+
def __init__(self, unitA, unitB):
96+
"""
97+
Class for length conversion
98+
99+
Examples
100+
--------
101+
>>> conv = LengthConversion("angstrom", "nm")
102+
>>> conv.value()
103+
0.1
104+
"""
105+
super().__init__(unitA, unitB)
106+
self.set_value(lconvs[unitA] / lconvs[unitB])
107+
108+
class ForceConversion(Conversion):
109+
def __init__(self, unitA, unitB):
110+
"""
111+
Class for force conversion
112+
113+
Parameters
114+
----------
115+
unitA, unitB : str, in format of "energy_unit/length_unit"
116+
117+
Examples
118+
--------
119+
>>> conv = ForceConversion("kJ_mol/nm", "eV/angstrom")
120+
>>> conv.value()
121+
0.0010364269656262175
122+
"""
123+
super().__init__(unitA, unitB)
124+
econv = EnergyConversion(unitA.split("/")[0], unitB.split("/")[0]).value()
125+
lconv = LengthConversion(unitA.split("/")[1], unitB.split("/")[1]).value()
126+
self.set_value(econv / lconv)
127+
128+
class PressureConversion(Conversion):
129+
def __init__(self, unitA, unitB):
130+
"""
131+
Class for pressure conversion
132+
133+
Parameters
134+
----------
135+
unitA, unitB : str, in format of "energy_unit/length_unit^3", or in `["Pa", "pa", "kPa", "kpa", "bar", "kbar"]`
136+
137+
Examples
138+
--------
139+
>>> conv = PressureConversion("kbar", "eV/angstrom^3")
140+
>>> conv.value()
141+
0.0006241509074460763
142+
"""
143+
super().__init__(unitA, unitB, check=False)
144+
unitA, factorA = self._convert_unit(unitA)
145+
unitB, factorB = self._convert_unit(unitB)
146+
eunitA, lunitA = self._split_unit(unitA)
147+
eunitB, lunitB = self._split_unit(unitB)
148+
econv = EnergyConversion(eunitA, eunitB).value() * factorA / factorB
149+
lconv = LengthConversion(lunitA, lunitB).value()
150+
self.set_value(econv / lconv**3)
151+
152+
def _convert_unit(self, unit):
153+
if unit == "Pa" or unit == "pa":
154+
return "J/m^3", 1
155+
elif unit == "kPa" or unit == "kpa":
156+
return "kJ/m^3", 1
157+
elif unit == "GPa" or unit == "Gpa":
158+
return "kJ/m^3", 1E6
159+
elif unit == "bar":
160+
return "J/m^3", 1E5
161+
elif unit == "kbar":
162+
return "kJ/m^3", 1E5
163+
else:
164+
return unit, 1
165+
166+
def _split_unit(self, unit):
167+
eunit = unit.split("/")[0]
168+
lunit = unit.split("/")[1][:-2]
169+
return eunit, lunit

0 commit comments

Comments
 (0)