1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
import re
import numpy as np
from ase import Atoms
from ase.geometry import cellpar_to_cell
from .parser import _define_pattern
# Geometry block parser
_geom = _define_pattern(
r'^[ \t]*geometry[ \t\S]*\n'
r'((?:^[ \t]*[\S]+[ \t\S]*\n)+)'
r'^[ \t]*end\n\n',
"""\
geometry units angstrom nocenter noautosym noautoz
system crystal units angstrom
lattice_vectors
4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
end
O 5.0000000000000000e-01 5.0000000000000011e-01 5.6486824536818558e-01
H 5.0000000000000000e-01 6.3810586054988372e-01 4.3513175463181430e-01
H 5.0000000000000000e-01 3.6189413945011639e-01 4.3513175463181430e-01
end
""", re.M)
# Finds crystal specification
_crystal = _define_pattern(
r'^[ \t]*system crystal[ \t\S]*\n'
r'((?:[ \t]*[\S]+[ \t\S]*\n)+?)'
r'^[ \t]*end[ \t]*\n',
"""\
system crystal units angstrom
lattice_vectors
4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
end
""", re.M)
# Finds 3d-periodic unit cell
_cell_3d = _define_pattern(
r'^[ \t]*lattice_vectors[ \t]*\n'
r'^((?:(?:[ \t]+[\S]+){3}\n){3})',
"""\
lattice_vectors
4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
""", re.M)
# Extracts chemical species from a geometry block
_species = _define_pattern(
r'^[ \t]*[A-Z][a-z]?(?:[ \t]+[\S]+){3}\n',
" O 0.0 0.0 0.0\n", re.M,
)
def read_nwchem_in(fobj, index=-1):
text = ''.join(fobj.readlines())
atomslist = []
for match in _geom.findall(text):
symbols = []
positions = []
for atom in _species.findall(match):
atom = atom.split()
symbols.append(atom[0])
positions.append([float(x) for x in atom[1:]])
positions = np.array(positions)
atoms = Atoms(symbols)
cell, pbc = _get_cell(text)
pos = np.zeros_like(positions)
for dim, ipbc in enumerate(pbc):
if ipbc:
pos += np.outer(positions[:, dim], cell[dim, :])
else:
pos[:, dim] = positions[:, dim]
atoms.set_cell(cell)
atoms.pbc = pbc
atoms.set_positions(pos)
atomslist.append(atoms)
return atomslist[index]
def _get_cell(text):
# first check whether there is a lattice definition
cell = np.zeros((3, 3))
lattice = _cell_3d.findall(text)
if lattice:
pbc = [True, True, True]
for i, row in enumerate(lattice[0].strip().split('\n')):
cell[i] = [float(x) for x in row.split()]
return cell, pbc
pbc = [False, False, False]
lengths = [None, None, None]
angles = [None, None, None]
for row in text.strip().split('\n'):
row = row.strip().lower()
for dim, vecname in enumerate(['a', 'b', 'c']):
if row.startswith('lat_{}'.format(vecname)):
pbc[dim] = True
lengths[dim] = float(row.split()[1])
for i, angle in enumerate(['alpha', 'beta', 'gamma']):
if row.startswith(angle):
angles[i] = float(row.split()[1])
if not np.any(pbc):
return None, pbc
for i in range(3):
a, b, c = np.roll(np.array([0, 1, 2]), i)
if pbc[a] and pbc[b]:
assert angles[c] is not None
if angles[c] is not None:
assert pbc[a] and pbc[b]
# The easiest case: all three lattice vectors and angles are specified
if np.all(pbc):
return cellpar_to_cell(lengths + angles), pbc
# Next easiest case: exactly one lattice vector has been specified
if np.sum(pbc) == 1:
dim = np.argmax(pbc)
cell[dim, dim] = lengths[dim]
return cell, pbc
# Hardest case: two lattice vectors are specified.
dim1, dim2 = [dim for dim, ipbc in enumerate(pbc) if ipbc]
angledim = np.argmin(pbc)
cell[dim1, dim1] = lengths[dim1]
cell[dim2, dim2] = lengths[dim2] * np.sin(angles[angledim])
cell[dim2, dim1] = lengths[dim2] * np.cos(angles[angledim])
return cell, pbc
|
