1# type: ignore 

2""" 

3This module defines an ASE interface to Turbomole: http://www.turbomole.com/ 

4 

5QMMM functionality provided by Markus Kaukonen <markus.kaukonen@iki.fi>. 

6 

7Please read the license file (../../LICENSE) 

8 

9Contact: Ivan Kondov <ivan.kondov@kit.edu> 

10""" 

11import os 

12import re 

13import warnings 

14from math import log10, floor 

15import numpy as np 

16from ase.units import Ha, Bohr 

17from ase.io import read, write 

18from ase.calculators.calculator import FileIOCalculator 

19from ase.calculators.calculator import PropertyNotImplementedError, ReadError 

20from ase.calculators.turbomole.executor import execute, get_output_filename 

21from ase.calculators.turbomole.writer import add_data_group, delete_data_group 

22from ase.calculators.turbomole import reader 

23from ase.calculators.turbomole.reader import read_data_group, read_convergence 

24from ase.calculators.turbomole.parameters import TurbomoleParameters 

25 

26 

27class TurbomoleOptimizer: 

28 def __init__(self, atoms, calc): 

29 self.atoms = atoms 

30 self.calc = calc 

31 self.atoms.calc = self.calc 

32 

33 def todict(self): 

34 return {'type': 'optimization', 

35 'optimizer': 'TurbomoleOptimizer'} 

36 

37 def run(self, fmax=None, steps=None): 

38 if fmax is not None: 

39 self.calc.parameters['force convergence'] = fmax 

40 self.calc.parameters.verify() 

41 if steps is not None: 

42 self.calc.parameters['geometry optimization iterations'] = steps 

43 self.calc.parameters.verify() 

44 self.calc.calculate() 

45 self.atoms.positions[:] = self.calc.atoms.positions 

46 self.calc.parameters['task'] = 'energy' 

47 

48 

49class Turbomole(FileIOCalculator): 

50 

51 """constants""" 

52 name = 'Turbomole' 

53 

54 implemented_properties = ['energy', 'forces', 'dipole', 'free_energy', 

55 'charges'] 

56 

57 tm_files = [ 

58 'control', 'coord', 'basis', 'auxbasis', 'energy', 'gradient', 'mos', 

59 'alpha', 'beta', 'statistics', 'GEO_OPT_CONVERGED', 'GEO_OPT_FAILED', 

60 'not.converged', 'nextstep', 'hessapprox', 'job.last', 'job.start', 

61 'optinfo', 'statistics', 'converged', 'vibspectrum', 

62 'vib_normal_modes', 'hessian', 'dipgrad', 'dscf_problem', 'pc.txt', 

63 'pc_gradients.txt' 

64 ] 

65 tm_tmp_files = [ 

66 'errvec', 'fock', 'oldfock', 'dens', 'ddens', 'diff_densmat', 

67 'diff_dft_density', 'diff_dft_oper', 'diff_fockmat', 'diis_errvec', 

68 'diis_oldfock' 

69 ] 

70 

71 # initialize attributes 

72 results = {} 

73 initialized = False 

74 pc_initialized = False 

75 converged = False 

76 updated = False 

77 update_energy = None 

78 update_forces = None 

79 update_geometry = None 

80 update_hessian = None 

81 atoms = None 

82 forces = None 

83 e_total = None 

84 dipole = None 

85 charges = None 

86 version = None 

87 runtime = None 

88 datetime = None 

89 hostname = None 

90 pcpot = None 

91 

92 def __init__(self, label=None, calculate_energy='dscf', 

93 calculate_forces='grad', post_HF=False, atoms=None, 

94 restart=False, define_str=None, control_kdg=None, 

95 control_input=None, reset_tolerance=1e-2, **kwargs): 

96 

97 super().__init__(label=label) 

98 self.parameters = TurbomoleParameters() 

99 

100 self.calculate_energy = calculate_energy 

101 self.calculate_forces = calculate_forces 

102 self.post_HF = post_HF 

103 self.restart = restart 

104 self.parameters.define_str = define_str 

105 self.control_kdg = control_kdg 

106 self.control_input = control_input 

107 self.reset_tolerance = reset_tolerance 

108 

109 if self.restart: 

