XYZ to OM2

xyz2om2.py

#!/usr/bin/env python2
import numpy as np
import sys

elements = dict()
elements["H"] = 1
elements["C"] = 6
elements["N"] = 7
elements["O"] = 8
elements["F"] = 9



if __name__ == "__main__":

    xyzfile = sys.argv[1]

    if xyzfile[-4:] != ".xyz":
        print "ERROR: Not an .xyz file!", xyzfile
        exit()

    f = open(xyzfile, "r")
    lines = f.readlines()
    f.close()

    charge = 0
    if "charge = " in lines[1]:
        tokens = lines[1].split()
        charge = float(tokens[2])

    output = """OM2 1SCF MULLIK PRECISE charge=%-i

""" % charge

    natoms = int(lines[0])

    if (natoms > 3):

        for line in lines[2:2+natoms]:

            tokens = line.split()

            output += "\n%-2s %s 0 %s 0 %s 0" % (tokens[0],
                    tokens[1], tokens[2], tokens[3])

        print output


    if (natoms == 3):

        atomtypes = []
        xyz = []
        for line in lines[2:2+natoms]:

            tokens = line.split()

            atomtypes.append(tokens[0])
            xyz.append([float(tokens[1]), float(tokens[2]), float(tokens[3])])

        xyz = np.array(xyz)
        # print xyz


        ba = xyz[1] - xyz[0]
        bc = xyz[1] - xyz[2]
        cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
        angle = np.arccos(cosine_angle) / np.pi * 180.0

        print output

        print atomtypes[0]
        print atomtypes[1], np.linalg.norm(ba), 0
        print atomtypes[2], np.linalg.norm(bc), 0, angle, 0

    elif (natoms == 2):

        atomtypes = []
        xyz = []

        for line in lines[2:2+natoms]:

            tokens = line.split()

            atomtypes.append(tokens[0])
            xyz.append([float(tokens[1]), float(tokens[2]), float(tokens[3])])

        xyz = np.array(xyz)

        ba = xyz[1] - xyz[0]

        print output

        print atomtypes[0]
        print atomtypes[1], np.linalg.norm(ba), 0

    elif (natoms == 1):

        atomtypes = []
        xyz = []

        for line in lines[2:2+natoms]:

            tokens = line.split()

            atomtypes.append(tokens[0])
            xyz.append([float(tokens[1]), float(tokens[2]), float(tokens[3])])

        xyz = np.array(xyz)

        print output
        print atomtypes[0]