Coverage for /builds/ase/ase/ase/calculators/turbomole/turbomole.py: 27.15%

1# fmt: off

3# type: ignore

4"""

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

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

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

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

12"""

13import os

14import re

15import warnings

16from math import floor, log10

18import numpy as np

20from ase.calculators.calculator import (

21 Calculator,

22 PropertyNotImplementedError,

23 ReadError,

24)

25from ase.calculators.turbomole import reader

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

27from ase.calculators.turbomole.parameters import TurbomoleParameters

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

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

30from ase.io import read, write

31from ase.units import Bohr, Ha

34class TurbomoleOptimizer:

35 def __init__(self, atoms, calc):

36 self.atoms = atoms

37 self.calc = calc

38 self.atoms.calc = self.calc

40 def todict(self):

41 return {'type': 'optimization',

42 'optimizer': 'TurbomoleOptimizer'}

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

45 if fmax is not None:

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

47 self.calc.parameters.verify()

48 if steps is not None:

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

50 self.calc.parameters.verify()

51 self.calc.calculate()

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

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

56class Turbomole(Calculator):

58 """constants"""

59 name = 'Turbomole'

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

62 'charges']

64 tm_files = [

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

66 'alpha', 'beta', 'statistics', 'GEO_OPT_CONVERGED', 'GEO_OPT_FAILED',

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

68 'optinfo', 'statistics', 'converged', 'vibspectrum',

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

70 'pc_gradients.txt'

71 ]

72 tm_tmp_files = [

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

74 'diff_dft_density', 'diff_dft_oper', 'diff_fockmat', 'diis_errvec',

75 'diis_oldfock', 'dh'

76 ]

78 # initialize attributes

79 results = {}

80 initialized = False

81 pc_initialized = False

82 converged = False

83 updated = False

84 update_energy = None

85 update_forces = None

86 update_geometry = None

87 update_hessian = None

88 atoms = None

89 forces = None

90 e_total = None

91 dipole = None

92 charges = None

93 version = None

94 runtime = None

95 datetime = None

96 hostname = None

97 pcpot = None

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

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

101 restart=False, define_str=None, control_kdg=None,

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

103

104 super().__init__(label=label)

105 self.parameters = TurbomoleParameters()

106

107 self.calculate_energy = calculate_energy

108 self.calculate_forces = calculate_forces

109 self.post_HF = post_HF

110 self.restart = restart

111 self.parameters.define_str = define_str

112 self.control_kdg = control_kdg

113 self.control_input = control_input

114 self.reset_tolerance = reset_tolerance

115

116 if self.restart:

117 self._set_restart(kwargs)

118 else:

119 self.parameters.update(kwargs)

120 self.set_parameters()

121 self.parameters.verify()

122 self.reset()

123

124 if atoms is not None:

125 atoms.calc = self

126 self.set_atoms(atoms)

127

128 def __getitem__(self, item):

129 return getattr(self, item)

130

131 def _set_restart(self, params_update):

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

133 calculation"""

134

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

136 self.read_restart()

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

138

139 # update and verify parameters

140 self.parameters.update_restart(params_update)

141 self.set_parameters()

142 self.parameters.verify()

143

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

145 if self.parameters.define_str:

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

147

148 self._set_post_define()

149

150 self.initialized = True

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

152 self.update_energy = True

153 self.update_forces = True

154 self.update_geometry = True

155 self.update_hessian = True

156

157 def _set_post_define(self):

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

159 self.parameters.update_no_define_parameters()

160

161 # delete user-specified data groups

162 if self.control_kdg:

163 for dg in self.control_kdg:

164 delete_data_group(dg)

165

166 # append user-defined input to control

167 if self.control_input:

168 for inp in self.control_input:

169 add_data_group(inp, raw=True)

170

171 # add point charges if pcpot defined:

172 if self.pcpot:

173 self.set_point_charges()

174

175 def set_parameters(self):

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

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

178 self.calculate_energy = 'ridft'

179 self.calculate_forces = 'rdgrad'

180

181 def reset(self):

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

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

184 self.atoms = None

185 self.results = {}

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

187 ase_files = [

188 f for f in os.listdir(

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

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

191 if os.path.exists(f):

192 os.remove(f)

193 self.initialized = False

194 self.pc_initialized = False

195 self.converged = False

196

197 def set_atoms(self, atoms):

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

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

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

201 version.

202 """

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

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

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

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

207 self.updated = False

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

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

210 return

211 else:

212 return

213

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

215 if 'positions' in changes:

216 # print(two atoms obj are different')

217 self.reset()

218 else:

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

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

221 self.reset()

222

223 write('coord', atoms)

224 self.atoms = atoms.copy()

225 self.update_energy = True

226 self.update_forces = True

227 self.update_geometry = True

228 self.update_hessian = True

229

230 def initialize(self):

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

232 if self.initialized:

233 return

234 self.parameters.verify()

235 if not self.atoms:

236 raise RuntimeError('atoms missing during initialization')

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

238 raise OSError('file coord not found')

239

240 # run define

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

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