110 self._set_restart(kwargs) 

111 else: 

112 self.parameters.update(kwargs) 

113 self.set_parameters() 

114 self.parameters.verify() 

115 self.reset() 

116 

117 if atoms is not None: 

118 atoms.calc = self 

119 self.set_atoms(atoms) 

120 

121 def __getitem__(self, item): 

122 return getattr(self, item) 

123 

124 def _set_restart(self, params_update): 

125 """constructs atoms, parameters and results from a previous 

126 calculation""" 

127 

128 # read results, key parameters and non-key parameters 

129 self.read_restart() 

130 self.parameters.read_restart(self.atoms, self.results) 

131 

132 # update and verify parameters 

133 self.parameters.update_restart(params_update) 

134 self.set_parameters() 

135 self.parameters.verify() 

136 

137 # if a define string is specified by user then run define 

138 if self.parameters.define_str: 

139 execute(['define'], input_str=self.parameters.define_str) 

140 

141 self._set_post_define() 

142 

143 self.initialized = True 

144 # more precise convergence tests are necessary to set these flags: 

145 self.update_energy = True 

146 self.update_forces = True 

147 self.update_geometry = True 

148 self.update_hessian = True 

149 

150 def _set_post_define(self): 

151 """non-define keys, user-specified changes in the control file""" 

152 self.parameters.update_no_define_parameters() 

153 

154 # delete user-specified data groups 

155 if self.control_kdg: 

156 for dg in self.control_kdg: 

157 delete_data_group(dg) 

158 

159 # append user-defined input to control 

160 if self.control_input: 

161 for inp in self.control_input: 

162 add_data_group(inp, raw=True) 

163 

164 # add point charges if pcpot defined: 

165 if self.pcpot: 

166 self.set_point_charges() 

167 

168 def set_parameters(self): 

169 """loads the default parameters and updates with actual values""" 

170 if self.parameters.get('use resolution of identity'): 

171 self.calculate_energy = 'ridft' 

172 self.calculate_forces = 'rdgrad' 

173 

174 def reset(self): 

175 """removes all turbomole input, output and scratch files, 

176 and deletes results dict and the atoms object""" 

177 self.atoms = None 

178 self.results = {} 

179 self.results['calculation parameters'] = {} 

180 ase_files = [ 

181 f for f in os.listdir( 

182 self.directory) if f.startswith('ASE.TM.')] 

183 for f in self.tm_files + self.tm_tmp_files + ase_files: 

184 if os.path.exists(f): 

185 os.remove(f) 

186 self.initialized = False 

187 self.pc_initialized = False 

188 self.converged = False 

189 

190 def set_atoms(self, atoms): 

191 """Create the self.atoms object and writes the coord file. If 

192 self.atoms exists a check for changes and an update of the atoms 

193 is performed. Note: Only positions changes are tracked in this 

194 version. 

195 """ 

196 changes = self.check_state(atoms, tol=1e-13) 

197 if self.atoms == atoms or 'positions' not in changes: 

198 # print('two atoms obj are (almost) equal') 

199 if self.updated and os.path.isfile('coord'): 

200 self.updated = False 

201 a = read('coord').get_positions() 

202 if np.allclose(a, atoms.get_positions(), rtol=0, atol=1e-13): 

203 return 

204 else: 

205 return 

206 

207 changes = self.check_state(atoms, tol=self.reset_tolerance) 

208 if 'positions' in changes: 

209 # print(two atoms obj are different') 

210 self.reset() 

211 else: 

212 # print('two atoms obj are slightly different') 

213 if self.parameters['use redundant internals']: 

214 self.reset() 

215 

216 write('coord', atoms) 

217 self.atoms = atoms.copy() 

218 self.update_energy = True 

219 self.update_forces = True 

220 self.update_geometry = True 

221 self.update_hessian = True 

222 

223 def initialize(self): 

224 """prepare turbomole control file by running module 'define'""" 

225 if self.initialized: 

226 return 

227 self.parameters.verify() 

228 if not self.atoms: 

229 raise RuntimeError('atoms missing during initialization') 

230 if not os.path.isfile('coord'): 

231 raise IOError('file coord not found') 

232 

233 # run define 

