Coverage for /builds/kinetik161/ase/ase/calculators/turbomole/turbomole.py: 27.31%

1# type: ignore

2"""

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

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

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

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

10"""

11import os

12import re

13import warnings

14from math import floor, log10

16import numpy as np

18from ase.calculators.calculator import (Calculator,

19 PropertyNotImplementedError, ReadError)

20from ase.calculators.turbomole import reader

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

22from ase.calculators.turbomole.parameters import TurbomoleParameters

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

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

25from ase.io import read, write

26from ase.units import Bohr, Ha

29class TurbomoleOptimizer:

30 def __init__(self, atoms, calc):

31 self.atoms = atoms

32 self.calc = calc

33 self.atoms.calc = self.calc

35 def todict(self):

36 return {'type': 'optimization',

37 'optimizer': 'TurbomoleOptimizer'}

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

40 if fmax is not None:

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

42 self.calc.parameters.verify()

43 if steps is not None:

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

45 self.calc.parameters.verify()

46 self.calc.calculate()

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

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

51class Turbomole(Calculator):

53 """constants"""

54 name = 'Turbomole'

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

57 'charges']

59 tm_files = [

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

61 'alpha', 'beta', 'statistics', 'GEO_OPT_CONVERGED', 'GEO_OPT_FAILED',

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

63 'optinfo', 'statistics', 'converged', 'vibspectrum',

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

65 'pc_gradients.txt'

66 ]

67 tm_tmp_files = [

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

69 'diff_dft_density', 'diff_dft_oper', 'diff_fockmat', 'diis_errvec',

70 'diis_oldfock', 'dh'

71 ]

73 # initialize attributes

74 results = {}

75 initialized = False

76 pc_initialized = False

77 converged = False

78 updated = False

79 update_energy = None

80 update_forces = None

81 update_geometry = None

82 update_hessian = None

83 atoms = None

84 forces = None

85 e_total = None

86 dipole = None

87 charges = None

88 version = None

89 runtime = None

90 datetime = None

91 hostname = None

92 pcpot = None

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

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

96 restart=False, define_str=None, control_kdg=None,

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

99 super().__init__(label=label)

100 self.parameters = TurbomoleParameters()

101

102 self.calculate_energy = calculate_energy

103 self.calculate_forces = calculate_forces

104 self.post_HF = post_HF

105 self.restart = restart

106 self.parameters.define_str = define_str

107 self.control_kdg = control_kdg

108 self.control_input = control_input

109 self.reset_tolerance = reset_tolerance

110

111 if self.restart:

112 self._set_restart(kwargs)

113 else:

114 self.parameters.update(kwargs)

115 self.set_parameters()

116 self.parameters.verify()

117 self.reset()

118

119 if atoms is not None:

120 atoms.calc = self

121 self.set_atoms(atoms)

122

123 def __getitem__(self, item):

124 return getattr(self, item)

125

126 def _set_restart(self, params_update):

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

128 calculation"""

129

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

131 self.read_restart()

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

133

134 # update and verify parameters

135 self.parameters.update_restart(params_update)

136 self.set_parameters()

137 self.parameters.verify()

138

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

140 if self.parameters.define_str:

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

142

143 self._set_post_define()

144

145 self.initialized = True

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

147 self.update_energy = True

148 self.update_forces = True

149 self.update_geometry = True

150 self.update_hessian = True

151

152 def _set_post_define(self):

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

154 self.parameters.update_no_define_parameters()

155

156 # delete user-specified data groups

157 if self.control_kdg:

158 for dg in self.control_kdg:

159 delete_data_group(dg)

160

161 # append user-defined input to control

162 if self.control_input:

163 for inp in self.control_input:

164 add_data_group(inp, raw=True)

165

166 # add point charges if pcpot defined:

167 if self.pcpot:

168 self.set_point_charges()

169

170 def set_parameters(self):

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

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

173 self.calculate_energy = 'ridft'

174 self.calculate_forces = 'rdgrad'

175

176 def reset(self):

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

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

179 self.atoms = None

180 self.results = {}

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

182 ase_files = [

183 f for f in os.listdir(

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

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

186 if os.path.exists(f):

187 os.remove(f)

188 self.initialized = False

189 self.pc_initialized = False

190 self.converged = False

191

192 def set_atoms(self, atoms):

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

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

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

196 version.

197 """

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

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

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

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

202 self.updated = False

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

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

205 return

206 else:

207 return

208

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

210 if 'positions' in changes:

211 # print(two atoms obj are different')

212 self.reset()

213 else:

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

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

216 self.reset()

217

218 write('coord', atoms)

219 self.atoms = atoms.copy()

220 self.update_energy = True

221 self.update_forces = True

222 self.update_geometry = True

223 self.update_hessian = True

224

225 def initialize(self):

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

227 if self.initialized:

228 return

229 self.parameters.verify()

230 if not self.atoms:

231 raise RuntimeError('atoms missing during initialization')

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

233 raise OSError('file coord not found')

234

235 # run define

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

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

238

239 # process non-default initial guess

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

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