243

244 # process non-default initial guess

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

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

247 # "use" initial guess

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

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

250 else:

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

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

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

254 # "hcore" initial guess

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

256 delete_data_group('uhfmo_alpha')

257 delete_data_group('uhfmo_beta')

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

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

260 else:

261 delete_data_group('scfmo')

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

263

264 self._set_post_define()

265

266 self.initialized = True

267 self.converged = False

268

269 def calculation_required(self, atoms, properties):

270 if self.atoms != atoms:

271 return True

272 for prop in properties:

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

274 return True

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

276 return True

277 return False

278

279 def calculate(self, atoms=None):

280 """execute the requested job"""

281 if atoms is None:

282 atoms = self.atoms

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

284 self.get_potential_energy(atoms)

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

286 self.get_forces(atoms)

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

288 self.relax_geometry(atoms)

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

290 self.normal_mode_analysis(atoms)

291 self.read_results()

292

293 def relax_geometry(self, atoms=None):

294 """execute geometry optimization with script jobex"""

295 if atoms is None:

296 atoms = self.atoms

297 self.set_atoms(atoms)

298 if self.converged and not self.update_geometry:

299 return

300 self.initialize()

301 jobex_command = ['jobex']

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

303 jobex_command.append('-trans')

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

305 jobex_command.append('-ri')

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

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

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

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

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

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

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

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

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

315 if geom_iter is not None:

316 assert isinstance(geom_iter, int)

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

318 self.converged = False

319 execute(jobex_command)

320 # check convergence

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

322 if self.converged:

323 self.update_energy = False

324 self.update_forces = False

325 self.update_geometry = False

326 self.update_hessian = True

327 # read results

328 new_struct = read('coord')

329 atoms.set_positions(new_struct.get_positions())

330 self.atoms = atoms.copy()

331 self.read_energy()

332

333 def normal_mode_analysis(self, atoms=None):

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

335 from ase.constraints import FixAtoms

336 if atoms is None:

337 atoms = self.atoms

338 self.set_atoms(atoms)

339 self.initialize()

340 if self.update_energy:

341 self.get_potential_energy(atoms)

342 if self.update_hessian:

343 fixatoms = []

344 for constr in atoms.constraints:

345 if isinstance(constr, FixAtoms):

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

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

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

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

350 if len(fixatoms) > 0:

351 define_str = '\n\ny\n'

352 for index in fixatoms:

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

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

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

356 dg = read_data_group('atoms')

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

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

359 dg += '\n'

360 delete_data_group('atoms')

361 add_data_group(dg, raw=True)

362 execute(['aoforce'])

363 else:

364 nforce_cmd = ['NumForce']

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

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

367 nforce_cmd.append('-ri')

368 if len(fixatoms) > 0:

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

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

371 nforce_cmd.append('-central')

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

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

374 if self.update_forces:

375 self.get_forces(atoms)

376 execute(nforce_cmd)

377 self.update_hessian = False

378

379 def read_restart(self):

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

381 self.atoms = read('coord')

382 self.atoms.calc = self

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

384 self.read_results()

385

386 def read_results(self):

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

388 read_methods = [

389 self.read_energy,

390 self.read_gradient,

391 self.read_forces,

392 self.read_basis_set,

393 self.read_ecps,

394 self.read_mos,

395 self.read_occupation_numbers,

396 self.read_dipole_moment,

397 self.read_ssquare,

398 self.read_hessian,

399 self.read_vibrational_reduced_masses,

400 self.read_normal_modes,

401 self.read_vibrational_spectrum,

402 self.read_charges,

403 self.read_point_charges,

404 self.read_run_parameters

405 ]

406 for method in read_methods:

407 try:

408 method()

409 except ReadError as err:

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

411

412 def read_run_parameters(self):

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

414 reader.read_run_parameters(self.results)

415

416 def read_energy(self):

417 """Read energy from Turbomole energy file."""

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

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

420

421 def read_forces(self):

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

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

424

425 def read_occupation_numbers(self):

426 """read occupation numbers"""

427 reader.read_occupation_numbers(self.results)

428

429 def read_mos(self):

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

