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 """ 

70 returns potential Energy 

71 """ 

72 self.update(atoms) 

73 return self.energy 

74 

75 def get_forces(self, atoms): 

76 self.update(atoms) 

77 return self.forces.copy() 

78 

79 def get_stress(self, atoms): 

80 raise PropertyNotImplementedError 

81 

82 def calculate(self, atoms): 

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

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

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

86 

87 self.initialize(atoms) 

88 from pathlib import Path 

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

90 assert xmlfile.is_file() 

91 print(xmlfile.read_text()) 

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

93 from subprocess import check_call 

94 check_call(argv, cwd=self.dir) 

95 

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

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

98 

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

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

101 #print(syscall) 

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

103 self.read() 

104 

105 def write(self, atoms): 

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

107 os.mkdir(self.dir) 

108 root = atoms2etree(atoms) 

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

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

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

112 

113 def prettify(elem): 

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

115 reparsed = minidom.parseString(rough_string) 

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

117 

118 if(self.paramdict): 

119 self.dicttoxml(self.paramdict, root) 

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

121 fd.write(prettify(root)) 

122 fd.close() 

123 else: 

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

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

126 if key == 'title': 

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

128 else: 

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

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

131 fd.write(prettify(root)) 

132 fd.close() 

133 

134 def dicttoxml(self, pdict, element): 

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

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

137 element.text = value 

138 elif (isinstance(value, str)): 

139 element.attrib[key] = value 

140 elif (isinstance(value, list)): 

141 for item in value: 

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

143 elif (isinstance(value, dict)): 

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

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

146 else: 

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

148 else: 

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

150 

151 def read(self): 

152 """ 

153 reads Total energy and forces from info.xml 

154 """ 

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

156 

157 try: 

158 fd = open(INFO_file) 

159 except IOError: 

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

161 info = ET.parse(fd) 

162 self.energy = float(info.findall( 

163 'groundstate/scl/iter/energies')[-1].attrib['totalEnergy']) * Hartree 

164 forces = [] 

165 forcesnodes = info.findall( 

166 'groundstate/scl/structure')[-1].findall('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')