Wrapper to get energy and gradients from MNDO99

mndo99_wrapper.py

from __future__ import print_function
import os
import numpy as np

__EXECUTABLE__ = "/home/andersx/opt/mndo/mndo99_20121112_intel64_ifort-11.1.080_mkl-10.3.12.361"


def write_om2_gradient(coords, species, filename):
    """ Writes an input file to calculate the
        OM2 energy and gradient with MNDO99.

        Coords: Numpy array of coordinates in angstrom
        Species: list of names of each elements
        filename: Filename of the inputfile
    """

    content = "OM2 JOP=-2 PRECISE\n\n"

    for i in range(len(species)):

        content += "\n%-2s %20.10f 0 %20.10f 0 %20.10f 0 " % \
            (species[i], coords[i,0], coords[i,1], coords[i,2])

    f = open(filename, "w")
    f.write(content)
    f.close()

def read_om2_gradient(filename):
    """ Reads the output file from a energy+gradient
        calculation and returns the energy in [kcal/mol]
        and gradient in [kcal/(mol*angstrom)]

        filename: Filename of the outputfile
    """

    f = open(filename)
    lines = f.readlines()
    f.close()
    start = -1

    energy = 1e20

    for i, line in enumerate(lines):
        if "SCF HEAT OF FORMATION" in line:
            tokens = line.split()
            energy = float(tokens[4])
            break

    if (energy > 1e19):
        print("QM calculation failed with bad energy")
        exit()


    for i, line in enumerate(lines):
        if "GRADIENTS" in line:
            start = i + 4
            break

    if start == -1:
        print("QM calculation failed without gradient")
        exit()

    grad = []

    for line in lines[start:]:
        tokens = line.split()
        if len(tokens) != 8:
            break

        x = float(tokens[5])
        y = float(tokens[6])
        z = float(tokens[7])

        grad.append([x,y,z])

    grad = np.array(grad)
    return energy, grad


def om2_gradient(coords, species):
    """ Calculates the OM2 energy and gradient using MNDO99.

        Coords: Numpy array of coordinates in angstrom
        Species: list of names of each elements
    """

    inp_filename = "temp.inp"
    out_filename = "temp.log"

    write_om2_gradient(coords, species, inp_filename)

    cmd = "%s  < %s > %s" % (__EXECUTABLE__, inp_filename, out_filename)
    os.system(cmd)

    energy, gradient = read_om2_gradient(out_filename)

    os.remove(inp_filename)
    os.remove(out_filename)

    return energy, gradient