1import os 

2 

3import numpy as np 

4import xml.etree.ElementTree as ET 

5from ase.io.exciting import atoms2etree 

6from ase.units import Bohr, Hartree 

7from ase.calculators.calculator import PropertyNotImplementedError 

8from xml.dom import minidom 

9 

10 

11class Exciting: 

12 def __init__(self, dir='calc', paramdict=None, 

13 speciespath=None, 

14 bin='excitingser', kpts=(1, 1, 1), 

15 autormt=False, tshift=True, **kwargs): 

16 """Exciting calculator object constructor 

17 

18 dir: string 

19 directory in which to execute exciting 

20 paramdict: dict 

21 Dictionary containing XML parameters. String values are 

22 translated to attributes, nested dictionaries are translated 

23 to sub elements. A list of dictionaries is translated to a 

24 list of sub elements named after the key of which the list 

25 is the value. Default: None 

26 speciespath: string 

27 Directory or URL to look up species files 

28 bin: string 

29 Path or executable name of exciting. Default: ``excitingser`` 

30 kpts: integer list length 3 

31 Number of k-points 

32 autormt: bool 

33 Bla bla? 

34 kwargs: dictionary like 

35 list of key value pairs to be converted into groundstate attributes 

36 

37 """ 

38 self.dir = dir 

39 self.energy = None 

40 

41 self.paramdict = paramdict 

42 if speciespath is None: 

43 speciespath = os.environ['EXCITINGROOT'] + '/species' 

44 self.speciespath = speciespath 

45 self.converged = False 

46 self.excitingbinary = bin 

47 self.autormt = autormt 

48 self.tshift = tshift 

49 self.groundstate_attributes = kwargs 

50 if ('ngridk' not in kwargs.keys() and (not (self.paramdict))): 

51 self.groundstate_attributes['ngridk'] = ' '.join(map(str, kpts)) 

52 

53 def update(self, atoms): 

54 if (not self.converged or 

55 len(self.numbers) != len(atoms) or 

56 (self.numbers != atoms.get_atomic_numbers()).any()): 

57 self.initialize(atoms) 

58 self.calculate(atoms) 

59 elif ((self.positions != atoms.get_positions()).any() or 

60 (self.pbc != atoms.get_pbc()).any() or 

61 (self.cell != atoms.get_cell()).any()): 

62 self.calculate(atoms) 

63 

64 def initialize(self, atoms): 

65 self.numbers = atoms.get_atomic_numbers().copy() 

66 self.write(atoms) 

67 

68 def get_potential_energy(self, atoms): 

69 """Returns potential Energy.""" 

70 self.update(atoms) 

71 return self.energy 

72 

73 def get_forces(self, atoms): 

74 self.update(atoms) 

75 return self.forces.copy() 

76 

77 def get_stress(self, atoms): 

78 raise PropertyNotImplementedError 

79 

80 def calculate(self, atoms): 

81 self.positions = atoms.get_positions().copy() 

82 self.cell = atoms.get_cell().copy() 

83 self.pbc = atoms.get_pbc().copy() 

84 

85 self.initialize(atoms) 

86 from pathlib import Path 

87 xmlfile = Path(self.dir) / 'input.xml' 

88 assert xmlfile.is_file() 

89 print(xmlfile.read_text()) 

90 argv = [self.excitingbinary, 'input.xml'] 

91 from subprocess import check_call 

92 check_call(argv, cwd=self.dir) 

93 

94 assert (Path(self.dir) / 'INFO.OUT').is_file() 

95 assert (Path(self.dir) / 'info.xml').exists() 

96 

97 # syscall = ('cd %(dir)s; %(bin)s;' % 

98 # {'dir': self.dir, 'bin': self.excitingbinary}) 

99 # print(syscall) 

100 # assert os.system(syscall) == 0 

101 self.read() 

102 

103 def write(self, atoms): 

104 if not os.path.isdir(self.dir): 

105 os.mkdir(self.dir) 

106 root = atoms2etree(atoms) 

107 root.find('structure').attrib['speciespath'] = self.speciespath 

108 root.find('structure').attrib['autormt'] = str(self.autormt).lower() 

109 root.find('structure').attrib['tshift'] = str(self.tshift).lower() 

110 

111 def prettify(elem): 

112 rough_string = ET.tostring(elem, 'utf-8') 

113 reparsed = minidom.parseString(rough_string) 

114 return reparsed.toprettyxml(indent="\t") 

115 

116 if self.paramdict: 

117 self.dicttoxml(self.paramdict, root) 

118 fd = open('%s/input.xml' % self.dir, 'w') 

119 fd.write(prettify(root)) 

120 fd.close() 

121 else: 

122 groundstate = ET.SubElement(root, 'groundstate', tforce='true') 

123 for key, value in self.groundstate_attributes.items(): 

124 if key == 'title': 

125 root.findall('title')[0].text = value 

126 else: 

127 groundstate.attrib[key] = str(value) 

128 fd = open('%s/input.xml' % self.dir, 'w') 

129 fd.write(prettify(root)) 

130 fd.close() 

131 

132 def dicttoxml(self, pdict, element): 

133 for key, value in pdict.items(): 

134 if (isinstance(value, str) and key == 'text()'): 

135 element.text = value 

136 elif (isinstance(value, str)): 

137 element.attrib[key] = value 

138 elif (isinstance(value, list)): 

139 for item in value: 

140 self.dicttoxml(item, ET.SubElement(element, key)) 

141 elif (isinstance(value, dict)): 

142 if element.findall(key) == []: 

143 self.dicttoxml(value, ET.SubElement(element, key)) 

144 else: 

145 self.dicttoxml(value, element.findall(key)[0]) 

146 else: 

147 print('cannot deal with', key, '=', value) 

148 

149 def read(self): 

150 """ 

151 reads Total energy and forces from info.xml 

152 """ 

153 INFO_file = '%s/info.xml' % self.dir 

154 

155 try: 

156 fd = open(INFO_file) 

157 except IOError: 

158 raise RuntimeError("output doesn't exist") 

159 info = ET.parse(fd) 

160 self.energy = float(info.findall( 

161 'groundstate/scl/iter/energies')[-1].attrib[ 

162 'totalEnergy']) * Hartree 

163 forces = [] 

164 forcesnodes = info.findall( 

165 'groundstate/scl/structure')[-1].findall( 

166 'species/atom/forces/totalforce') 

167 for force in forcesnodes: 

168 forces.append(np.array(list(force.attrib.values())).astype(float)) 

169 self.forces = np.reshape(forces, (-1, 3)) * Hartree / Bohr 

170 

171 if str(info.find('groundstate').attrib['status']) == 'finished': 

172 self.converged = True 

173 else: 

174 raise RuntimeError('calculation did not finish correctly')