242 # "use" initial guess

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

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

245 else:

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

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

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

249 # "hcore" initial guess

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

251 delete_data_group('uhfmo_alpha')

252 delete_data_group('uhfmo_beta')

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

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

255 else:

256 delete_data_group('scfmo')

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

258

259 self._set_post_define()

260

261 self.initialized = True

262 self.converged = False

263

264 def calculation_required(self, atoms, properties):

265 if self.atoms != atoms:

266 return True

267 for prop in properties:

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

269 return True

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

271 return True

272 return False

273

274 def calculate(self, atoms=None):

275 """execute the requested job"""

276 if atoms is None:

277 atoms = self.atoms

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

279 self.get_potential_energy(atoms)

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

281 self.get_forces(atoms)

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

283 self.relax_geometry(atoms)

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

285 self.normal_mode_analysis(atoms)

286 self.read_results()

287

288 def relax_geometry(self, atoms=None):

289 """execute geometry optimization with script jobex"""

290 if atoms is None:

291 atoms = self.atoms

292 self.set_atoms(atoms)

293 if self.converged and not self.update_geometry:

294 return

295 self.initialize()

296 jobex_command = ['jobex']

297 if self.parameters['transition vector']:

298 jobex_command.append('-trans')

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

300 jobex_command.append('-ri')

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

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

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

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

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

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

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

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

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

310 if geom_iter is not None:

311 assert isinstance(geom_iter, int)

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

313 self.converged = False

314 execute(jobex_command)

315 # check convergence

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

317 if self.converged:

318 self.update_energy = False

319 self.update_forces = False

320 self.update_geometry = False

321 self.update_hessian = True

322 # read results

323 new_struct = read('coord')

324 atoms.set_positions(new_struct.get_positions())

325 self.atoms = atoms.copy()

326 self.read_energy()

327

328 def normal_mode_analysis(self, atoms=None):

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

330 from ase.constraints import FixAtoms

331 if atoms is None:

332 atoms = self.atoms

333 self.set_atoms(atoms)

334 self.initialize()

335 if self.update_energy:

336 self.get_potential_energy(atoms)

337 if self.update_hessian:

338 fixatoms = []

339 for constr in atoms.constraints:

340 if isinstance(constr, FixAtoms):

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

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

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

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

345 if len(fixatoms) > 0:

346 define_str = '\n\ny\n'

347 for index in fixatoms:

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

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

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

351 dg = read_data_group('atoms')

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

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

354 dg += '\n'

355 delete_data_group('atoms')

356 add_data_group(dg, raw=True)

357 execute(['aoforce'])

358 else:

359 nforce_cmd = ['NumForce']

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

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

362 nforce_cmd.append('-ri')

363 if len(fixatoms) > 0:

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

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

366 nforce_cmd.append('-central')

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

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

369 if self.update_forces:

370 self.get_forces(atoms)

371 execute(nforce_cmd)

372 self.update_hessian = False

373

374 def read_restart(self):

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

376 self.atoms = read('coord')

377 self.atoms.calc = self

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

379 self.read_results()

380

381 def read_results(self):

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

383 read_methods = [

384 self.read_energy,

385 self.read_gradient,

386 self.read_forces,

387 self.read_basis_set,

388 self.read_ecps,

389 self.read_mos,

390 self.read_occupation_numbers,

391 self.read_dipole_moment,

392 self.read_ssquare,

393 self.read_hessian,

394 self.read_vibrational_reduced_masses,

395 self.read_normal_modes,

396 self.read_vibrational_spectrum,

397 self.read_charges,

398 self.read_point_charges,

399 self.read_run_parameters

400 ]

401 for method in read_methods:

402 try:

403 method()

404 except ReadError as err:

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

406

407 def read_run_parameters(self):

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

409 reader.read_run_parameters(self.results)

410

411 def read_energy(self):

412 """Read energy from Turbomole energy file."""

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

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

415

416 def read_forces(self):

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

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

419

420 def read_occupation_numbers(self):

421 """read occupation numbers"""

422 reader.read_occupation_numbers(self.results)

423

424 def read_mos(self):

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

