# -*- coding: utf-8 -*-
"""Contains a class for semiempirical calculations using Mopac."""
import os
import numpy as np
from .iface import iface
from .. import constants as c
[docs]class iface_mopac(iface):
"""
Interface to Mopac for semiempirical 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
-----------------
keywords : str
keywords to be written to the mopac inputfile additional to 1SCF GRAD
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)
keywords : str
keywords to be written to the mopac inputfile additional to 1SCF GRAD
See Also
--------
iface.set_symbols()
iface.set_coords()
iface.get_energy()
iface.get_gradient()
"""
def __init__(self, symbols, coords, **kwargs):
"""Construct Mopac interface object."""
iface.__init__(self, symbols, coords)
self.keywords = kwargs["keywords"]
[docs] def calculate(self):
"""Calculate energies and gradients for the current set of coordinates."""
with open("mopac.mop", "w") as f:
f.write(self.keywords.upper() + " 1SCF GRAD \n\n\n")
for sym, coord in zip(self.symbols, self.coords*c.a0):
f.write("%s %20.12f %20.12f %20.12f \n" %(sym, coord[0], coord[1], coord[2]))
os.system("mopac mopac.mop")
self.gradient = []
with open("mopac.out") as f:
for line in f:
if "FINAL HEAT OF FORMATION" in line:
self.energy = float(line.split()[5])*c.kcalmol2au
if "FINAL POINT AND DERIVATIVES" in line:
f.next()
f.next()
if "mozyme" in self.keywords.lower():
f.next()
f.next()
grad_list = [f.next().split()[2:5] for coord in self.coords]
self.gradient = np.array(grad_list, dtype=float)
else:
grad_list = [float(f.next().split()[6]) for i in
range(np.prod(self.coords.shape))]
self.gradient = np.array(grad_list).reshape(self.coords.shape)
self.gradient *= c.kcalmol2au*c.a0
self.converged = True
if not len(self.gradient):
self.energy = 0.0
self.gradient = np.zeros_like(self.coords)
self.converged = False