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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
import os
import numpy as np
import xml.etree.ElementTree as ET
from ase.io.exciting import atoms2etree
from ase.units import Bohr, Hartree
from ase.calculators.calculator import PropertyNotImplementedError
from xml.dom import minidom
class Exciting:
def __init__(self, dir='calc', paramdict=None,
speciespath=None,
bin='excitingser', kpts=(1, 1, 1),
autormt=False, tshift=True, **kwargs):
"""Exciting calculator object constructor
dir: string
directory in which to execute exciting
paramdict: dict
Dictionary containing XML parameters. String values are
translated to attributes, nested dictionaries are translated
to sub elements. A list of dictionaries is translated to a
list of sub elements named after the key of which the list
is the value. Default: None
speciespath: string
Directory or URL to look up species files
bin: string
Path or executable name of exciting. Default: ``excitingser``
kpts: integer list length 3
Number of k-points
autormt: bool
Bla bla?
kwargs: dictionary like
list of key value pairs to be converted into groundstate attributes
"""
self.dir = dir
self.energy = None
self.paramdict = paramdict
if speciespath is None:
speciespath = os.environ['EXCITINGROOT'] + '/species'
self.speciespath = speciespath
self.converged = False
self.excitingbinary = bin
self.autormt = autormt
self.tshift = tshift
self.groundstate_attributes = kwargs
if ('ngridk' not in kwargs.keys() and (not (self.paramdict))):
self.groundstate_attributes['ngridk'] = ' '.join(map(str, kpts))
def update(self, atoms):
if (not self.converged or
len(self.numbers) != len(atoms) or
(self.numbers != atoms.get_atomic_numbers()).any()):
self.initialize(atoms)
self.calculate(atoms)
elif ((self.positions != atoms.get_positions()).any() or
(self.pbc != atoms.get_pbc()).any() or
(self.cell != atoms.get_cell()).any()):
self.calculate(atoms)
def initialize(self, atoms):
self.numbers = atoms.get_atomic_numbers().copy()
self.write(atoms)
def get_potential_energy(self, atoms):
"""
returns potential Energy
"""
self.update(atoms)
return self.energy
def get_forces(self, atoms):
self.update(atoms)
return self.forces.copy()
def get_stress(self, atoms):
raise PropertyNotImplementedError
def calculate(self, atoms):
self.positions = atoms.get_positions().copy()
self.cell = atoms.get_cell().copy()
self.pbc = atoms.get_pbc().copy()
self.initialize(atoms)
from pathlib import Path
xmlfile = Path(self.dir) / 'input.xml'
assert xmlfile.is_file()
print(xmlfile.read_text())
argv = [self.excitingbinary, 'input.xml']
from subprocess import check_call
check_call(argv, cwd=self.dir)
assert (Path(self.dir) / 'INFO.OUT').is_file()
assert (Path(self.dir) / 'info.xml').exists()
#syscall = ('cd %(dir)s; %(bin)s;' %
# {'dir': self.dir, 'bin': self.excitingbinary})
#print(syscall)
#assert os.system(syscall) == 0
self.read()
def write(self, atoms):
if not os.path.isdir(self.dir):
os.mkdir(self.dir)
root = atoms2etree(atoms)
root.find('structure').attrib['speciespath'] = self.speciespath
root.find('structure').attrib['autormt'] = str(self.autormt).lower()
root.find('structure').attrib['tshift'] = str(self.tshift).lower()
def prettify(elem):
rough_string = ET.tostring(elem, 'utf-8')
reparsed = minidom.parseString(rough_string)
return reparsed.toprettyxml(indent="\t")
if(self.paramdict):
self.dicttoxml(self.paramdict, root)
fd = open('%s/input.xml' % self.dir, 'w')
fd.write(prettify(root))
fd.close()
else:
groundstate = ET.SubElement(root, 'groundstate', tforce='true')
for key, value in self.groundstate_attributes.items():
if key == 'title':
root.findall('title')[0].text = value
else:
groundstate.attrib[key] = str(value)
fd = open('%s/input.xml' % self.dir, 'w')
fd.write(prettify(root))
fd.close()
def dicttoxml(self, pdict, element):
for key, value in pdict.items():
if (isinstance(value, str) and key == 'text()'):
element.text = value
elif (isinstance(value, str)):
element.attrib[key] = value
elif (isinstance(value, list)):
for item in value:
self.dicttoxml(item, ET.SubElement(element, key))
elif (isinstance(value, dict)):
if(element.findall(key) == []):
self.dicttoxml(value, ET.SubElement(element, key))
else:
self.dicttoxml(value, element.findall(key)[0])
else:
print('cannot deal with', key, '=', value)
def read(self):
"""
reads Total energy and forces from info.xml
"""
INFO_file = '%s/info.xml' % self.dir
try:
fd = open(INFO_file)
except IOError:
raise RuntimeError("output doesn't exist")
info = ET.parse(fd)
self.energy = float(info.findall(
'groundstate/scl/iter/energies')[-1].attrib['totalEnergy']) * Hartree
forces = []
forcesnodes = info.findall(
'groundstate/scl/structure')[-1].findall('species/atom/forces/totalforce')
for force in forcesnodes:
forces.append(np.array(list(force.attrib.values())).astype(float))
self.forces = np.reshape(forces, (-1, 3)) * Hartree / Bohr
if str(info.find('groundstate').attrib['status']) == 'finished':
self.converged = True
else:
raise RuntimeError('calculation did not finish correctly')
|