234 define_str = self.parameters.get_define_str(len(self.atoms)) 

235 execute(['define'], input_str=define_str) 

236 

237 # process non-default initial guess 

238 iguess = self.parameters['initial guess'] 

239 if isinstance(iguess, dict) and 'use' in iguess.keys(): 

240 # "use" initial guess 

241 if self.parameters['multiplicity'] != 1 or self.parameters['uhf']: 

242 define_str = '\n\n\ny\nuse ' + iguess['use'] + '\nn\nn\nq\n' 

243 else: 

244 define_str = '\n\n\ny\nuse ' + iguess['use'] + '\nn\nq\n' 

245 execute(['define'], input_str=define_str) 

246 elif self.parameters['initial guess'] == 'hcore': 

247 # "hcore" initial guess 

248 if self.parameters['multiplicity'] != 1 or self.parameters['uhf']: 

249 delete_data_group('uhfmo_alpha') 

250 delete_data_group('uhfmo_beta') 

251 add_data_group('uhfmo_alpha', 'none file=alpha') 

252 add_data_group('uhfmo_beta', 'none file=beta') 

253 else: 

254 delete_data_group('scfmo') 

255 add_data_group('scfmo', 'none file=mos') 

256 

257 self._set_post_define() 

258 

259 self.initialized = True 

260 self.converged = False 

261 

262 def calculation_required(self, atoms, properties): 

263 if self.atoms != atoms: 

264 return True 

265 for prop in properties: 

266 if prop == 'energy' and self.e_total is None: 

267 return True 

268 elif prop == 'forces' and self.forces is None: 

269 return True 

270 return False 

271 

272 def calculate(self, atoms=None): 

273 """execute the requested job""" 

274 if atoms is None: 

275 atoms = self.atoms 

276 if self.parameters['task'] in ['energy', 'energy calculation']: 

277 self.get_potential_energy(atoms) 

278 if self.parameters['task'] in ['gradient', 'gradient calculation']: 

279 self.get_forces(atoms) 

280 if self.parameters['task'] in ['optimize', 'geometry optimization']: 

281 self.relax_geometry(atoms) 

282 if self.parameters['task'] in ['frequencies', 'normal mode analysis']: 

283 self.normal_mode_analysis(atoms) 

284 self.read_results() 

285 

286 def relax_geometry(self, atoms=None): 

287 """execute geometry optimization with script jobex""" 

288 if atoms is None: 

289 atoms = self.atoms 

290 self.set_atoms(atoms) 

291 if self.converged and not self.update_geometry: 

292 return 

293 self.initialize() 

294 jobex_command = ['jobex'] 

295 if self.parameters['use resolution of identity']: 

296 jobex_command.append('-ri') 

297 if self.parameters['force convergence']: 

298 par = self.parameters['force convergence'] 

299 conv = floor(-log10(par / Ha * Bohr)) 

300 jobex_command.extend(['-gcart', str(int(conv))]) 

301 if self.parameters['energy convergence']: 

302 par = self.parameters['energy convergence'] 

303 conv = floor(-log10(par / Ha)) 

304 jobex_command.extend(['-energy', str(int(conv))]) 

305 geom_iter = self.parameters['geometry optimization iterations'] 

306 if geom_iter is not None: 

307 assert isinstance(geom_iter, int) 

308 jobex_command.extend(['-c', str(geom_iter)]) 

309 self.converged = False 

310 execute(jobex_command) 

311 # check convergence 

312 self.converged = read_convergence(self.restart, self.parameters) 

313 if self.converged: 

314 self.update_energy = False 

315 self.update_forces = False 

316 self.update_geometry = False 

317 self.update_hessian = True 

318 # read results 

319 new_struct = read('coord') 

320 atoms.set_positions(new_struct.get_positions()) 

321 self.atoms = atoms.copy() 

322 self.read_energy() 

323 

324 def normal_mode_analysis(self, atoms=None): 

325 """execute normal mode analysis with modules aoforce or NumForce""" 

326 from ase.constraints import FixAtoms 

327 if atoms is None: 

328 atoms = self.atoms 

329 self.set_atoms(atoms) 

330 self.initialize() 

331 if self.update_energy: 

332 self.get_potential_energy(atoms) 

