Skip to content

Commit 2bec372

Browse files
felix5572felix5572
committed
cell bug fix 2
1 parent ba28695 commit 2bec372

File tree

4 files changed

+82
-59
lines changed

4 files changed

+82
-59
lines changed

‎dpdata/cp2k/__init__.py

Whitespace-only changes.

‎dpdata/cp2k/cell.py

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
2+
#%%
3+
import numpy as np
4+
from collections import OrderedDict
5+
import re
6+
7+
def cell_to_low_triangle(A,B,C,alpha,beta,gamma):
8+
"""
9+
Convert cell to low triangle matrix.
10+
11+
Parameters
12+
----------
13+
A : float
14+
cell length A
15+
B : float
16+
cell length B
17+
C : float
18+
cell length C
19+
alpha : float
20+
radian. The angle between vector B and vector C.
21+
beta : float
22+
radian. The angle between vector A and vector C.
23+
gamma : float
24+
radian. The angle between vector B and vector C.
25+
26+
Returns
27+
-------
28+
cell : list
29+
The cell matrix used by dpdata in low triangle form.
30+
"""
31+
if not np.pi*5/180<alpha< np.pi*175/180:
32+
raise RuntimeError("alpha=={}: must be a radian, and \
33+
must be in np.pi*5/180 < alpha < np.pi*175/180".format(alpha))
34+
if not np.pi*5/180<beta< np.pi*175/180:
35+
raise RuntimeError("beta=={}: must be a radian, and \
36+
must be in np.pi*5/180 < beta < np.pi*175/180".format(beta))
37+
if not np.pi*5/180<gamma< np.pi*175/180:
38+
raise RuntimeError("gamma=={}: must be a radian, and \
39+
must be in np.pi*5/180 < gamma < np.pi*175/180".format(gamma))
40+
if not A > 0.2:
41+
raise RuntimeError("A=={}, must be greater than 0.2".format(A))
42+
if not B > 0.2:
43+
raise RuntimeError("B=={}, must be greater than 0.2".format(B))
44+
if not C > 0.2:
45+
raise RuntimeError("C=={}, must be greater than 0.2".format(C))
46+
47+
lx = A
48+
xy = B * np.cos(gamma)
49+
xz = C * np.cos(beta)
50+
ly = B* np.sin(gamma)
51+
if not ly > 0.1:
52+
raise RuntimeError("ly:=B* np.sin(gamma)=={}, must be greater than 0.1",format(ly))
53+
yz = (B*C*np.cos(alpha)-xy*xz)/ly
54+
if not C**2-xz**2-yz**2 > 0.01:
55+
raise RuntimeError("lz^2:=C**2-xz**2-yz**2=={}, must be greater than 0.01",format(C**2-xz**2-yz**2))
56+
lz = np.sqrt(C**2-xz**2-yz**2)
57+
cell = np.asarray([[lx, 0 , 0],
58+
[xy, ly, 0 ],
59+
[xz, yz, lz]]).astype('float32')
60+
return cell
61+

‎dpdata/cp2k/output.py

Lines changed: 2 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
1-
21
#%%
32
import numpy as np
4-
from collections import OrderedDict
53
import re
4+
from collections import OrderedDict
5+
from .cell import cell_to_low_triangle
66

77
#%%
88
AU_TO_ANG = 5.29177208590000E-01
@@ -18,58 +18,6 @@
1818
avail_patterns.append(re.compile(r'^ INITIAL POTENTIAL ENERGY'))
1919
avail_patterns.append(re.compile(r'^ ENSEMBLE TYPE'))
2020

21-
def cell_to_low_triangle(A,B,C,alpha,beta,gamma):
22-
"""
23-
Convert cell to low triangle matrix.
24-
25-
Parameters
26-
----------
27-
A : float
28-
cell length A
29-
B : float
30-
cell length B
31-
C : float
32-
cell length C
33-
alpha : float
34-
radian. The angle between vector B and vector C.
35-
beta : float
36-
radian. The angle between vector A and vector C.
37-
gamma : float
38-
radian. The angle between vector B and vector C.
39-
40-
Returns
41-
-------
42-
cell : list
43-
The cell matrix used by dpdata in low triangle form.
44-
"""
45-
if not np.pi*5/180<alpha< np.pi*175/180:
46-
raise RuntimeError("alpha=={}: must be a radian, and \
47-
must be in np.pi*5/180 < alpha < np.pi*175/180".format(alpha))
48-
if not np.pi*5/180<beta< np.pi*175/180:
49-
raise RuntimeError("beta=={}: must be a radian, and \
50-
must be in np.pi*5/180 < beta < np.pi*175/180".format(beta))
51-
if not np.pi*5/180<gamma< np.pi*175/180:
52-
raise RuntimeError("gamma=={}: must be a radian, and \
53-
must be in np.pi*5/180 < gamma < np.pi*175/180".format(gamma))
54-
if not A > 0.2:
55-
raise RuntimeError("A=={}, must be greater than 0.2".format(A))
56-
if not B > 0.2:
57-
raise RuntimeError("B=={}, must be greater than 0.2".format(B))
58-
if not C > 0.2:
59-
raise RuntimeError("C=={}, must be greater than 0.2".format(C))
60-
61-
lx = A
62-
xy = B * np.cos(gamma)
63-
xz = C * np.cos(beta)
64-
ly = B* np.sin(gamma)
65-
if not ly > 0.1:
66-
raise RuntimeError("ly:=B* np.sin(gamma)=={}, must be greater than 0.1",format(ly))
67-
yz = (B*C*np.cos(alpha)-xy*xz)/ly
68-
lz = np.sqrt(C**2-xz**2-yz**2)
69-
cell = np.asarray([[lx, 0 , 0],
70-
[xy, ly, 0 ],
71-
[xz, yz, lz]]).astype('float32')
72-
return cell
7321
class Cp2kSystems(object):
7422
"""
7523
deal with cp2k outputfile

‎tests/test_cell_to_low_triangle.py

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,14 @@
55

66
class TestCellToLowTriangle(unittest.TestCase):
77
def test_func1(self):
8-
cell_1 = dpdata.cp2k.output.cell_to_low_triangle(6,6,6,np.pi*1/2, np.pi*1/2, np.pi*1/2)
8+
cell_1 = dpdata.cp2k.cell.cell_to_low_triangle(6,6,6,np.pi*1/2, np.pi*1/2, np.pi*1/2)
99
cell_2 = np.asarray([[6,0,0],[0,6,0],[0,0,6]])
1010
for ii in range(3):
1111
for jj in range(3):
1212
self.assertAlmostEqual(cell_1[ii,jj], cell_2[ii,jj], places=6)
1313

1414
def test_func2(self):
15-
cell_1 = dpdata.cp2k.output.cell_to_low_triangle(6,6,6,np.pi*1/3, np.pi*1/3, np.pi*1/3)
15+
cell_1 = dpdata.cp2k.cell.cell_to_low_triangle(6,6,6,np.pi*1/3, np.pi*1/3, np.pi*1/3)
1616
cell_2 = np.asarray([
1717
[6,0,0],
1818
[3,3*np.sqrt(3),0],
@@ -22,14 +22,28 @@ def test_func2(self):
2222
self.assertAlmostEqual(cell_1[ii,jj], cell_2[ii,jj], places=6)
2323

2424
def test_func3(self):
25+
cell_1 = dpdata.cp2k.cell.cell_to_low_triangle(6,7,8,np.pi*133/180,np.pi*84/180,np.pi*69/180)
26+
cell_2 = np.asarray([[ 6.0, 0.0, 0.0],
27+
[ 2.5085757, 6.535063 , 0.0],
28+
[ 0.8362277, -6.1651506, 5.0290794]], dtype='float32')
29+
for ii in range(3):
30+
for jj in range(3):
31+
self.assertAlmostEqual(cell_1[ii,jj], cell_2[ii,jj], places=6)
32+
33+
def test_func4(self):
2534
with self.assertRaises(Exception) as c:
26-
dpdata.cp2k.output.cell_to_low_triangle(0.1,6,6,np.pi*1/2,np.pi*1/2,np.pi*1/2)
35+
dpdata.cp2k.cell.cell_to_low_triangle(0.1,6,6,np.pi*1/2,np.pi*1/2,np.pi*1/2)
2736
self.assertTrue("A==0.1" in str(c.exception))
2837

29-
def test_func4(self):
38+
def test_func5(self):
3039
with self.assertRaises(Exception) as c:
31-
dpdata.cp2k.output.cell_to_low_triangle(6,6,6,np.pi*3/180,np.pi*1/2,np.pi*1/2)
40+
dpdata.cp2k.cell.cell_to_low_triangle(6,6,6,np.pi*3/180,np.pi*1/2,np.pi*1/2)
3241
self.assertTrue("alpha" in str(c.exception))
3342

43+
def test_func6(self):
44+
with self.assertRaises(Exception) as c:
45+
dpdata.cp2k.cell.cell_to_low_triangle(6,7,8,np.pi*153/180,np.pi*84/180,np.pi*69/180)
46+
self.assertTrue("lz^2" in str(c.exception))
47+
3448
if __name__ == '__main__':
3549
unittest.main()

0 commit comments

Comments
 (0)