1"""This module defines an ASE interface to FLAPW code FLEUR. 

2 

3http://www.flapw.de 

4""" 

5 

6import os 

7 

8from subprocess import Popen, PIPE 

9 

10import re 

11 

12import numpy as np 

13 

14from ase.units import Hartree, Bohr 

15from ase.calculators.calculator import PropertyNotImplementedError 

16 

17 

18class FLEUR: 

19 """Class for doing FLEUR calculations. 

20 

21 In order to use fleur one has to define the following environment 

22 variables: 

23 

24 FLEUR_INPGEN path to the input generator (inpgen.x) of fleur 

25 

26 FLEUR path to the fleur executable. Note that fleur uses different 

27 executable for real and complex cases (systems with/without inversion 

28 symmetry), so FLEUR must point to the correct executable. 

29 

30 The initialize_density step can be performed in parallel 

31 only if run on one compute node. FLEUR_SERIAL is used for this step. 

32 

33 It is probable that user needs to tune manually the input file before 

34 the actual calculation, so in addition to the standard 

35 get_potential_energy function this class defines the following utility 

36 functions: 

37 

38 write_inp 

39 generate the input file *inp* 

40 initialize_density 

41 creates the initial density after possible manual edits of *inp* 

42 calculate 

43 convergence the total energy. With fleur, one specifies always 

44 only the number of SCF-iterations so this function launches 

45 the executable several times and monitors the convergence. 

46 relax 

47 Uses fleur's internal algorithm for structure 

48 optimization. Requires that the proper optimization parameters 

49 (atoms to optimize etc.) are specified by hand in *inp* 

50 

51 """ 

52 def __init__(self, xc='LDA', kpts=None, nbands=None, convergence=None, 

53 width=None, kmax=None, mixer=None, maxiter=None, 

54 maxrelax=20, workdir=None, equivatoms=True, rmt=None, 

55 lenergy=None): 

56 """Construct FLEUR-calculator object. 

57 

58 Parameters 

59 ========== 

60 xc: str 

61 Exchange-correlation functional. Must be one of LDA, PBE, 

62 RPBE. 

63 kpts: list of three int 

64 Monkhost-Pack sampling. 

65 nbands: int 

66 Number of bands. (not used at the moment) 

67 convergence: dictionary 

68 Convergence parameters (currently only energy in eV) 

69 {'energy' : float} 

70 width: float 

71 Fermi-distribution width in eV. 

72 kmax: float 

73 Plane wave cutoff in a.u. If kmax is set then: 

74 gmax = 3.0 * kmax 

75 gmaxxc = int(2.5 * kmax * 10)/10. (from set_inp.f) 

76 mixer: dictionary 

77 Mixing parameters imix, alpha, spinf 

78 {'imix' : int, 'alpha' : float, 'spinf' : float} 

79 maxiter: int 

80 Maximum number of SCF iterations (name in the code: itmax) 

81 maxrelax: int 

82 Maximum number of relaxation steps 

83 workdir: str 

84 Working directory for the calculation 

85 equivatoms: bool 

86 If False: generate inequivalent atoms (default is True). 

87 Setting to False allows one for example to calculate spin-polarized dimers. 

88 See http://www.flapw.de/pm/index.php?n=User-Documentation.InputFileForTheInputGenerator. 

89 rmt: dictionary 

90 rmt values in Angstrom., e.g: {'O': 1.1 * Bohr, 'N': -0.1} 

91 Negative number with respect to the rmt set by FLEUR. 

92 lenergy: float 

93 Lower energy in eV. Default -1.8 * Hartree. 

94 """ 

95 

96 self.xc = xc 

97 self.kpts = kpts 

98 self.nbands = nbands 

99 self.width = width 

100 self.kmax = kmax 

101 self.itmax_step_default = 9 # SCF steps per run (default) 

102 self.itmax_step = 5 # SCF steps per run 

103 assert self.itmax_step_default <= 9 

104 assert self.itmax_step <= self.itmax_step_default 

105 self.itmax_default = 40 

106 if maxiter is None: 

107 self.itmax = self.itmax_default 

108 else: 

109 self.itmax = maxiter 

110 self.maxrelax = maxrelax 

111 self.mixer = mixer 

112 

113 if convergence: 

114 self.convergence = convergence 

