Coverage for /builds/debichem-team/python-ase/ase/calculators/aims.py: 42.28%

123 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-03-06 04:00 +0000

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.calculators.genericfileio import ( 

16 BaseProfile, 

17 CalculatorTemplate, 

18 GenericFileIOCalculator, 

19 read_stdout, 

20) 

21from ase.io.aims import write_aims, write_control 

22 

23 

24def get_aims_version(string): 

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

26 return match.group(1) 

27 

28 

29class AimsProfile(BaseProfile): 

30 configvars = {'default_species_directory'} 

31 

32 def __init__(self, command, default_species_directory=None, **kwargs): 

33 super().__init__(command, **kwargs) 

34 self.default_species_directory = default_species_directory 

35 

36 def get_calculator_command(self, inputfile): 

37 return [] 

38 

39 def version(self): 

40 return get_aims_version(read_stdout(self._split_command)) 

41 

42 

43class AimsTemplate(CalculatorTemplate): 

44 _label = 'aims' 

45 

46 def __init__(self): 

47 super().__init__( 

48 'aims', 

49 [ 

50 'energy', 

51 'free_energy', 

52 'forces', 

53 'stress', 

54 'stresses', 

55 'dipole', 

56 'magmom', 

57 ], 

58 ) 

59 

60 self.outputname = f'{self._label}.out' 

61 self.errorname = f'{self._label}.err' 

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, profile, 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 = ( 

130 'k_grid' in parameters 

131 or 'kpts' in parameters 

132 or 'k_grid_density' in parameters 

133 ) 

134 if have_lattice_vectors and not have_k_grid: 

135 raise RuntimeError('Found lattice vectors but no k-grid!') 

136 if not have_lattice_vectors and have_k_grid: 

137 raise RuntimeError('Found k-grid but no lattice vectors!') 

138 

139 geometry_in = directory / 'geometry.in' 

140 

141 write_aims( 

142 geometry_in, 

143 atoms, 

144 scaled, 

145 geo_constrain, 

146 write_velocities=write_velocities, 

147 ghosts=ghosts, 

148 ) 

149 

150 control = directory / 'control.in' 

151 

152 if ( 

153 'species_dir' not in parameters 

154 and profile.default_species_directory is not None 

155 ): 

156 parameters['species_dir'] = profile.default_species_directory 

157 

158 write_control(control, atoms, parameters) 

159 

160 def execute(self, directory, profile): 

161 profile.run(directory, None, self.outputname, 

162 errorfile=self.errorname) 

163 

164 def read_results(self, directory): 

165 from ase.io.aims import read_aims_results 

166 

167 dst = directory / self.outputname 

168 return read_aims_results(dst, index=-1) 

169 

170 def load_profile(self, cfg, **kwargs): 

171 return AimsProfile.from_config(cfg, self.name, **kwargs) 

172 

173 def socketio_argv(self, profile, unixsocket, port): 

174 return [profile.command] 

175 

176 def socketio_parameters(self, unixsocket, port): 

177 if port: 

178 use_pimd_wrapper = ('localhost', port) 

179 else: 

180 # (INET port number should be unused.) 

181 use_pimd_wrapper = (f'UNIX:{unixsocket}', 31415) 

182 

183 return dict(use_pimd_wrapper=use_pimd_wrapper, compute_forces=True) 

184 

185 

186class Aims(GenericFileIOCalculator): 

187 def __init__( 

188 self, 

189 profile=None, 

190 directory='.', 

191 **kwargs, 

192 ): 

193 """Construct the FHI-aims calculator. 

194 

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

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

197 native keywords. 

198 

199 

200 Arguments: 

201 

202 cubes: AimsCube object 

203 Cube file specification. 

204 

205 tier: int or array of ints 

206 Set basis set tier for all atomic species. 

207 

208 plus_u : dict 

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

210 

211 kwargs : dict 

212 Any of the base class arguments. 

213 

214 """ 

215 

216 super().__init__( 

217 template=AimsTemplate(), 

218 profile=profile, 

219 parameters=kwargs, 

220 directory=directory, 

221 ) 

222 

223 

224class AimsCube: 

225 'Object to ensure the output of cube files, can be attached to Aims object' 

226 

227 def __init__( 

228 self, 

229 origin=(0, 0, 0), 

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

231 points=(50, 50, 50), 

232 plots=(), 

233 ): 

234 """parameters: 

235 

236 origin, edges, points: 

237 Same as in the FHI-aims output 

238 plots: 

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

240 

241 self.name = 'AimsCube' 

242 self.origin = origin 

243 self.edges = edges 

244 self.points = points 

245 self.plots = plots 

246 

247 def ncubes(self): 

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

249 return len(self.plots) 

250 

251 def move_to_base_name(self, basename): 

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

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

254 easier tracking""" 

255 for plot in self.plots: 

256 found = False 

257 cube = plot.split() 

258 if ( 

259 cube[0] == 'total_density' 

260 or cube[0] == 'spin_density' 

261 or cube[0] == 'delta_density' 

262 ): 

263 found = True 

264 old_name = cube[0] + '.cube' 

265 new_name = basename + '.' + old_name 

266 if cube[0] == 'eigenstate' or cube[0] == 'eigenstate_density': 

267 found = True 

268 state = int(cube[1]) 

269 s_state = cube[1] 

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

271 if state < i: 

272 s_state = '0' + s_state 

273 old_name = cube[0] + '_' + s_state + '_spin_1.cube' 

274 new_name = basename + '.' + old_name 

275 if found: 

276 # XXX Should not use platform dependent commands! 

277 os.system('mv ' + old_name + ' ' + new_name) 

278 

279 def add_plot(self, name): 

280 """in case you forgot one ...""" 

281 self.plots += [name] 

282 

283 def write(self, file): 

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

285 file.write('output cube ' + self.plots[0] + '\n') 

286 file.write(' cube origin ') 

287 for ival in self.origin: 

288 file.write(str(ival) + ' ') 

289 file.write('\n') 

290 for i in range(3): 

291 file.write(' cube edge ' + str(self.points[i]) + ' ') 

292 for ival in self.edges[i]: 

293 file.write(str(ival) + ' ') 

294 file.write('\n') 

295 if self.ncubes() > 1: 

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

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