333 if self.update_hessian: 

334 fixatoms = [] 

335 for constr in atoms.constraints: 

336 if isinstance(constr, FixAtoms): 

337 ckwargs = constr.todict()['kwargs'] 

338 if 'indices' in ckwargs.keys(): 

339 fixatoms.extend(ckwargs['indices']) 

340 if self.parameters['numerical hessian'] is None: 

341 if len(fixatoms) > 0: 

342 define_str = '\n\ny\n' 

343 for index in fixatoms: 

344 define_str += 'm ' + str(index + 1) + ' 999.99999999\n' 

345 define_str += '*\n*\nn\nq\n' 

346 execute(['define'], input_str=define_str) 

347 dg = read_data_group('atoms') 

348 regex = r'(mass\s*=\s*)999.99999999' 

349 dg = re.sub(regex, r'\g<1>9999999999.9', dg) 

350 dg += '\n' 

351 delete_data_group('atoms') 

352 add_data_group(dg, raw=True) 

353 execute(['aoforce']) 

354 else: 

355 nforce_cmd = ['NumForce'] 

356 pdict = self.parameters['numerical hessian'] 

357 if self.parameters['use resolution of identity']: 

358 nforce_cmd.append('-ri') 

359 if len(fixatoms) > 0: 

360 nforce_cmd.extend(['-frznuclei', '-central', '-c']) 

361 if 'central' in pdict.keys(): 

362 nforce_cmd.append('-central') 

363 if 'delta' in pdict.keys(): 

364 nforce_cmd.extend(['-d', str(pdict['delta'] / Bohr)]) 

365 if self.update_forces: 

366 self.get_forces(atoms) 

367 execute(nforce_cmd) 

368 self.update_hessian = False 

369 

370 def read_restart(self): 

371 """read a previous calculation from control file""" 

372 self.atoms = read('coord') 

373 self.atoms.calc = self 

374 self.converged = read_convergence(self.restart, self.parameters) 

375 self.read_results() 

376 

377 def read_results(self): 

378 """read all results and load them in the results entity""" 

379 read_methods = [ 

380 self.read_energy, 

381 self.read_gradient, 

382 self.read_forces, 

383 self.read_basis_set, 

384 self.read_ecps, 

385 self.read_mos, 

386 self.read_occupation_numbers, 

387 self.read_dipole_moment, 

388 self.read_ssquare, 

389 self.read_hessian, 

390 self.read_vibrational_reduced_masses, 

391 self.read_normal_modes, 

392 self.read_vibrational_spectrum, 

393 self.read_charges, 

394 self.read_point_charges, 

395 self.read_run_parameters 

396 ] 

397 for method in read_methods: 

398 try: 

399 method() 

400 except ReadError as err: 

401 warnings.warn(err.args[0]) 

402 

403 def read_run_parameters(self): 

404 """read parameters set by define and not in self.parameters""" 

405 reader.read_run_parameters(self.results) 

406 

407 def read_energy(self): 

408 """Read energy from Turbomole energy file.""" 

409 reader.read_energy(self.results, self.post_HF) 

410 self.e_total = self.results['total energy'] 

411 

412 def read_forces(self): 

413 """Read forces from Turbomole gradient file.""" 

414 self.forces = reader.read_forces(self.results, len(self.atoms)) 

415 

416 def read_occupation_numbers(self): 

417 """read occupation numbers""" 

418 reader.read_occupation_numbers(self.results) 

419 

420 def read_mos(self): 

421 """read the molecular orbital coefficients and orbital energies 

422 from files mos, alpha and beta""" 

423 

424 ans = reader.read_mos(self.results) 

425 if ans is not None: 

426 self.converged = ans 

427 

428 def read_basis_set(self): 

429 """read the basis set""" 

430 reader.read_basis_set(self.results) 

431 

432 def read_ecps(self): 

433 """read the effective core potentials""" 

434 reader.read_ecps(self.results) 

435 

436 def read_gradient(self): 

437 """read all information in file 'gradient'""" 

438 reader.read_gradient(self.results) 

439 

440 def read_hessian(self): 

441 """Read in the hessian matrix""" 

442 reader.read_hessian(self.results) 

443 

444 def read_normal_modes(self): 

