# -*- coding: utf-8 -*-
"""Contains a class for DFT and RI-DFT calculations using Turbomole."""
from __future__ import print_function
import os
import numpy as np
import pickle
from .iface import iface
[docs]class iface_turbomole(iface):
"""
Interface to Turbomole for DFT and RI-DFT calculations, inherits from
`metafalcon.interfaces.iface.iface`.
Parameters
----------
symbols : list of str
atomic symbols of the molecule
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
**kwargs :
See below
Keyword Arguments
-----------------
directory : str
path to the working directory of the turbomole calculation (default: tm)
take_lower : bool
whether to choose the lower energy of the two calculations in `directory` and `directory2`
(default: False)
directory2 : str
path to the second working directory, only required if `take_lower` is used
Attributes
----------
symbols : list of str
atomic symbols of the molecule
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
energy : float
electronic energy of the molecule calculated most recently (hartrees)
gradient : np.2darray (shape: (N, 3))
electronic energy gradient of the molecule with current coordinates (hartree/bohr)
ridft : bool
True, if the ri-approximation is used
See Also
--------
iface.set_symbols()
iface.set_coords()
iface.get_energy()
iface.get_gradient()
"""
def __init__(self, symbols, coords, **kwargs):
"""Construct Turbomole interface object."""
iface.__init__(self, symbols, coords)
if "directory" in kwargs.keys():
self.wd = kwargs["directory"]
else:
self.wd = "tm"
if "directory2" in kwargs.keys():
self.wd2 = kwargs["directory2"]
else:
self.wd2 = ""
if "take_lower" in kwargs.keys():
self.lower = kwargs["take_lower"]
else:
self.lower = False
self.cwd = os.getcwd()
self._check_ridft()
[docs] def calculate(self):
"""Calculate energies and gradients for the current set of coordinates."""
if not self.lower:
self.energy, self.gradient, self.converged = self._run_turbomole(self.wd)
else:
energy1, gradient1, converged1 = self._run_turbomole(self.wd)
energy2, gradient2, converged2 = self._run_turbomole(self.wd2)
if converged1 and converged2:
self.converged = True
if energy1 < energy2:
self.energy = energy1
self.gradient = gradient1
else:
self.energy = energy2
self.gradient = gradient2
else:
self.converged = False
self.energy = 0.0
self.gradient = np.zeros_like(self.coords)
print(self.converged, self.energy)
def _run_turbomole(self, wd):
"""Run turbomole calulation in the given working directory."""
os.chdir(wd)
self._write_input()
if self.ridft:
os.system("ridft > tm.out; rdgrad > rdgrad.out")
else:
os.system("dscf > tm.out; grad > grad.out")
try:
gradient = self._read_gradient()
energy = self._read_energy()
converged = True
except IOError as e:
energy = 0.0
gradient = np.zeros_like(self.coords)
converged = False
print(e)
with open("energies.pickle", "w") as f:
pickle.dump(energy, f, protocol=pickle.HIGHEST_PROTOCOL)
with open("gradients.pickle", "w") as f:
pickle.dump(gradient, f, protocol=pickle.HIGHEST_PROTOCOL)
os.chdir(self.cwd)
return energy, gradient, converged
def _check_ridft(self):
self.ridft = False
with open(self.wd+"/control") as f:
for line in f:
if "$rij" in line:
self.ridft = True
break
def _write_input(self):
fmt = "%20.12f %20.12f %20.12f %s \n"
with open("coord", "w") as f:
f.write("$coord\n")
for sym, coord in zip(self.symbols, self.coords):
f.write(fmt % (coord[0], coord[1], coord[2], sym))
f.write("$end")
def _read_gradient(self):
gradient = []
nat = self.coords.shape[0]
with open("gradient") as f:
for i, line in enumerate(f):
if i > nat+1 and len(gradient) < nat*3:
gradient += line.split()
os.system("rm gradient")
gradient = [grad.replace("D", "E") for grad in gradient]
return np.array(gradient, dtype=float).reshape(self.coords.shape)
def _read_energy(self):
with open("energy") as f:
energies = [line.split()[1] for line in f.readlines() if "$" not in line]
energy = float(energies[-1])
os.system("rm energy")
return energy