115 self.convergence['energy'] /= Hartree 

116 else: 

117 self.convergence = {'energy': 0.0001} 

118 

119 self.start_dir = None 

120 self.workdir = workdir 

121 if self.workdir: 

122 self.start_dir = os.getcwd() 

123 if not os.path.isdir(workdir): 

124 os.mkdir(workdir) 

125 else: 

126 self.workdir = '.' 

127 self.start_dir = '.' 

128 

129 self.equivatoms = equivatoms 

130 

131 self.rmt = rmt 

132 self.lenergy = lenergy 

133 

134 self.converged = False 

135 

136 def run_executable(self, mode='fleur', executable='FLEUR'): 

137 

138 assert executable in ['FLEUR', 'FLEUR_SERIAL'] 

139 

140 executable_use = executable 

141 if executable == 'FLEUR_SERIAL' and not os.environ.get(executable, ''): 

142 executable_use = 'FLEUR' # use FLEUR if FLEUR_SERIAL not set 

143 try: 

144 code_exe = os.environ[executable_use] 

145 except KeyError: 

146 raise RuntimeError('Please set ' + executable_use) 

147 p = Popen(code_exe, shell=True, stdin=PIPE, stdout=PIPE, 

148 stderr=PIPE) 

149 stat = p.wait() 

150 out = p.stdout.read() 

151 err = p.stderr.read() 

152 print(mode, ': stat= ', stat, ' out= ', out, ' err=', err) 

153 # special handling of exit status from density generation and regular fleur.x 

154 if mode in ['density']: 

155 if '!' in err: 

156 os.chdir(self.start_dir) 

157 raise RuntimeError(executable_use + ' exited with a code %s' % err) 

158 else: 

159 if stat != 0: 

160 os.chdir(self.start_dir) 

161 raise RuntimeError(executable_use + ' exited with a code %d' % stat) 

162 

163 def update(self, atoms): 

164 """Update a FLEUR calculation.""" 

165 

