1"""This module defines an ASE interface to FHI-aims. 

2 

3Felix Hanke hanke@liverpool.ac.uk 

4Jonas Bjork j.bjork@liverpool.ac.uk 

5Simon P. Rittmeyer simon.rittmeyer@tum.de 

6 

7Edits on (24.11.2021) by Thomas A. R. Purcell purcell@fhi-berlin.mpg.de 

8""" 

9 

10import os 

11import re 

12 

13import numpy as np 

14 

15from ase.io.aims import write_aims, write_control 

16from ase.calculators.genericfileio import (GenericFileIOCalculator, 

17 CalculatorTemplate) 

18 

19 

20def get_aims_version(string): 

21 match = re.search(r"\s*FHI-aims version\s*:\s*(\S+)", string, re.M) 

22 return match.group(1) 

23 

24 

25class AimsProfile: 

26 def __init__(self, argv): 

27 if isinstance(argv, str): 

28 argv = argv.split() 

29 

30 self.argv = argv 

31 

32 def run(self, directory, outputname): 

33 from subprocess import check_call 

34 

35 with open(directory / outputname, "w") as fd: 

36 check_call(self.argv, stdout=fd, cwd=directory, 

37 env=os.environ) 

38 

39 def socketio_argv_unix(self, socket): 

40 return list(self.argv) 

41 

42 def socketio_argv_inet(self, port): 

43 return list(self.argv) 

44 

45 

46class AimsTemplate(CalculatorTemplate): 

47 def __init__(self): 

48 super().__init__( 

49 "aims", 

50 [ 

51 "energy", 

52 "free_energy", 

53 "forces", 

54 "stress", 

55 "stresses", 

56 "dipole", 

57 "magmom", 

58 ], 

59 ) 

60 

61 self.outputname = "aims.out" 

62 

63 def update_parameters(self, properties, parameters): 

64 """Check and update the parameters to match the desired calculation 

65 

66 Parameters 

67 ---------- 

68 properties: list of str 

69 The list of properties to calculate 

70 parameters: dict 

71 The parameters used to perform the calculation. 

72 

73 Returns 

74 ------- 

75 dict 

76 The updated parameters object 

77 """ 

78 parameters = dict(parameters) 

79 property_flags = { 

80 "forces": "compute_forces", 

81 "stress": "compute_analytical_stress", 

82 "stresses": "compute_heat_flux", 

83 } 

84 # Ensure FHI-aims will calculate all desired properties 

85 for property in properties: 

86 aims_name = property_flags.get(property, None) 

87 if aims_name is not None: 

88 parameters[aims_name] = True 

89 

90 if "dipole" in properties: 

91 if "output" in parameters and "dipole" not in parameters["output"]: 

92 parameters["output"] = list(parameters["output"]) 

93 parameters["output"].append("dipole") 

94 elif "output" not in parameters: 

95 parameters["output"] = ["dipole"] 

96 

97 return parameters 

98 

99 def write_input(self, directory, atoms, parameters, properties): 

100 """Write the geometry.in and control.in files for the calculation 

101 

102 Parameters 

103 ---------- 

104 directory : Path 

105 The working directory to store the input files. 

106 atoms : atoms.Atoms 

107 The atoms object to perform the calculation on. 

108 parameters: dict 

109 The parameters used to perform the calculation. 

110 properties: list of str 

111 The list of properties to calculate 

112 """ 

113 parameters = self.update_parameters(properties, parameters) 

114 

115 ghosts = parameters.pop("ghosts", None) 

116 geo_constrain = parameters.pop("geo_constrain", None) 

117 scaled = parameters.pop("scaled", None) 

118 write_velocities = parameters.pop("write_velocities", None) 

119 

120 if scaled is None: 

121 scaled = np.all(atoms.pbc) 

122 if write_velocities is None: 

123 write_velocities = atoms.has("momenta") 

124 

125 if geo_constrain is None: 

126 geo_constrain = scaled and "relax_geometry" in parameters 

127 

128 have_lattice_vectors = atoms.pbc.any() 

129 have_k_grid = ("k_grid" in parameters or "kpts" in parameters 

130 or "k_grid_density" in parameters) 

131 if have_lattice_vectors and not have_k_grid: 

132 raise RuntimeError("Found lattice vectors but no k-grid!") 

133 if not have_lattice_vectors and have_k_grid: 

134 raise RuntimeError("Found k-grid but no lattice vectors!") 

135 