445 """Read in vibrational normal modes""" 

446 reader.read_normal_modes(self.results) 

447 

448 def read_vibrational_reduced_masses(self): 

449 """Read vibrational reduced masses""" 

450 reader.read_vibrational_reduced_masses(self.results) 

451 

452 def read_vibrational_spectrum(self): 

453 """Read the vibrational spectrum""" 

454 reader.read_vibrational_spectrum(self.results) 

455 

456 def read_ssquare(self): 

457 """Read the expectation value of S^2 operator""" 

458 reader.read_ssquare(self.results) 

459 

460 def read_dipole_moment(self): 

461 """Read the dipole moment""" 

462 reader.read_dipole_moment(self.results) 

463 dip_vec = self.results['electric dipole moment']['vector']['array'] 

464 self.dipole = np.array(dip_vec) * Bohr 

465 

466 def read_charges(self): 

467 """read partial charges on atoms from an ESP fit""" 

468 filename = get_output_filename(self.calculate_energy) 

469 self.charges = reader.read_charges(filename, len(self.atoms)) 

470 

471 def get_version(self): 

472 """get the version of the installed turbomole package""" 

473 if not self.version: 

474 self.version = reader.read_version(self.directory) 

475 return self.version 

476 

477 def get_datetime(self): 

478 """get the timestamp of most recent calculation""" 

479 if not self.datetime: 

480 self.datetime = reader.read_datetime(self.directory) 

481 return self.datetime 

482 

483 def get_runtime(self): 

484 """get the total runtime of calculations""" 

485 if not self.runtime: 

486 self.runtime = reader.read_runtime(self.directory) 

487 return self.runtime 

488 

489 def get_hostname(self): 

490 """get the hostname of the computer on which the calc has run""" 

491 if not self.hostname: 

492 self.hostname = reader.read_hostname(self.directory) 

493 return self.hostname 

494 

495 def get_optimizer(self, atoms, trajectory=None, logfile=None): 

496 """returns a TurbomoleOptimizer object""" 

497 self.parameters['task'] = 'optimize' 

498 self.parameters.verify() 

499 return TurbomoleOptimizer(atoms, self) 

500 

501 def get_results(self): 

502 """returns the results dictionary""" 

503 return self.results 

504 

505 def get_potential_energy(self, atoms, force_consistent=True): 

506 # update atoms 

507 self.updated = self.e_total is None 

508 self.set_atoms(atoms) 

509 self.initialize() 

510 # if update of energy is necessary 

511 if self.update_energy: 

512 # calculate energy 

513 execute([self.calculate_energy]) 

514 # check convergence 

515 self.converged = read_convergence(self.restart, self.parameters) 

516 if not self.converged: 

517 return None 

518 # read energy 

519 self.read_energy() 

520 

521 self.update_energy = False 

522 return self.e_total 

523 

524 def get_forces(self, atoms): 

525 # update atoms 

526 self.updated = self.forces is None 

527 self.set_atoms(atoms) 

528 # complete energy calculations 

529 if self.update_energy: 

530 self.get_potential_energy(atoms) 

531 # if update of forces is necessary 

532 if self.update_forces: 

533 # calculate forces 

534 execute([self.calculate_forces]) 

535 # read forces 

536 self.read_forces() 

537 

538 self.update_forces = False 

539 return self.forces.copy() 

540 

541 def get_dipole_moment(self, atoms): 

542 self.get_potential_energy(atoms) 

543 self.read_dipole_moment() 

544 return self.dipole 

545 

546 def get_property(self, name, atoms=None, allow_calculation=True): 

547 """return the value of a property""" 

548 

549 if name not in self.implemented_properties: 

550 # an ugly work around; the caller should test the raised error 

551 # if name in ['magmom', 'magmoms', 'charges', 'stress']: 

552 # return None 

553 raise PropertyNotImplementedError(name) 

554 

555 if atoms is None: 

556 atoms = self.atoms.copy() 

557 

558 persist_property = { 

559 'energy': 'e_total', 

560 'forces': 'forces', 

561 'dipole': 'dipole', 

562 'free_energy': 'e_total', 

563 'charges': 'charges' 

564 } 

