Coverage for /builds/debichem-team/python-ase/ase/calculators/turbomole/turbomole.py: 27.15%

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 (

19 Calculator,

20 PropertyNotImplementedError,

21 ReadError,

22)

23from ase.calculators.turbomole import reader

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

25from ase.calculators.turbomole.parameters import TurbomoleParameters

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

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

28from ase.io import read, write

29from ase.units import Bohr, Ha

32class TurbomoleOptimizer:

33 def __init__(self, atoms, calc):

34 self.atoms = atoms

35 self.calc = calc

36 self.atoms.calc = self.calc

38 def todict(self):

39 return {'type': 'optimization',

40 'optimizer': 'TurbomoleOptimizer'}

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

43 if fmax is not None:

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

45 self.calc.parameters.verify()

46 if steps is not None:

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

48 self.calc.parameters.verify()

49 self.calc.calculate()

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

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

54class Turbomole(Calculator):

56 """constants"""

57 name = 'Turbomole'

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

60 'charges']

62 tm_files = [

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

64 'alpha', 'beta', 'statistics', 'GEO_OPT_CONVERGED', 'GEO_OPT_FAILED',

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

66 'optinfo', 'statistics', 'converged', 'vibspectrum',

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

68 'pc_gradients.txt'

69 ]

70 tm_tmp_files = [

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

72 'diff_dft_density', 'diff_dft_oper', 'diff_fockmat', 'diis_errvec',

73 'diis_oldfock', 'dh'

74 ]

76 # initialize attributes

77 results = {}

78 initialized = False

79 pc_initialized = False

80 converged = False

81 updated = False

82 update_energy = None

83 update_forces = None

84 update_geometry = None

85 update_hessian = None

86 atoms = None

87 forces = None

88 e_total = None

89 dipole = None

90 charges = None

91 version = None

92 runtime = None

93 datetime = None

94 hostname = None

95 pcpot = None

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

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

99 restart=False, define_str=None, control_kdg=None,

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

101

102 super().__init__(label=label)

103 self.parameters = TurbomoleParameters()

104

105 self.calculate_energy = calculate_energy

106 self.calculate_forces = calculate_forces

107 self.post_HF = post_HF

108 self.restart = restart

109 self.parameters.define_str = define_str

110 self.control_kdg = control_kdg

111 self.control_input = control_input

112 self.reset_tolerance = reset_tolerance

113

114 if self.restart:

115 self._set_restart(kwargs)

116 else:

117 self.parameters.update(kwargs)

118 self.set_parameters()

119 self.parameters.verify()

120 self.reset()

121

122 if atoms is not None:

123 atoms.calc = self

124 self.set_atoms(atoms)

125

126 def __getitem__(self, item):

127 return getattr(self, item)

128

129 def _set_restart(self, params_update):

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

131 calculation"""

132

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

134 self.read_restart()

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

136

137 # update and verify parameters

138 self.parameters.update_restart(params_update)

139 self.set_parameters()

140 self.parameters.verify()

141

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

143 if self.parameters.define_str:

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

145

146 self._set_post_define()

147

148 self.initialized = True

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

150 self.update_energy = True

151 self.update_forces = True

152 self.update_geometry = True

153 self.update_hessian = True

154

155 def _set_post_define(self):

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

157 self.parameters.update_no_define_parameters()

158

159 # delete user-specified data groups

160 if self.control_kdg:

161 for dg in self.control_kdg:

162 delete_data_group(dg)

163

164 # append user-defined input to control

165 if self.control_input:

166 for inp in self.control_input:

167 add_data_group(inp, raw=True)

168

169 # add point charges if pcpot defined:

170 if self.pcpot:

171 self.set_point_charges()

172

173 def set_parameters(self):

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

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

176 self.calculate_energy = 'ridft'

177 self.calculate_forces = 'rdgrad'

178

179 def reset(self):

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

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

182 self.atoms = None

183 self.results = {}

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

185 ase_files = [

186 f for f in os.listdir(

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

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

189 if os.path.exists(f):

190 os.remove(f)

191 self.initialized = False

192 self.pc_initialized = False

193 self.converged = False

194

195 def set_atoms(self, atoms):

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

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

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

199 version.

200 """

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

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

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

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

205 self.updated = False

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

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

208 return

209 else:

210 return

211

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

213 if 'positions' in changes:

214 # print(two atoms obj are different')

215 self.reset()

216 else:

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

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

219 self.reset()

220

221 write('coord', atoms)

222 self.atoms = atoms.copy()

223 self.update_energy = True

224 self.update_forces = True

225 self.update_geometry = True

226 self.update_hessian = True

227

228 def initialize(self):

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

230 if self.initialized:

231 return

232 self.parameters.verify()

233 if not self.atoms:

234 raise RuntimeError('atoms missing during initialization')

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

236 raise OSError('file coord not found')

237

238 # run define

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

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

241

242 # process non-default initial guess

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

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

245 # "use" initial guess

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

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

248 else:

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

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

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