426 from files mos, alpha and beta"""

427

428 ans = reader.read_mos(self.results)

429 if ans is not None:

430 self.converged = ans

431

432 def read_basis_set(self):

433 """read the basis set"""

434 reader.read_basis_set(self.results)

435

436 def read_ecps(self):

437 """read the effective core potentials"""

438 reader.read_ecps(self.results)

439

440 def read_gradient(self):

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

442 reader.read_gradient(self.results)

443

444 def read_hessian(self):

445 """Read in the hessian matrix"""

446 reader.read_hessian(self.results)

447

448 def read_normal_modes(self):

449 """Read in vibrational normal modes"""

450 reader.read_normal_modes(self.results)

451

452 def read_vibrational_reduced_masses(self):

453 """Read vibrational reduced masses"""

454 reader.read_vibrational_reduced_masses(self.results)

455

456 def read_vibrational_spectrum(self):

457 """Read the vibrational spectrum"""

458 reader.read_vibrational_spectrum(self.results)

459

460 def read_ssquare(self):

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

462 reader.read_ssquare(self.results)

463

464 def read_dipole_moment(self):

465 """Read the dipole moment"""

466 reader.read_dipole_moment(self.results)

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

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

469

470 def read_charges(self):

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

472 filename = get_output_filename(self.calculate_energy)

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

474

475 def get_version(self):

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

477 if not self.version:

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

479 return self.version

480

481 def get_datetime(self):

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

483 if not self.datetime:

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

485 return self.datetime

486

487 def get_runtime(self):

488 """get the total runtime of calculations"""

489 if not self.runtime:

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

491 return self.runtime

492

493 def get_hostname(self):

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

495 if not self.hostname:

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

497 return self.hostname

498

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

500 """returns a TurbomoleOptimizer object"""

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

502 self.parameters.verify()

503 return TurbomoleOptimizer(atoms, self)

504

505 def get_results(self):

506 """returns the results dictionary"""

507 return self.results

508

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

510 # update atoms

511 self.updated = self.e_total is None

512 self.set_atoms(atoms)

513 self.initialize()

514 # if update of energy is necessary

515 if self.update_energy:

516 # calculate energy

517 execute([self.calculate_energy])

518 # check convergence

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

520 if not self.converged:

521 return None

522 # read energy

523 self.read_energy()

524

525 self.update_energy = False

526 return self.e_total

527

528 def get_forces(self, atoms):

529 # update atoms

530 self.updated = self.forces is None

531 self.set_atoms(atoms)

532 # complete energy calculations

533 if self.update_energy:

534 self.get_potential_energy(atoms)

535 # if update of forces is necessary

536 if self.update_forces:

537 # calculate forces

538 execute([self.calculate_forces])

539 # read forces

540 self.read_forces()

541

542 self.update_forces = False

543 return self.forces.copy()

544

545 def get_dipole_moment(self, atoms):

546 self.get_potential_energy(atoms)

547 self.read_dipole_moment()

548 return self.dipole

549

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

551 """return the value of a property"""

552

553 if name not in self.implemented_properties:

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

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

556 # return None

557 raise PropertyNotImplementedError(name)

558

559 if atoms is None:

560 atoms = self.atoms.copy()

561

562 persist_property = {

563 'energy': 'e_total',

564 'forces': 'forces',

565 'dipole': 'dipole',

566 'free_energy': 'e_total',

567 'charges': 'charges'

568 }

569 property_getter = {

570 'energy': self.get_potential_energy,

571 'forces': self.get_forces,

572 'dipole': self.get_dipole_moment,

573 'free_energy': self.get_potential_energy,

574 'charges': self.get_charges

575 }

576 getter_args = {

577 'energy': [atoms],

578 'forces': [atoms],

579 'dipole': [atoms],

580 'free_energy': [atoms, True],

581 'charges': [atoms]

582 }

583

584 if allow_calculation:

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

586 else:

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

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

589 else:

590 result = None

591

592 if isinstance(result, np.ndarray):

593 result = result.copy()

594 return result

595

596 def get_charges(self, atoms):

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

598 self.get_potential_energy(atoms)

599 self.read_charges()

600 return self.charges

601

602 def get_forces_on_point_charges(self):

603 """return forces acting on point charges"""

604 self.get_forces(self.atoms)

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

606 forces = []

607 for line in lines:

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

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

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

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

612

613 def set_point_charges(self, pcpot=None):

614 """write external point charges to control"""

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

616 self.pcpot = pcpot

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

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

619

620 if not self.pc_initialized:

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

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

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

624 add_data_group(

625 'point_charge_gradients',

626 'file=pc_gradients.txt'

627 )

628 drvopt = read_data_group('drvopt')

629 if 'point charges' not in drvopt:

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

631 delete_data_group('drvopt')

632 add_data_group(drvopt, raw=True)

633 self.pc_initialized = True

634

635 if self.pcpot.updated:

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

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

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

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

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

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

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

643 self.pcpot.updated = False

644

645 def read_point_charges(self):

646 """read point charges from previous calculation"""

647 charges, positions = reader.read_point_charges()

648 if len(charges) > 0:

649 self.pcpot = PointChargePotential(charges, positions)

650

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

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

653 qmmm calculations."""

654 self.pcpot = PointChargePotential(charges, positions)

655 return self.pcpot

656

657

658class PointChargePotential:

659 """Point-charge potential for Turbomole"""

660

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

662 self.mmcharges = mmcharges

663 self.mmpositions = mmpositions

664 self.mmforces = None

665 self.updated = True

666

667 def set_positions(self, mmpositions):

668 """set the positions of point charges"""

669 self.mmpositions = mmpositions

670 self.updated = True

671

672 def set_charges(self, mmcharges):

673 """set the values of point charges"""

674 self.mmcharges = mmcharges

675 self.updated = True

676

677 def get_forces(self, calc):

678 """forces acting on point charges"""

679 self.mmforces = calc.get_forces_on_point_charges()

680 return self.mmforces