431 from files mos, alpha and beta"""

432

433 ans = reader.read_mos(self.results)

434 if ans is not None:

435 self.converged = ans

436

437 def read_basis_set(self):

438 """read the basis set"""

439 reader.read_basis_set(self.results)

440

441 def read_ecps(self):

442 """read the effective core potentials"""

443 reader.read_ecps(self.results)

444

445 def read_gradient(self):

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

447 reader.read_gradient(self.results)

448

449 def read_hessian(self):

450 """Read in the hessian matrix"""

451 reader.read_hessian(self.results)

452

453 def read_normal_modes(self):

454 """Read in vibrational normal modes"""

455 reader.read_normal_modes(self.results)

456

457 def read_vibrational_reduced_masses(self):

458 """Read vibrational reduced masses"""

459 reader.read_vibrational_reduced_masses(self.results)

460

461 def read_vibrational_spectrum(self):

462 """Read the vibrational spectrum"""

463 reader.read_vibrational_spectrum(self.results)

464

465 def read_ssquare(self):

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

467 reader.read_ssquare(self.results)

468

469 def read_dipole_moment(self):

470 """Read the dipole moment"""

471 reader.read_dipole_moment(self.results)

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

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

474

475 def read_charges(self):

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

477 filename = get_output_filename(self.calculate_energy)

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

479

480 def get_version(self):

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

482 if not self.version:

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

484 return self.version

485

486 def get_datetime(self):

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

488 if not self.datetime:

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

490 return self.datetime

491

492 def get_runtime(self):

493 """get the total runtime of calculations"""

494 if not self.runtime:

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

496 return self.runtime

497

498 def get_hostname(self):

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

500 if not self.hostname:

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

502 return self.hostname

503

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

505 """returns a TurbomoleOptimizer object"""

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

507 self.parameters.verify()

508 return TurbomoleOptimizer(atoms, self)

509

510 def get_results(self):

511 """returns the results dictionary"""

512 return self.results

513

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

515 # update atoms

516 self.updated = self.e_total is None

517 self.set_atoms(atoms)

518 self.initialize()

519 # if update of energy is necessary

520 if self.update_energy:

521 # calculate energy

522 execute([self.calculate_energy])

523 # check convergence

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

525 if not self.converged:

526 return None

527 # read energy

528 self.read_energy()

529

530 self.update_energy = False

531 return self.e_total

532

533 def get_forces(self, atoms):

534 # update atoms

535 self.updated = self.forces is None

536 self.set_atoms(atoms)

537 # complete energy calculations

538 if self.update_energy:

539 self.get_potential_energy(atoms)

540 # if update of forces is necessary

541 if self.update_forces:

542 # calculate forces

543 execute([self.calculate_forces])

544 # read forces

545 self.read_forces()

546

547 self.update_forces = False

548 return self.forces.copy()

549

550 def get_dipole_moment(self, atoms):

551 self.get_potential_energy(atoms)

552 self.read_dipole_moment()

553 return self.dipole

554

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

556 """return the value of a property"""

557

558 if name not in self.implemented_properties:

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

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

561 # return None

562 raise PropertyNotImplementedError(name)

563

564 if atoms is None:

565 atoms = self.atoms.copy()

566

567 persist_property = {

568 'energy': 'e_total',

569 'forces': 'forces',

570 'dipole': 'dipole',

571 'free_energy': 'e_total',

572 'charges': 'charges'

573 }

574 property_getter = {

575 'energy': self.get_potential_energy,

576 'forces': self.get_forces,

577 'dipole': self.get_dipole_moment,

578 'free_energy': self.get_potential_energy,

579 'charges': self.get_charges

580 }

581 getter_args = {

582 'energy': [atoms],

583 'forces': [atoms],

584 'dipole': [atoms],

585 'free_energy': [atoms, True],

586 'charges': [atoms]

587 }

588

589 if allow_calculation:

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

591 else:

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

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

594 else:

595 result = None

596

597 if isinstance(result, np.ndarray):

598 result = result.copy()

599 return result

600

601 def get_charges(self, atoms):

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

603 self.get_potential_energy(atoms)

604 self.read_charges()

605 return self.charges

606

607 def get_forces_on_point_charges(self):

608 """return forces acting on point charges"""

609 self.get_forces(self.atoms)

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

611 forces = []

612 for line in lines:

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

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

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

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

617

618 def set_point_charges(self, pcpot=None):

619 """write external point charges to control"""

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

621 self.pcpot = pcpot

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

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

624

625 if not self.pc_initialized:

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

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

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

629 add_data_group(

630 'point_charge_gradients',

631 'file=pc_gradients.txt'

632 )

633 drvopt = read_data_group('drvopt')

634 if 'point charges' not in drvopt:

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

636 delete_data_group('drvopt')

637 add_data_group(drvopt, raw=True)

638 self.pc_initialized = True

639

640 if self.pcpot.updated:

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

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

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

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

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

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

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

648 self.pcpot.updated = False

649

650 def read_point_charges(self):

651 """read point charges from previous calculation"""

652 charges, positions = reader.read_point_charges()

653 if len(charges) > 0:

654 self.pcpot = PointChargePotential(charges, positions)

655

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

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

658 qmmm calculations."""

659 self.pcpot = PointChargePotential(charges, positions)

660 return self.pcpot

661

662

663class PointChargePotential:

664 """Point-charge potential for Turbomole"""

665

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

667 self.mmcharges = mmcharges

668 self.mmpositions = mmpositions

669 self.mmforces = None

670 self.updated = True

671

672 def set_positions(self, mmpositions):

673 """set the positions of point charges"""

674 self.mmpositions = mmpositions

675 self.updated = True

676

677 def set_charges(self, mmcharges):

678 """set the values of point charges"""

679 self.mmcharges = mmcharges

680 self.updated = True

681

682 def get_forces(self, calc):

683 """forces acting on point charges"""

684 self.mmforces = calc.get_forces_on_point_charges()

685 return self.mmforces