136 geometry_in = directory / "geometry.in" 

137 

138 write_aims( 

139 geometry_in, 

140 atoms, 

141 scaled, 

142 geo_constrain, 

143 write_velocities=write_velocities, 

144 ghosts=ghosts, 

145 ) 

146 

147 control = directory / "control.in" 

148 write_control(control, atoms, parameters) 

149 

150 def execute(self, directory, profile): 

151 profile.run(directory, self.outputname) 

152 

153 def read_results(self, directory): 

154 from ase.io.aims import read_aims_results 

155 

156 dst = directory / self.outputname 

157 return read_aims_results(dst, index=-1) 

158 

159 

160class Aims(GenericFileIOCalculator): 

161 def __init__(self, profile=None, directory='.', **kwargs): 

162 """Construct the FHI-aims calculator. 

163 

164 The keyword arguments (kwargs) can be one of the ASE standard 

165 keywords: 'xc', 'kpts' and 'smearing' or any of FHI-aims' 

166 native keywords. 

167 

168 

169 Arguments: 

170 

171 cubes: AimsCube object 

172 Cube file specification. 

173 

174 tier: int or array of ints 

175 Set basis set tier for all atomic species. 

176 

177 plus_u : dict 

178 For DFT+U. Adds a +U term to one specific shell of the species. 

179 

180 kwargs : dict 

181 Any of the base class arguments. 

182 

183 """ 

184 

185 if profile is None: 

186 profile = AimsProfile( 

187 kwargs.pop( 

188 "run_command", 

189 os.getenv("ASE_AIMS_COMMAND", "aims.x") 

190 ) 

191 ) 

192 

193 super().__init__(template=AimsTemplate(), 

194 profile=profile, 

195 parameters=kwargs, 

196 directory=directory) 

197 

198 

199class AimsCube: 

200 "Object to ensure the output of cube files, can be attached to Aims object" 

201 

202 def __init__( 

203 self, 

204 origin=(0, 0, 0), 

205 edges=[(0.1, 0.0, 0.0), (0.0, 0.1, 0.0), (0.0, 0.0, 0.1)], 

206 points=(50, 50, 50), 

207 plots=tuple(), 

208 ): 

209 """parameters: 

210 

211 origin, edges, points: 

212 Same as in the FHI-aims output 

213 plots: 

214 what to print, same names as in FHI-aims""" 

215 

216 self.name = "AimsCube" 

217 self.origin = origin 

218 self.edges = edges 

219 self.points = points 

220 self.plots = plots 

221 

222 def ncubes(self): 

223 """returns the number of cube files to output """ 

224 return len(self.plots) 

225 

226 def move_to_base_name(self, basename): 

227 """when output tracking is on or the base namem is not standard, 

228 this routine will rename add the base to the cube file output for 

229 easier tracking""" 

230 for plot in self.plots: 

231 found = False 

232 cube = plot.split() 

233 if ( 

234 cube[0] == "total_density" 

235 or cube[0] == "spin_density" 

236 or cube[0] == "delta_density" 

237 ): 

238 found = True 

239 old_name = cube[0] + ".cube" 

240 new_name = basename + "." + old_name 

241 if cube[0] == "eigenstate" or cube[0] == "eigenstate_density": 

242 found = True 

243 state = int(cube[1]) 

244 s_state = cube[1] 

245 for i in [10, 100, 1000, 10000]: 

246 if state < i: 

247 s_state = "0" + s_state 

248 old_name = cube[0] + "_" + s_state + "_spin_1.cube" 

249 new_name = basename + "." + old_name 

250 if found: 

251 # XXX Should not use platform dependent commands! 

252 os.system("mv " + old_name + " " + new_name) 

253 

254 def add_plot(self, name): 

255 """ in case you forgot one ... """ 

256 self.plots += [name] 

257 

258 def write(self, file): 

259 """ write the necessary output to the already opened control.in """ 

260 file.write("output cube " + self.plots[0] + "\n") 

261 file.write(" cube origin ") 

262 for ival in self.origin: 

263 file.write(str(ival) + " ") 

264 file.write("\n") 

265 for i in range(3): 

266 file.write(" cube edge " + str(self.points[i]) + " ") 

267 for ival in self.edges[i]: 

268 file.write(str(ival) + " ") 

269 file.write("\n") 

270 if self.ncubes() > 1: 

271 for i in range(self.ncubes() - 1): 

272 file.write("output cube " + self.plots[i + 1] + "\n")