565 property_getter = { 

566 'energy': self.get_potential_energy, 

567 'forces': self.get_forces, 

568 'dipole': self.get_dipole_moment, 

569 'free_energy': self.get_potential_energy, 

570 'charges': self.get_charges 

571 } 

572 getter_args = { 

573 'energy': [atoms], 

574 'forces': [atoms], 

575 'dipole': [atoms], 

576 'free_energy': [atoms, True], 

577 'charges': [atoms] 

578 } 

579 

580 if allow_calculation: 

581 result = property_getter[name](*getter_args[name]) 

582 else: 

583 if hasattr(self, persist_property[name]): 

584 result = getattr(self, persist_property[name]) 

585 else: 

586 result = None 

587 

588 if isinstance(result, np.ndarray): 

589 result = result.copy() 

590 return result 

591 

592 def get_charges(self, atoms): 

593 """return partial charges on atoms from an ESP fit""" 

594 self.get_potential_energy(atoms) 

595 self.read_charges() 

596 return self.charges 

597 

598 def get_forces_on_point_charges(self): 

599 """return forces acting on point charges""" 

600 self.get_forces(self.atoms) 

601 lines = read_data_group('point_charge_gradients').split('\n')[1:] 

602 forces = [] 

603 for line in lines: 

604 linef = line.strip().replace('D', 'E') 

605 forces.append([float(x) for x in linef.split()]) 

606 # Note the '-' sign for turbomole, to get forces 

607 return -np.array(forces) * Ha / Bohr 

608 

609 def set_point_charges(self, pcpot=None): 

610 """write external point charges to control""" 

611 if pcpot is not None and pcpot != self.pcpot: 

612 self.pcpot = pcpot 

613 if self.pcpot.mmcharges is None or self.pcpot.mmpositions is None: 

614 raise RuntimeError('external point charges not defined') 

615 

616 if not self.pc_initialized: 

617 if len(read_data_group('point_charges')) == 0: 

618 add_data_group('point_charges', 'file=pc.txt') 

619 if len(read_data_group('point_charge_gradients')) == 0: 

620 add_data_group( 

621 'point_charge_gradients', 

622 'file=pc_gradients.txt' 

623 ) 

624 drvopt = read_data_group('drvopt') 

625 if 'point charges' not in drvopt: 

626 drvopt += '\n point charges\n' 

627 delete_data_group('drvopt') 

628 add_data_group(drvopt, raw=True) 

629 self.pc_initialized = True 

630 

631 if self.pcpot.updated: 

632 with open('pc.txt', 'w') as pcfile: 

633 pcfile.write('$point_charges nocheck list\n') 

634 for (x, y, z), charge in zip( 

635 self.pcpot.mmpositions, self.pcpot.mmcharges): 

636 pcfile.write('%20.14f %20.14f %20.14f %20.14f\n' 

637 % (x / Bohr, y / Bohr, z / Bohr, charge)) 

638 pcfile.write('$end \n') 

639 self.pcpot.updated = False 

640 

641 def read_point_charges(self): 

642 """read point charges from previous calculation""" 

643 charges, positions = reader.read_point_charges() 

644 if len(charges) > 0: 

645 self.pcpot = PointChargePotential(charges, positions) 

646 

647 def embed(self, charges=None, positions=None): 

648 """embed atoms in an array of point-charges; function used in 

649 qmmm calculations.""" 

650 self.pcpot = PointChargePotential(charges, positions) 

651 return self.pcpot 

652 

653 

654class PointChargePotential: 

655 """Point-charge potential for Turbomole""" 

656 

657 def __init__(self, mmcharges, mmpositions=None): 

658 self.mmcharges = mmcharges 

659 self.mmpositions = mmpositions 

660 self.mmforces = None 

661 self.updated = True 

662 

663 def set_positions(self, mmpositions): 

664 """set the positions of point charges""" 

665 self.mmpositions = mmpositions 

666 self.updated = True 

667 

668 def set_charges(self, mmcharges): 

669 """set the values of point charges""" 

670 self.mmcharges = mmcharges 

671 self.updated = True 

672 

673 def get_forces(self, calc): 

674 """forces acting on point charges""" 

675 self.mmforces = calc.get_forces_on_point_charges() 

676 return self.mmforces