Coverage for /builds/hweiske/ase/ase/calculators/aims.py: 41.46%

123 statements  

« prev     ^ index     » next       coverage.py v7.2.7, created at 2024-04-22 11:22 +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 (BaseProfile, CalculatorTemplate, 

16 GenericFileIOCalculator, 

17 read_stdout) 

18from ase.io.aims import write_aims, write_control 

19 

20 

21def get_aims_version(string): 

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

23 return match.group(1) 

24 

25 

26class AimsProfile(BaseProfile): 

27 def __init__(self, binary, default_species_directory=None, **kwargs): 

28 super().__init__(**kwargs) 

29 self.binary = binary 

30 self.default_species_directory = default_species_directory 

31 

32 def get_calculator_command(self, inputfile): 

33 return [self.binary] 

34 

35 def version(self): 

36 return get_aims_version(read_stdout(self.binary)) 

37 

38 

39class AimsTemplate(CalculatorTemplate): 

40 _label = 'aims' 

41 

42 def __init__(self): 

43 super().__init__( 

44 'aims', 

45 [ 

46 'energy', 

47 'free_energy', 

48 'forces', 

49 'stress', 

50 'stresses', 

51 'dipole', 

52 'magmom', 

53 ], 

54 ) 

55 

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

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

58 

59 def update_parameters(self, properties, parameters): 

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

61 

62 Parameters 

63 ---------- 

64 properties: list of str 

65 The list of properties to calculate 

66 parameters: dict 

67 The parameters used to perform the calculation. 

68 

69 Returns 

70 ------- 

71 dict 

72 The updated parameters object 

73 """ 

74 parameters = dict(parameters) 

75 property_flags = { 

76 'forces': 'compute_forces', 

77 'stress': 'compute_analytical_stress', 

78 'stresses': 'compute_heat_flux', 

79 } 

80 # Ensure FHI-aims will calculate all desired properties 

81 for property in properties: 

82 aims_name = property_flags.get(property, None) 

83 if aims_name is not None: 

84 parameters[aims_name] = True 

85 

86 if 'dipole' in properties: 

87 if 'output' in parameters and 'dipole' not in parameters['output']: 

88 parameters['output'] = list(parameters['output']) 

89 parameters['output'].append('dipole') 

90 elif 'output' not in parameters: 

91 parameters['output'] = ['dipole'] 

92 

93 return parameters 

94 

95 def write_input(self, profile, directory, atoms, parameters, properties): 

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

97 

98 Parameters 

99 ---------- 

100 directory : Path 

101 The working directory to store the input files. 

102 atoms : atoms.Atoms 

103 The atoms object to perform the calculation on. 

104 parameters: dict 

105 The parameters used to perform the calculation. 

106 properties: list of str 

107 The list of properties to calculate 

108 """ 

109 parameters = self.update_parameters(properties, parameters) 

110 

111 ghosts = parameters.pop('ghosts', None) 

112 geo_constrain = parameters.pop('geo_constrain', None) 

113 scaled = parameters.pop('scaled', None) 

114 write_velocities = parameters.pop('write_velocities', None) 

115 

116 if scaled is None: 

117 scaled = np.all(atoms.pbc) 

118 if write_velocities is None: 

119 write_velocities = atoms.has('momenta') 

120 

121 if geo_constrain is None: 

122 geo_constrain = scaled and 'relax_geometry' in parameters 

123 

124 have_lattice_vectors = atoms.pbc.any() 

125 have_k_grid = ( 

126 'k_grid' in parameters 

127 or 'kpts' in parameters 

128 or 'k_grid_density' in parameters 

129 ) 

130 if have_lattice_vectors and not have_k_grid: 

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

132 if not have_lattice_vectors and have_k_grid: 

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

134 

135 geometry_in = directory / 'geometry.in' 

136 

137 write_aims( 

138 geometry_in, 

139 atoms, 

140 scaled, 

141 geo_constrain, 

142 write_velocities=write_velocities, 

143 ghosts=ghosts, 

144 ) 

145 

146 control = directory / 'control.in' 

147 

148 if ( 

149 'species_dir' not in parameters 

150 and profile.default_species_directory is not None 

151 ): 

152 parameters['species_dir'] = profile.default_species_directory 

153 

154 write_control(control, atoms, parameters) 

155 

156 def execute(self, directory, profile): 

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

158 errorfile=self.errorname) 

159 

160 def read_results(self, directory): 

161 from ase.io.aims import read_aims_results 

162 

163 dst = directory / self.outputname 

164 return read_aims_results(dst, index=-1) 

165 

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

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

168 

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

170 return [profile.binary] 

171 

172 def socketio_parameters(self, unixsocket, port): 

173 if port: 

174 use_pimd_wrapper = ('localhost', port) 

175 else: 

176 # (INET port number should be unused.) 

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

178 

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

180 

181 

182class Aims(GenericFileIOCalculator): 

183 def __init__( 

184 self, 

185 profile=None, 

186 directory='.', 

187 parallel_info=None, 

188 parallel=True, 

189 **kwargs, 

190 ): 

191 """Construct the FHI-aims calculator. 

192 

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

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

195 native keywords. 

196 

197 

198 Arguments: 

199 

200 cubes: AimsCube object 

201 Cube file specification. 

202 

203 tier: int or array of ints 

204 Set basis set tier for all atomic species. 

205 

206 plus_u : dict 

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

208 

209 kwargs : dict 

210 Any of the base class arguments. 

211 

212 """ 

213 

214 super().__init__( 

215 template=AimsTemplate(), 

216 profile=profile, 

217 parameters=kwargs, 

218 parallel_info=parallel_info, 

219 parallel=parallel, 

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')