252 # "hcore" initial guess

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

254 delete_data_group('uhfmo_alpha')

255 delete_data_group('uhfmo_beta')

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

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

258 else:

259 delete_data_group('scfmo')

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

261

262 self._set_post_define()

263

264 self.initialized = True

265 self.converged = False

266

267 def calculation_required(self, atoms, properties):

268 if self.atoms != atoms:

269 return True

270 for prop in properties:

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

272 return True

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

274 return True

275 return False

276

277 def calculate(self, atoms=None):

278 """execute the requested job"""

279 if atoms is None:

280 atoms = self.atoms

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

282 self.get_potential_energy(atoms)

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

284 self.get_forces(atoms)

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

286 self.relax_geometry(atoms)

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

288 self.normal_mode_analysis(atoms)

289 self.read_results()

290

291 def relax_geometry(self, atoms=None):

292 """execute geometry optimization with script jobex"""

293 if atoms is None:

294 atoms = self.atoms

295 self.set_atoms(atoms)

296 if self.converged and not self.update_geometry:

297 return

298 self.initialize()

299 jobex_command = ['jobex']

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

301 jobex_command.append('-trans')

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

303 jobex_command.append('-ri')

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

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

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

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

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

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

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

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

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

313 if geom_iter is not None:

314 assert isinstance(geom_iter, int)

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

316 self.converged = False

317 execute(jobex_command)

318 # check convergence

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

320 if self.converged:

321 self.update_energy = False

322 self.update_forces = False

323 self.update_geometry = False

324 self.update_hessian = True

325 # read results

326 new_struct = read('coord')

327 atoms.set_positions(new_struct.get_positions())

328 self.atoms = atoms.copy()

329 self.read_energy()

330

331 def normal_mode_analysis(self, atoms=None):

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

333 from ase.constraints import FixAtoms

334 if atoms is None:

335 atoms = self.atoms

336 self.set_atoms(atoms)

337 self.initialize()

338 if self.update_energy:

339 self.get_potential_energy(atoms)

340 if self.update_hessian:

341 fixatoms = []

342 for constr in atoms.constraints:

343 if isinstance(constr, FixAtoms):

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

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

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

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

348 if len(fixatoms) > 0:

349 define_str = '\n\ny\n'

350 for index in fixatoms:

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

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

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

354 dg = read_data_group('atoms')

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

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

357 dg += '\n'

358 delete_data_group('atoms')

359 add_data_group(dg, raw=True)

360 execute(['aoforce'])

361 else:

362 nforce_cmd = ['NumForce']

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

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

365 nforce_cmd.append('-ri')

366 if len(fixatoms) > 0:

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

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

369 nforce_cmd.append('-central')

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

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

372 if self.update_forces:

373 self.get_forces(atoms)

374 execute(nforce_cmd)

375 self.update_hessian = False

376

377 def read_restart(self):

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

379 self.atoms = read('coord')

380 self.atoms.calc = self

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

382 self.read_results()

383

384 def read_results(self):

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

386 read_methods = [

387 self.read_energy,

388 self.read_gradient,

389 self.read_forces,

390 self.read_basis_set,

391 self.read_ecps,

392 self.read_mos,

393 self.read_occupation_numbers,

394 self.read_dipole_moment,

395 self.read_ssquare,

396 self.read_hessian,

397 self.read_vibrational_reduced_masses,

398 self.read_normal_modes,

399 self.read_vibrational_spectrum,

400 self.read_charges,

401 self.read_point_charges,

402 self.read_run_parameters

403 ]

404 for method in read_methods:

405 try:

406 method()

407 except ReadError as err:

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

409

410 def read_run_parameters(self):

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

412 reader.read_run_parameters(self.results)

413

414 def read_energy(self):

415 """Read energy from Turbomole energy file."""

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

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

418

419 def read_forces(self):

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

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

422

423 def read_occupation_numbers(self):

424 """read occupation numbers"""

425 reader.read_occupation_numbers(self.results)

426

427 def read_mos(self):

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