166 if (not self.converged or 

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

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

169 self.initialize(atoms) 

170 self.calculate(atoms) 

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

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

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

174 self.converged = False 

175 self.initialize(atoms) 

176 self.calculate(atoms) 

177 

178 def initialize(self, atoms): 

179 """Create an input file inp and generate starting density.""" 

180 

181 self.converged = False 

182 self.initialize_inp(atoms) 

183 self.initialize_density(atoms) 

184 

185 def initialize_inp(self, atoms): 

186 """Create a inp file""" 

187 os.chdir(self.workdir) 

188 

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

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

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

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

193 

194 # create the input 

195 self.write_inp(atoms) 

196 

197 os.chdir(self.start_dir) 

198 

199 def initialize_density(self, atoms): 

200 """Creates a new starting density.""" 

201 

202 os.chdir(self.workdir) 

203 # remove possible conflicting files 

204 files2remove = ['cdn1', 'fl7para', 'stars', 'wkf2', 'enpara', 

205 'kpts', 'broyd', 'broyd.7', 'tmat', 'tmas'] 

206 if 0: 

207 # avoid STOP bzone3 error by keeping the kpts file 

208 files2remove.remove('kpts') 

209 

210 for f in files2remove: 

211 if os.path.isfile(f): 

212 os.remove(f) 

213 

214 # generate the starting density 

215 os.system("sed -i -e 's/strho=./strho=T/' inp") 

216 self.run_executable(mode='density', executable='FLEUR_SERIAL') 

217 os.system("sed -i -e 's/strho=./strho=F/' inp") 

218 

219 os.chdir(self.start_dir) 

220 # generate spin-polarized density 

221 # http://www.flapw.de/pm/index.php?n=User-Documentation.Magnetism 

222 if atoms.get_initial_magnetic_moments().sum() > 0.0: 

223 os.chdir(self.workdir) 

224 # generate cdnc file (1 SCF step: swsp=F - non-magnetic) 

225 os.system("sed -i -e 's/itmax=.*,maxiter/itmax= 1,maxiter/' inp") 

226 self.run_executable(mode='cdnc', executable='FLEUR') 

227 sedline = "'s/itmax=.*,maxiter/itmax= '" 

228 sedline += str(self.itmax_step_default) + "',maxiter/'" 

229 os.system("sed -i -e " + sedline + " inp") 

230 # generate spin polarized density (swsp=T) 

231 os.system("sed -i -e 's/swsp=./swsp=T/' inp") 

232 self.run_executable(mode='swsp', executable='FLEUR_SERIAL') 

233 # restore swsp=F 

234 os.system("sed -i -e 's/swsp=./swsp=F/' inp") 

235 os.chdir(self.start_dir) 

236 

237 def get_potential_energy(self, atoms, force_consistent=False): 

238 self.update(atoms) 

239 

240 if force_consistent: 

241 return self.efree * Hartree 

242 else: 

243 # Energy extrapolated to zero Kelvin: 

244 return (self.etotal + self.efree) / 2 * Hartree 

245 

246 def get_number_of_iterations(self, atoms): 

247 self.update(atoms) 

248 return self.niter 

249 

250 def get_forces(self, atoms): 

251 self.update(atoms) 

252 # electronic structure is converged, so let's calculate forces: 

253 # TODO 

254 return np.array((0.0, 0.0, 0.0)) 

255 

256 def get_stress(self, atoms): 

257 raise PropertyNotImplementedError 

258 

259 def get_dipole_moment(self, atoms): 

260 """Returns total dipole moment of the system.""" 

261 raise PropertyNotImplementedError 

262 

263 def calculate(self, atoms): 

264 """Converge a FLEUR calculation to self-consistency. 

265 

266 Input files should be generated before calling this function 

267 FLEUR performs always fixed number of SCF steps. This function 

268 reduces the number of iterations gradually, however, a minimum 

269 of five SCF steps is always performed. 

270 """ 

271 

272 os.chdir(self.workdir) 

273 

274 self.niter = 0 

275 out = '' 

276 err = '' 

277 while not self.converged: 

278 if self.niter > self.itmax: 

279 os.chdir(self.start_dir) 

280 raise RuntimeError('FLEUR failed to convergence in %d iterations' % self.itmax) 

281 

282 self.run_executable(mode='fleur', executable='FLEUR') 

283 

284 # catenate new output with the old one 

285 os.system('cat out >> out.old') 

286 self.read() 

287 self.check_convergence() 

288 

289 if os.path.exists('out.old'): 

290 os.rename('out.old', 'out') 

291 # After convergence clean up broyd* files 

292 os.system('rm -f broyd*') 

293 os.chdir(self.start_dir) 

294 return out, err 

295 

296 def relax(self, atoms): 

297 """Currently, user has to manually define relaxation parameters 

298 (atoms to relax, relaxation directions, etc.) in inp file 

299 before calling this function.""" 

300 

301 nrelax = 0 

302 relaxed = False 

303 while not relaxed: 

304 # Calculate electronic structure 

305 self.calculate(atoms) 

306 # Calculate the Pulay forces 

307 os.system("sed -i -e 's/l_f=./l_f=T/' inp") 

308 while True: 

309 self.converged = False 

310 out, err = self.calculate(atoms) 

311 if 'GEO new' in err: 

312 os.chdir(self.workdir) 

313 os.rename('inp_new', 'inp') 

314 os.chdir(self.start_dir) 

315 break 

316 if 'GEO: Des woas' in err: 

317 relaxed = True 

318 break 

319 nrelax += 1 

320 # save the out and cdn1 files 

321 os.system('cp out out_%d' % nrelax) 

322 os.system('cp cdn1 cdn1_%d' % nrelax) 

323 if nrelax > self.maxrelax: 

324 os.chdir(self.start_dir) 

325 raise RuntimeError('Failed to relax in %d iterations' % self.maxrelax) 

326 self.converged = False 

327 

328 def write_inp(self, atoms): 

329 """Write the *inp* input file of FLEUR. 

330 

331 First, the information from Atoms is written to the simple input 

332 file and the actual input file *inp* is then generated with the 

333 FLEUR input generator. The location of input generator is specified 

334 in the environment variable FLEUR_INPGEN. 

335 

336 Finally, the *inp* file is modified according to the arguments of 

337 the FLEUR calculator object. 

338 """ 

339 

340 with open('inp_simple', 'w') as fh: 

341 self._write_inp(atoms, fh) 

342 

343 def _write_inp(self, atoms, fh): 

344 fh.write('FLEUR input generated with ASE\n') 

345 fh.write('\n') 

346 

347 if atoms.pbc[2]: 

348 film = 'f' 

349 else: 

350 film = 't' 

351 fh.write('&input film=%s /' % film) 

352 fh.write('\n') 

353 

354 for vec in atoms.get_cell(): 

355 fh.write(' ') 

356 for el in vec: 

357 fh.write(' %21.16f' % (el/Bohr)) 

358 fh.write('\n') 

359 fh.write(' %21.16f\n' % 1.0) 

360 fh.write(' %21.16f %21.16f %21.16f\n' % (1.0, 1.0, 1.0)) 

361 fh.write('\n') 

362 

363 natoms = len(atoms) 

364 fh.write(' %6d\n' % natoms) 

365 positions = atoms.get_scaled_positions() 

366 if not atoms.pbc[2]: 

367 # in film calculations z position has to be in absolute 

368 # coordinates and symmetrical 

369 cart_pos = atoms.get_positions() 

370 cart_pos[:, 2] -= atoms.get_cell()[2, 2]/2.0 

371 positions[:, 2] = cart_pos[:, 2] / Bohr 

372 atomic_numbers = atoms.get_atomic_numbers() 

373 for n, (Z, pos) in enumerate(zip(atomic_numbers, positions)): 

374 if self.equivatoms: 

375 fh.write('%3d' % Z) 

376 else: 

377 # generate inequivalent atoms, by using non-integer Z 

378 # (only the integer part will be used as Z of the atom) 

379 # see http://www.flapw.de/pm/index.php?n=User-Documentation.InputFileForTheInputGenerator 

380 fh.write('%3d.%04d' % (Z, n)) # MDTMP don't think one can calculate more that 10**4 atoms 

381 for el in pos: 

382 fh.write(' %21.16f' % el) 

383 fh.write('\n') 

384 

385 # avoid "STOP read_record: ERROR reading input" 

386 fh.write('&end /') 

387 

388 try: 

389 inpgen = os.environ['FLEUR_INPGEN'] 

390 except KeyError: 

391 raise RuntimeError('Please set FLEUR_INPGEN') 

392 

393 # rename the previous inp if it exists 

394 if os.path.isfile('inp'): 

395 os.rename('inp', 'inp.bak') 

396 os.system('%s -old < inp_simple' % inpgen) 

397 

398 # read the whole inp-file for possible modifications 

399 with open('inp', 'r') as fh: 

400 lines = fh.readlines() 

401 

402 window_ln = -1 

403 for ln, line in enumerate(lines): 

404 # XC potential 

405 if line.startswith('pbe'): 

406 if self.xc == 'PBE': 

407 pass 

408 elif self.xc == 'RPBE': 

409 lines[ln] = 'rpbe non-relativi\n' 

410 elif self.xc == 'LDA': 

411 lines[ln] = 'mjw non-relativic\n' 

412 del lines[ln+1] 

413 else: 

414 raise RuntimeError('XC-functional %s is not supported' % self.xc) 

415 if line.startswith('Window'): 

416 # few things are set around this line 

417 window_ln = ln 

418 # kmax 

419 if self.kmax and ln == window_ln: 

420 line = '%10.5f\n' % self.kmax 

421 lines[ln+2] = line 

422 # lower energy 

423 if self.lenergy is not None and ln == window_ln: 

424 l0 = lines[ln+1].split()[0] 

425 l = lines[ln+1].replace(l0, '%8.5f' % (self.lenergy / Hartree)) 

426 lines[ln+1] = l 

427 # gmax cutoff for PW-expansion of potential & density ( > 2*kmax) 

428 # gmaxxc cutoff for PW-expansion of XC-potential ( > 2*kmax, < gmax) 

429 if self.kmax and line.startswith('vchk'): 

430 gmax = 3. * self.kmax 

431 line = ' %10.6f %10.6f\n' % (gmax, int(2.5 * self.kmax * 10)/10.) 

432 lines[ln-1] = line 

433 # Fermi width 

434 if self.width and line.startswith('gauss'): 

435 line = 'gauss=F %7.5ftria=F\n' % (self.width / Hartree) 

436 lines[ln] = line 

437 # kpts 

438 if self.kpts and line.startswith('nkpt'): 

439 line = 'nkpt= nx=%2d,ny=%2d,nz=%2d\n' % (self.kpts[0], 

440 self.kpts[1], 

441 self.kpts[2]) 

442 lines[ln] = line 

443 # itmax 

444 if self.itmax < self.itmax_step_default and line.startswith('itmax'): 

445 # decrease number of SCF steps; increasing is done by 'while not self.converged:' 

446 lsplit = line.split(',') 

447 if lsplit[0].find('itmax') != -1: 

448 lsplit[0] = 'itmax=' + ('%2d' % self.itmax) 

449 lines[ln] = ",".join(lsplit) 

450 # Mixing 

451 if self.mixer and line.startswith('itmax'): 

452 imix = self.mixer['imix'] 

453 alpha = self.mixer['alpha'] 

454 spinf = self.mixer['spinf'] 

455 line_end = 'imix=%2d,alpha=%6.2f,spinf=%6.2f\n' % (imix, 

456 alpha, 

457 spinf) 

458 line = line[:21] + line_end 

459 lines[ln] = line 

460 # jspins and swsp 

461 if atoms.get_initial_magnetic_moments().sum() > 0.0: 

462 assert not self.equivatoms, 'equivatoms currently not allowed in magnetic systems' 

463 if line.find('jspins=1') != -1: 

464 lines[ln] = line.replace('jspins=1', 'jspins=2') 

465 if line.startswith('swsp=F'): 

466 # setting initial magnetic moments for all atom types 

467 lines[ln] = 'swsp=F' 

468 for m in atoms.get_initial_magnetic_moments(): 

469 lines[ln] += (' %5.2f' % m) 

470 lines[ln] += '\n' 

471 # inpgen produces incorrect symbol 'J' for Iodine 

472 if line.startswith(' J 53'): 

473 lines[ln] = lines[ln].replace(' J 53', ' I 53') 

474 # rmt 

475 if self.rmt is not None: 

476 for s in list(set(atoms.get_chemical_symbols())): # unique 

477 if s in self.rmt: 

478 # set the requested rmt 

479 for ln, line in enumerate(lines): 

480 ls = line.split() 

481 if len(ls) == 7 and ls[0].strip() == s: 

482 rorig = ls[5].strip() 

483 if self.rmt[s] < 0.0: 

484 r = float(rorig) + self.rmt[s] / Bohr 

485 else: 

486 r = self.rmt[s] / Bohr 

487 print(s, rorig, r) 

488 lines[ln] = lines[ln].replace(rorig, ("%.6f" % r)) 

489 

490 # write everything back to inp 

491 with open('inp', 'w') as fh: 

492 for line in lines: 

493 fh.write(line) 

494 

495 def read(self): 

496 """Read results from FLEUR's text-output file `out`.""" 

497 

498 with open('out', 'r') as fd: 

499 lines = fd.readlines() 

500 

501 # total energies 

502 self.total_energies = [] 

503 pat = re.compile(r'(.*total energy=)(\s)*([-0-9.]*)') 

504 for line in lines: 

505 m = pat.match(line) 

506 if m: 

507 self.total_energies.append(float(m.group(3))) 

508 self.etotal = self.total_energies[-1] 

509 

510 # free_energies 

511 self.free_energies = [] 

512 pat = re.compile(r'(.*free energy=)(\s)*([-0-9.]*)') 

513 for line in lines: 

514 m = pat.match(line) 

515 if m: 

516 self.free_energies.append(float(m.group(3))) 

517 self.efree = self.free_energies[-1] 

518 

519 # TODO forces, charge density difference... 

520 

521 def check_convergence(self): 

522 """Check the convergence of calculation""" 

523 energy_error = np.ptp(self.total_energies[-3:]) 

524 self.converged = energy_error < self.convergence['energy'] 

525 

526 # TODO check charge convergence 

527 

528 # reduce the itmax in inp 

529 with open('inp', 'r') as fh: 

530 lines = fh.readlines() 

531 pat = re.compile('(itmax=)([ 0-9]*)') 

532 with open('inp', 'w') as fh: 

533 for line in lines: 

534 m = pat.match(line) 

535 if m: 

536 itmax = int(m.group(2)) 

537 self.niter += itmax 

538 itmax_new = itmax // 2 

539 itmax = max(self.itmax_step, itmax_new) 

540 line = 'itmax=%2d' % itmax + line[8:] 

541 fh.write(line)