# -*- coding: utf-8 -*-
"""Contains a class for semiempirical calculations using Mopac."""
import os
import numpy as np
from .. import xyz
from .iface import iface
[docs]class iface_dftbplus(iface):
"""
Interface to DFTB+ for semiempirical tight-bindind 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
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)
See Also
--------
iface.set_symbols()
iface.set_coords()
iface.get_energy()
iface.get_gradient()
"""
def __init__(self, symbols, coords, **kwargs):
"""Construct DFTB+ interface object."""
iface.__init__(self, symbols, coords)
[docs] def initialCalculation(self):
"""Calulates everything the first time from the scratch to construct partial charges
as an initial guess for the next steps"""
with open("dftb_in.hsd", "r") as f:
lines = f.readlines()
for idx, line in enumerate(lines):
if "ReadInitialCharges" in line:
break
lines[idx] = "ReadInitialCharges = No \n"
with open("dftb_in.hsd", "w") as f:
f.writelines(lines)
os.system("dftb+")
lines[idx] = "ReadInitialCharges = Yes \n"
with open("dftb_in.hsd", "w") as f:
f.writelines(lines)
return None
[docs] def calculate(self):
"""Calculate energies and gradients for the current set of coordinates."""
xyz.write_xyz("input.xyz", self.symbols, self.coords)
os.system("xyz2gen input.xyz")
if not os.path.exists("charges.bin"):
self.initialCalculation()
os.system("dftb+")
self.gradient = []
with open("results.tag") as f:
lines = f.readlines()
for iline, line in enumerate(lines):
if "forces" in line:
index_force_begin = iline + 1
line1 = line.replace(':', ',')
index_force_end = iline + 1 + int(line1.split(',')[-1])
break
for j in range(index_force_begin, index_force_end):
word = lines[j].split()
self.gradient.append([float(word[k]) for k in range(0, 3)])
self.gradient = -np.array(self.gradient)
with open("detailed.out") as f:
lines = f.readlines()
for line in lines:
if line.strip().startswith('Dipole moment:'):
dipole = np.array([float(line.split()[2]), float(line.split()[3]), float(line.split()[4])])
break
for line in lines:
if line.strip().startswith('Total energy:'):
self.energy = float(line.split()[2])
break
with open("dipole.dat", "a") as f:
f.write("%12.6f %12.6f %12.6f \n" % (dipole[0], dipole[1], dipole[2]))
self.converged = True
if not len(self.gradient):
self.energy = 0.0
self.gradient = np.zeros_like(self.coords)
self.converged = False
os.remove('results.tag')