# -*- coding: utf-8 -*-
"""Contains a class for DFT, AM1 or CASSCF calculations using Gaussian09."""
import os
import numpy as np
from .iface import iface
from .. import constants as c
[docs]class iface_gaussian(iface):
"""
Interface to Gaussian09 for DFT, AM1 or CASSCF 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
-----------------
none
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)
head : list of str
lines to be used as head for the g09 inputfile
tail : list of str
lines to be used as tail for the g09 inputfile
head1 : list of str
lines to be used as head for the optional excited state g09 inputfile
tail1 : list of str
lines to be used as tail for the optional excited state g09 inputfile
cas : bool
whether to do CASSCF calculation
See Also
--------
iface.set_symbols()
iface.set_coords()
iface.get_energy()
iface.get_gradient()
"""
def __init__(self, symbols, coords, **kwargs):
"""Construct Gaussian09 interface object."""
iface.__init__(self, symbols, coords)
# read in head and tail
if os.path.exists("0/head"):
with open("0/head") as f:
self.head = f.readlines()
with open("0/tail") as f:
self.tail = f.readlines()
with open("1/head") as f:
self.head1 = f.readlines()
with open("1/tail") as f:
self.tail1 = f.readlines()
else:
with open("head") as f:
self.head = f.readlines()
with open("tail") as f:
self.tail = f.readlines()
self.append_to_head("fchk")
self.append_to_head("nosymm")
self.append_to_head("force")
if "cas" in " ".join(self.head).lower():
self.cas = True
if "#p" not in " ".join(self.head).lower():
self.head = [line.replace("#", "#p") for line in self.head]
else:
self.cas = False
if os.path.exists("inp.chk"):
os.remove("inp.chk")
if os.path.exists("0/inp.chk"):
os.remove("0/inp.chk")
if os.path.exists("1/inp.chk"):
os.remove("1/inp.chk")
[docs] def append_to_head(self, task):
for i, line in enumerate(self.head):
if "#" in line and task not in line:
self.head[i] = line + " " + task + "\n"
[docs] def calculate(self):
"""Calculate energies and gradients for the current set of coordinates."""
try:
if os.path.exists("0/head"):
fchkpath = "0/Test.FChk"
outpath = "0/out"
if os.path.exists(fchkpath):
os.system("rm " + fchkpath)
if os.path.exists("0/inp.chk") and os.path.exists("1/inp.chk"):
for i, line in enumerate(self.head):
if "#" in line and "guess" not in line:
self.head[i] = line + " guess=read\n"
for i, line in enumerate(self.head1):
if "#" in line and "guess" not in line:
self.head1[i] = line + " guess=read\n"
self._write_input("0/inp", self.head, self.tail)
self._write_input("1/inp", self.head1, self.tail1)
run_g09(path=str(0))
os.system("cp 0/inp.chk 1/inp.chk")
run_g09(path=str(1))
else:
fchkpath = "Test.FChk"
outpath = "out"
if os.path.exists(fchkpath):
os.system("rm " + fchkpath)
if os.path.exists("inp.chk"):
for i, line in enumerate(self.head):
if "#" in line and "guess" not in line:
self.head[i] = line + " guess=read\n"
self._write_input("inp", self.head, self.tail)
run_g09()
gradient = []
with open(fchkpath) as f:
for line in f:
if "Cartesian Gradient" in line:
while len(gradient) < np.prod(self.coords.shape):
gradient += next(f).split()
break
self.gradient = np.array(gradient,dtype=float).reshape(self.coords.shape)
with open(outpath) as f:
for line in f:
if "SCF Done:" in line:
self.energy = float(line.split()[4])
break
elif self.cas and "EIGENVALUE " in line:
tmp = line.split()
if tmp[1] == "1)":
self.energy = float(tmp[-1])
self.converged = True
except IOError as e:
self.energy = 0.0
self.gradient = np.zeros_like(self.coords)
self.converged = False
print(e)
def _write_input(self, path, head, tail):
fmt = "%s %20.12f %20.12f %20.12f \n"
with open(path, "w") as f:
f.writelines(head)
for sym, coord in zip(self.symbols, self.coords*c.a0):
f.write(fmt % (sym, coord[0], coord[1], coord[2]))
f.write("\n")
f.writelines(tail)
[docs]def run_g09(path="."):
"""
Execute Gaussian09 in order to get outputfile *out* from inputfile *inp*.
Parameters
----------
path : str
path to the calculation directory to be used (default: .)
"""
os.system("cd %s; g09 < inp > out" % path)