429 from files mos, alpha and beta"""

430

431 ans = reader.read_mos(self.results)

432 if ans is not None:

433 self.converged = ans

434

435 def read_basis_set(self):

436 """read the basis set"""

437 reader.read_basis_set(self.results)

438

439 def read_ecps(self):

440 """read the effective core potentials"""

441 reader.read_ecps(self.results)

442

443 def read_gradient(self):

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

445 reader.read_gradient(self.results)

446

447 def read_hessian(self):

448 """Read in the hessian matrix"""

449 reader.read_hessian(self.results)

450

451 def read_normal_modes(self):

452 """Read in vibrational normal modes"""

453 reader.read_normal_modes(self.results)

454

455 def read_vibrational_reduced_masses(self):

456 """Read vibrational reduced masses"""

457 reader.read_vibrational_reduced_masses(self.results)

458

459 def read_vibrational_spectrum(self):

460 """Read the vibrational spectrum"""

461 reader.read_vibrational_spectrum(self.results)

462

463 def read_ssquare(self):

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

465 reader.read_ssquare(self.results)

466

467 def read_dipole_moment(self):

468 """Read the dipole moment"""

469 reader.read_dipole_moment(self.results)

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

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

472

473 def read_charges(self):

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

475 filename = get_output_filename(self.calculate_energy)

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

477

478 def get_version(self):

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

480 if not self.version:

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

482 return self.version

483

484 def get_datetime(self):

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

486 if not self.datetime:

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

488 return self.datetime

489

490 def get_runtime(self):

491 """get the total runtime of calculations"""

492 if not self.runtime:

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

494 return self.runtime

495

496 def get_hostname(self):

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

498 if not self.hostname:

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

500 return self.hostname

501

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

503 """returns a TurbomoleOptimizer object"""

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

505 self.parameters.verify()

506 return TurbomoleOptimizer(atoms, self)

507

508 def get_results(self):

509 """returns the results dictionary"""

510 return self.results

511

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

513 # update atoms

514 self.updated = self.e_total is None

515 self.set_atoms(atoms)

516 self.initialize()

517 # if update of energy is necessary

518 if self.update_energy:

519 # calculate energy

520 execute([self.calculate_energy])

521 # check convergence

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

523 if not self.converged:

524 return None

525 # read energy

526 self.read_energy()

527

528 self.update_energy = False

529 return self.e_total

530

531 def get_forces(self, atoms):

532 # update atoms

533 self.updated = self.forces is None

534 self.set_atoms(atoms)

535 # complete energy calculations

536 if self.update_energy:

537 self.get_potential_energy(atoms)

538 # if update of forces is necessary

539 if self.update_forces:

540 # calculate forces

541 execute([self.calculate_forces])

542 # read forces

543 self.read_forces()

544

545 self.update_forces = False

546 return self.forces.copy()

547

548 def get_dipole_moment(self, atoms):

549 self.get_potential_energy(atoms)

550 self.read_dipole_moment()

551 return self.dipole

552

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

554 """return the value of a property"""

555

556 if name not in self.implemented_properties:

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

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

559 # return None

560 raise PropertyNotImplementedError(name)

561

562 if atoms is None:

563 atoms = self.atoms.copy()

564

565 persist_property = {

566 'energy': 'e_total',

567 'forces': 'forces',

568 'dipole': 'dipole',

569 'free_energy': 'e_total',

570 'charges': 'charges'

571 }

572 property_getter = {

573 'energy': self.get_potential_energy,

574 'forces': self.get_forces,

575 'dipole': self.get_dipole_moment,

576 'free_energy': self.get_potential_energy,

577 'charges': self.get_charges

578 }

579 getter_args = {

580 'energy': [atoms],

581 'forces': [atoms],

582 'dipole': [atoms],

583 'free_energy': [atoms, True],

584 'charges': [atoms]

585 }

586

587 if allow_calculation:

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

589 else:

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

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

592 else:

593 result = None

594

595 if isinstance(result, np.ndarray):

596 result = result.copy()

597 return result

598

599 def get_charges(self, atoms):

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

601 self.get_potential_energy(atoms)

602 self.read_charges()

603 return self.charges

604

605 def get_forces_on_point_charges(self):

606 """return forces acting on point charges"""

607 self.get_forces(self.atoms)

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

609 forces = []

610 for line in lines:

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

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

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

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

615

616 def set_point_charges(self, pcpot=None):

617 """write external point charges to control"""

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

619 self.pcpot = pcpot

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

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

622

623 if not self.pc_initialized:

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

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

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

627 add_data_group(

628 'point_charge_gradients',

629 'file=pc_gradients.txt'

630 )

631 drvopt = read_data_group('drvopt')

632 if 'point charges' not in drvopt:

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

634 delete_data_group('drvopt')

635 add_data_group(drvopt, raw=True)

636 self.pc_initialized = True

637

638 if self.pcpot.updated:

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

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

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

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

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

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

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

646 self.pcpot.updated = False

647

648 def read_point_charges(self):

649 """read point charges from previous calculation"""

650 charges, positions = reader.read_point_charges()

651 if len(charges) > 0:

652 self.pcpot = PointChargePotential(charges, positions)

653

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

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

656 qmmm calculations."""

657 self.pcpot = PointChargePotential(charges, positions)

658 return self.pcpot

659

660

661class PointChargePotential:

662 """Point-charge potential for Turbomole"""

663

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

665 self.mmcharges = mmcharges

666 self.mmpositions = mmpositions

667 self.mmforces = None

668 self.updated = True

669

670 def set_positions(self, mmpositions):

671 """set the positions of point charges"""

672 self.mmpositions = mmpositions

673 self.updated = True

674

675 def set_charges(self, mmcharges):

676 """set the values of point charges"""

677 self.mmcharges = mmcharges

678 self.updated = True

679

680 def get_forces(self, calc):

681 """forces acting on point charges"""

682 self.mmforces = calc.get_forces_on_point_charges()

683 return self.mmforces