#!/usr/bin/env python
"""
Handles importing data from the various filetypes that Q2MM uses.

Schrodinger
-----------
When importing Schrodinger files, if the atom.typ file isn't in the directory
where you execute the Q2MM Python scripts, you may see this warning:

  WARNING mmat_get_atomic_num x is not a valid atom type
  WARNING mmat_get_mmod_name x is not a valid atom type

In this example, x is the number of a custom atom type defined and added to
atom.typ. The warning can be ignored. If it's bothersome, copy atom.typ into
the directory where you execute the Q2MM Python scripts.

Note that the atom.typ must be located with your structure files, else the
Schrodinger jobs will fail.
"""
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division

from argparse import RawTextHelpFormatter
from string import digits
import logging
import mmap
import numpy as np
import math
import os
import re
import subprocess as sp
import time
import sys

try:
    from schrodinger import structure as sch_str
    from schrodinger.application.jaguar import input as jag_in
except:
    print("Schrodinger not installed, limited functionality")
    pass

import constants as co
import datatypes

logger = logging.getLogger(__name__)
# Print out full matrices rather than having Numpy truncate them.
# np.nan seems to no longer be supported for untruncated printing
# of arrays. The suggestion is to use sys.maxsize but I haven't checked
# that this works for python2 so leaving the commented code for now.
# np.set_printoptions(threshold=np.nan)
np.set_printoptions(threshold=sys.maxsize)

class File(object):
    """
    Base for every other filetype class.
    """
    def __init__(self, path):
        self._lines = None
        self.path = os.path.abspath(path)
        # self.path = path
        self.directory = os.path.dirname(self.path)
        self.filename = os.path.basename(self.path)
        # self.name = os.path.splitext(self.filename)[0]
    @property
    def lines(self):
        if self._lines is None:
            with open(self.path, 'r') as f:
                self._lines = f.readlines()
        return self._lines
    def write(self, path, lines=None):
        if lines is None:
            lines = self.lines
        with open(path, 'w') as f:
            for line in lines:
                f.write(line)


# Currently only for 1 system.
# 
class AmberHess(File):
    def __init__(self, path):
        super(AmberHess, self).__init__(path)
        self._hessian = None
        self.natoms = None
    @property
    def hessian(self):
        if self._hessian is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            with open("./calc/"+self.filename, 'r') as f:
                lines = f.readlines()
            for i,line in enumerate(lines):
                if i == 0:
                    self.natoms = int(line.split()[1])
                    hessian = np.zeros([self.natoms * 3, self.natoms * 3], dtype=float)
                else:
                    row = np.array(line.split()).astype(np.float)
                    hessian[:,i-1] = row
            # Convert hessian units to use kJ/mol instead of kcal/mol.

            # kcal/mol for energy in AMBER
            # E(kcal/mol -> cm**-1) = 349.75
            # freq = sqrt(lambda(kcal/mol)) / (2 pi c)
            
            w, v = np.linalg.eigh(hessian)
            eigval = np.zeros([self.natoms * 3],dtype=float)
            for i,eig in enumerate(w):
                if eig < 0:
                    eigval[i] = -np.sqrt(-eig)
                else:
                    eigval[i] = np.sqrt(eig)
            eigval *= 108.587 # freq in cm**-1
            self._hessian = hessian / co.HARTREE_TO_KCALMOL \
                * co.HARTREE_TO_KJMOL
            logger.log(5, '  -- Finished Creating {} Hessian matrix.'.format(
                hessian.shape))
            return self._hessian
class AmberEne(File):
    """
        Amber .ene file to read either current energy or optimized energy
    """
    def __init__(self, path):
        super(AmberEne, self).__init__(path)
        self._structures = None
        self.name = None
    @property
    def structures(self):
        if self._structures == None:
            logger.log(10, 'READING: {}'.format(self.filename))
            self._structures = []
            flag = 0
            with open('./calc/'+self.filename, 'r') as f:
                sections = {'sp':1, 'minimization':2}
                calc_section = 'sp'
                count_previous = 0
                    
                for line in f:
                    count_current = sections[calc_section]
                    if count_current != count_previous:
                        current_structure = Structure()
                        self._structures.append(current_structure)
                        count_previous += 1
                    if 'FINAL RESULTS' in line:
                        flag = 1
                    elif flag == 1 and "NSTEP" in line:
                        flag = 2
                    elif flag == 2:
                        energy = self.read_line_for_energy(line)
                        if energy is not None:
                            current_structure.props['energy']=energy
                        flag = 0
            logger.log(5, '  -- Imported {} structure(s)'.format(
                len(self._structures)))
        return self._structures


    def read_line_for_energy(self, line):
        # The Amber Energy is in units of kcal/mol, so we have to convert them to kJ/mol
        # for consistency purposes.
        # don't know how to use match = re.compile 
        linesplit = line.split()
        energy = float(linesplit[1])
        energy *= co.HARTREE_TO_KJMOL / co.HARTREE_TO_KCALMOL
        return energy
#        if match:
#            energy = float(match.group(1))
#            energy *= co.HARTREE_TO_KJMOL / co.HARTREE_TO_KCALMOL
#            return energy
#        else:
#            return None

class AmberGeo(File):
    """
        .geo file to be used for bond,angles,dihedral sets
        .out file for the current value
    """
    def __init__(self, path):
        super(AmberGeo, self).__init__(path)
        self._structures = None
        self.name = None
    @property
    def structures(self):
        if self._structures == None:
            logger.log(10, 'READING: {}'.format(self.filename))
            self._structures = []
            with open("./calc/"+self.filename, 'r') as f:
                sections = {'sp':1, 'minimization':2, 'hessian':2}
                count_previous = 0
                calc_section = 'sp'
                b = 0
                a = 0  
                t = 0
                for line in f:
                    count_current = sections[calc_section]
                    if count_current != count_previous:
                        bonds = []
                        angles = []
                        torsions = []
                        current_structure = Structure()
                        self._structures.append(current_structure)
                        count_previous += 1
                    section = None
                    if "END" in line:
                        t = 0
                        calc_section = 'minimization'
                        for bond in bonds:
                            bond.atom_nums.sort()
                        bonds.sort(key=lambda x: (x.atom_nums[0],
                                                  x.atom_nums[1]))
                        for angle in angles:
                            if angle.atom_nums[0] > angle.atom_nums[2]:
                                angle.atom_nums = [angle.atom_nums[2],
                                                   angle.atom_nums[1],
                                                   angle.atom_nums[0]]
                        for torsion in torsions:
                            if torsion.atom_nums[1] > torsion.atom_nums[2]:
                                torsion.atom_nums = [torsion.atom_nums[3],
                                                     torsion.atom_nums[2],
                                                     torsion.atom_nums[1],
                                                     torsion.atom_nums[0]]
                        angles.sort(key=lambda x: (x.atom_nums[1],
                                                   x.atom_nums[0],
                                                   x.atom_nums[2]))
                        torsions.sort(key=lambda x: (x.atom_nums[1],
                                                     x.atom_nums[2],
                                                     x.atom_nums[0],
                                                     x.atom_nums[3]))
                        current_structure.bonds.extend(bonds)
                        current_structure.angles.extend(angles)
                        current_structure.torsions.extend(torsions)

                    if t == 1:
                        torsion = self.read_line_for_torsion(line)
                        if torsion is not None:
                            torsions.append(torsion)
                    elif 'TORSIONS' in line:
                        t = 1
                        a = 0
                    if a == 1:
                        angle = self.read_line_for_angle(line)
                        if angle is not None:
                            angles.append(angle)
                    elif 'ANGLES' in line:
                        a = 1
                        b = 0
                    if b == 1:
                        bond = self.read_line_for_bond(line)
                        if bond is not None:
                            bonds.append(bond)
                    elif 'BONDS' in line:
                        b = 1


            logger.log(5, '  -- Imported {} structure(s)'.format(
                len(self._structures)))
        return self._structures

    def read_line_for_bond(self, line):
        # All bond data starts with the string "Bond" and then the rest of the
        # interaction information.
        a,b,z = line.split()
        atom_nums = [int(x) for x in [a,b]]
        value = float(z)
        return Bond(atom_nums=atom_nums, value=value)

    def read_line_for_angle(self, line):
        a,b,c,z = line.split()
        atom_nums = [int(x) for x in [a,b,c]]
        value = float(z)
        return Angle(atom_nums=atom_nums, value=value)

    def read_line_for_torsion(self, line):
        a,b,c,d,z = line.split()
        atom_nums = [int(x) for x in [a,b,c,d]]
        value = float(z)
        return Torsion(atom_nums=atom_nums, value=value)

    def read_line_for_energy(self, line):
        # The TPE is in units of kcal/mol, so we have to convert them to kJ/mol
        # for consistency purposes.
        match = re.compile('Total Potential Energy :\s+({0})'.format(
            co.RE_FLOAT)).search(line)
        if match:
            energy = float(match.group(1))
            energy *= co.HARTREE_TO_KJMOL / co.HARTREE_TO_KCALMOL
            return energy
        else:
            return None
class AmberLeap_Gaus(File):
    def __init__(self, path):
        """
            run -> gaus to amber -> sp -> traj -> cpptraj -> cpptraj -> AmberGeo
            path = leap.in
        """
        super(AmberLeap_Gaus, self).__init__(path)
        self._index_output_log = None
        self._structures = None
        self.commands = None
        self.name = os.path.splitext(self.filename)[0]
        self.filename = self.name + '.in' # .log file to .in (.in file is never replaced. so using .in should have original coordinate)
        self.name_log = 'gaus.' + self.name + '.log'
        self.name_prm = 'gaus.' + self.name + '.parm7' #topology
        self.name_rst = 'gaus.' + self.name + '.rst7' # coordinate
        self.name_min = 'gaus.' + self.name + '.min' # sander min input
        self.name_ene = 'gaus.' + self.name + '.ene'
        self.name_dyn = 'gaus.' + self.name + '.dyn' # sander dyn input
        self.name_int = 'gaus.' + self.name + '.int' # interaction input for cpptraj
        self.name_geo = 'gaus.' + self.name + '.geo' # cpptraj output for all interactions (to be read by AmberGeo)
        self.min_script = """Comments
 &cntrl
  imin      = 1,
  ntx       = 1,
  maxcyc    = 0,
  ncyc      = 0,
  cut       = 15.0,
  ntpr      = 10000,
  ntwx      = 0,
  ntb       = 0
 /
"""
        self.dyn_script = """Comments
 &cntrl
  imin      = 0,
  ntx       = 1,
  irest     = 0,
  nstlim    = 0,
  ntwx      = 1,
  cut       = 15.0,
  ntb       = 0,
  ntpr      = 1
 /
"""
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            struct = Structure()
            self._structures = [struct]
            with open(self.filename, 'r') as f:
                for line in f:
                    line = line.split()
                    if len(line) == 2:
                        struct.props['total atoms'] = int(line[0])
                        struct.props['title'] = line[1]
                        logger.log(5, '  -- Read {} atoms.'.format(
                            struct.props['total atoms']))
                    if len(line) > 2:
                        indx, ele, x, y, z, at, bonded_atom = line[0], \
                            line[1], line[2], line[3], line[4], \
                            line[5], line[6:]
                        struct.atoms.append(Atom(index=indx,
                            element=ele,
                            x=float(x),
                            y=float(y),
                            z=float(z),
                            atom_type=at,
                            atom_type_name=at,
                            bonded_atom_indices=bonded_atom))
            return self._structures
    def get_com_opts(self):
        com_opts = {'freq': False,
                    'opt': False,
                    'sp':True,
                    'tors': False,
                    'geo':True}
        return com_opts

    def extract(self,log):
        script="""
trajin calc/gaus.NAME.nc
fixatomorder
AA
run
write
exit
"""
        script = script.replace("NAME",self.name)
        geo = ""

        # read .geo file and store all possible interaction
        bonds = []
        angles = []
        torsions = []
        ref = open('./calc/'+self.name_geo,'r').readlines()
        count = 0
        for line in ref:
            # Bonds
            if "[angles]" in line:
                count = 0
            elif count == 1:
                bonds.append(line.split()[-4:-2])
            elif "Atom2" in line:
                count = 1

            # Angles
            if "[dihedrals]" in line:
                count = 0
            if count == 2:
                angles.append(line.split()[-6:-3])
            elif "Atom3" in line:
                count = 2

            # Dihedral
            # store the columns as negtive since there is unexpected "B" or E in front of column
            if "TIME" in line:
                count = 0
            if count == 3:
                torsions.append(line.split()[-8:-4])
            elif "Atom4" in line:
                count = 3
        
        for a,b in bonds:
            geo += "distance @{} @{} out calc/gaus.bonds".format(a,b) + '\n'
        for a,b,c in angles:
            geo += "angle @{} @{} @{} out calc/gaus.angles".format(a,b,c) + '\n'
        for a,b,c,d in torsions:
            geo += "dihedral @{} @{} @{} @{} out calc/gaus.torsions".format(a,b,c,d) + '\n'
        
        script = script.replace("AA",geo)
        script_f = './calc/' + self.name + '.temp'
        with open(script_f, 'w') as f:
            f.write(script)
        sp.call("cpptraj -p calc/prmtop < {}".format(script_f), shell=True, stderr=log, stdin=log, stdout=log)
        summary = ""
        if os.path.isfile("calc/gaus.bonds"):
            bond_file = open("calc/gaus.bonds","r").readlines()
            bond_line = bond_file[-1].split()[1:]
            summary += "BONDS\n"
            i = 0
            for a,b in bonds:
                summary += "{} {} {} \n".format(a,b,bond_line[i])
                i += 1
        if os.path.isfile("calc/gaus.angles"):
            angle_file = open("calc/gaus.angles","r").readlines()
            angle_line = angle_file[-1].split()[1:]
            summary += "ANGLES\n"
            i = 0
            for a,b,c in angles:
                summary += "{} {} {} {} \n".format(a,b,c,angle_line[i])
                i += 1
        if os.path.isfile("calc/gaus.torsions"):
            tors_file = open("calc/gaus.torsions","r").readlines()
            tors_line = tors_file[-1].split()[1:]
            summary += "TORSIONS\n"
            i = 0
            for a,b,c,d in torsions:
                summary += "{} {} {} {} {} \n".format(a,b,c,d,tors_line[i])
                i += 1
        summary += "END"
        # replace name_geo with summary
        with open('./calc/'+self.name_geo,'w') as f:
            f.write(summary)
        return

    def geometry(self,log):
        # Run Trajectory (Required for cpptraj)
        with open("./calc/"+self.name_dyn, 'w') as f:
            f.write(self.dyn_script)
        sp.call("sander -O -i calc/{} -o calc/traj.out -p calc/prmtop -c calc/gaus.{}.rst -x calc/gaus.{}.nc".format(self.name_dyn,self.name,self.name),shell=True)
        # Generate All geometry
        int_script = "bonds\nangles\ndihedrals\n"
        with open('./calc/'+self.name_int, 'w') as f:
            f.write(int_script)
        sp.call("cpptraj -p calc/prmtop < calc/{} > calc/{} \n".format(self.name_int,self.name_geo),shell=True)
        self.extract(log)
        return
    def run(self,check_tokens=False):
        logger.log(5, 'RUNNING: {}'.format(self.filename))
        self._index_output_log = []
        com_opts = self.get_com_opts()
        current_directory = os.getcwd()
        os.chdir(self.directory)
        log = open(self.name_log,'w')
        os.chdir(self.directory)
        if os.path.isfile('calc'):
            os.remove('calc')
        sp.call("mkdir calc",shell=True, stderr=log, stdin=log, stdout=log)
        if com_opts['sp']:
            logger.log(1, '  CALCULATE: {}'.format(self.filename))
            # Run leap
            sp.call("tleap -f {}".format(self.filename),shell=True, stderr=log, stdin=log, stdout=log) # parm7 rst7 files made
            # Run Min
            with open("./calc/"+self.name_min, 'w') as f:
                f.write(self.min_script)
            sp.call("sander -O -i calc/{} -o calc/{} -p calc/prmtop -c calc/inpcrd -r calc/gaus.{}.rst".format(self.name_min,self.name_ene,self.name),shell=True, stderr=log, stdin=log, stdout=log)
        if com_opts['geo']:
            self.geometry(log)
        os.chdir(current_directory)

class AmberLeap(File):
    def __init__(self, path):
        """
            path = leap.in
        """
        super(AmberLeap, self).__init__(path)
        self._index_output_log = None
        self._structures = None
        self.commands = None
        self.name = os.path.splitext(self.filename)[0]
        self.name_log = 'amber.' + self.name + '.log'
        self.name_prm = 'amber.' + self.name + '.parm7' #topology
        self.name_rst = 'amber.' + self.name + '.rst7' # coordinate
        self.name_min = 'amber.' + self.name + '.min' # sander min input
        self.name_ene = 'amber.' + self.name + '.ene'
        self.name_dyn = 'amber.' + self.name + '.dyn' # sander dyn input
        self.name_int = 'amber.' + self.name + '.int' # interaction input for cpptraj
        self.name_geo = 'amber.' + self.name + '.geo' # cpptraj output for all interactions
        self.name_hes = 'amber.' + self.name + '.hes'
        self.geo = None
        self.min_script = """Comments
 &cntrl
  imin      = 1,
  ntx       = 1,
  maxcyc    = aa,
  ncyc      = bb,
  cut       = 15.0,
  ntpr      = 10000,
  ntwx      = 0,
  ntb       = 0
 /
"""
        self.dyn_script = """Comments
 &cntrl
  imin      = 0,
  ntx       = 1,
  irest     = 0,
  nstlim    = 0,
  ntwx      = 1,
  cut       = 15.0,
  ntb       = 0,
  ntpr      = 1
 /
"""
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            struct = Structure()
            self._structures = [struct]
            with open(self.filename, 'r') as f:
                for line in f:
                    line = line.split()
                    if len(line) == 2:
                        struct.props['total atoms'] = int(line[0])
                        struct.props['title'] = line[1]
                        logger.log(5, '  -- Read {} atoms.'.format(
                            struct.props['total atoms']))
                    if len(line) > 2:
                        indx, ele, x, y, z, at, bonded_atom = line[0], \
                            line[1], line[2], line[3], line[4], \
                            line[5], line[6:]
                        struct.atoms.append(Atom(index=indx,
                            element=ele,
                            x=float(x),
                            y=float(y),
                            z=float(z),
                            atom_type=at,
                            atom_type_name=at,
                            bonded_atom_indices=bonded_atom))
            return self._structures
    def get_com_opts(self):
        com_opts = {'freq': False,
                    'opt': False,
                    'sp': False,
                    'tors': False,
                    'geo':False}
        if any(x in ['ab','aa','at','abo','aao','ato'] for x in self.commands):
            com_opts['geo'] = True
        if any(x in ['abo','aao','ato','aeo','ae1o','aeao'] for x in self.commands):
            com_opts['opt'] = True
            com_opts['sp'] = True
        if any(x in ['ah', 'ajeig', 'ageig'] for x in self.commands):
            com_opts['geo'] = True
            com_opts['freq'] = True
            com_opts['opt'] = True
            com_opts['sp'] = True
        if any(x in ['at', 'ato'] for x in self.commands):
            com_opts['tors'] = True
        return com_opts
    def extract(self,log):
        script="""
trajin calc/amber.NAME.nc
fixatomorder
AA
run
write
exit
"""
        script = script.replace("NAME",self.name)
        geo = ""

        # read .geo file and store all possible interaction
        bonds = []
        angles = []
        torsions = []
        ref = open('./calc/'+self.name_geo,'r').readlines()
        self.geo = ref
        count = 0
        for line in ref:
            # Bonds
            if "[angles]" in line:
                count = 0
            elif count == 1:
                bonds.append(line.split()[-4:-2])
            elif "Atom2" in line:
                count = 1

            # Angles
            if "[dihedrals]" in line:
                count = 0
            if count == 2:
                angles.append(line.split()[-6:-3])
            elif "Atom3" in line:
                count = 2

            # Dihedral
            # store the columns as negtive since there is unexpected "B" or E in front of column
            if "TIME" in line:
                count = 0
            if count == 3:
                torsions.append(line.split()[-8:-4])
            elif "Atom4" in line:
                count = 3

        for a,b in bonds:
            geo += "distance @{} @{} out calc/amber.bonds".format(a,b) + '\n'
        for a,b,c in angles:
            geo += "angle @{} @{} @{} out calc/amber.angles".format(a,b,c) + '\n'
        for a,b,c,d in torsions:
            geo += "dihedral @{} @{} @{} @{} out calc/amber.torsions".format(a,b,c,d) + '\n'
        script = script.replace("AA",geo)
        script_f = './calc/' + self.name + '.temp'
        with open(script_f, 'w') as f:
            f.write(script)
        sp.call("cpptraj -p calc/prmtop < {}".format(script_f), shell=True, stderr=log, stdin=log, stdout=log)
        summary = ""
        if os.path.isfile("calc/amber.bonds"):
            bond_file = open("calc/amber.bonds","r").readlines()
            bond_line = bond_file[-1].split()[1:]
            summary += "BONDS\n"
            i = 0
            for a,b in bonds:
                summary += "{} {} {} \n".format(a,b,bond_line[i])
                i += 1
        if os.path.isfile("calc/amber.angles"):
            angle_file = open("calc/amber.angles","r").readlines()
            angle_line = angle_file[-1].split()[1:]
            summary += "ANGLES\n"
            i = 0
            for a,b,c in angles:
                summary += "{} {} {} {} \n".format(a,b,c,angle_line[i])
                i += 1
        if os.path.isfile("calc/amber.torsions"):
            tors_file = open("calc/amber.torsions","r").readlines()
            tors_line = tors_file[-1].split()[1:]
            summary += "TORSIONS\n"
            i = 0
            for a,b,c,d in torsions:
                summary += "{} {} {} {} {} \n".format(a,b,c,d,tors_line[i])
                i += 1
        summary += "END"
        # replace name_geo with summary
        with open('./calc/'+self.name_geo,'w') as f:
            f.write(summary)
        return

    def hessian(self,log):
        # if pdb file does not exit, then convert mol2 to pdb
        os.chdir(self.directory)
        if os.path.isfile(self.name+".pdb"):
            0
        else:
            sp.call("antechamber -dr no -i {} -fi mol2 -o {} -fo pdb".format(self.name+".mol2",self.name+".pdb"),shell=True)
        # nab input file
        # dielectric constant = 80.4 for water.
        # currently manual change required
        script = """molecule m;
float x[4000], fret;

m = getpdb("{}.pdb");
readparm(m, "./calc/prmtop");

mm_options( "cut=15., ntpr=1, nsnb=99999, diel = C, dielc = 80.40" );
mme_init( m, NULL, "::Z", x, NULL);
setxyz_from_mol( m, NULL, x );

nmode( x, 3*m.natoms, mme2, 0, 0, 0.0, 0.0, 0);""".format(self.name)
        with open('./calc/'+self.name+'.nab','w') as f:
            f.write(script)
        # nab compile
        sp.call("nab calc/{}.nab".format(self.name),shell=True)
        # nab run
        sp.call("./calc/{}".format(self.name),shell=True,stderr=log, stdin=log, stdout=log)
        # hessian.mat formed
        # rename to .hess
        sp.call("mv ./calc/hessian.mat ./calc/{}".format(self.name_hes),shell = True)
        return
    def geo_extract(self):
    
        bonds = []
        angles = []
        torsions = []
    
        ref = self.geo
        count = 0
        for line in ref:
            # Bonds
            if "[angles]" in line:
                count = 0
            elif count == 1:
                bonds.append(line.split()[-4:-2])
            elif "Atom2" in line:
                count = 1

            # Angles
            if "[dihedrals]" in line:
                count = 0
            if count == 2:
                angles.append(line.split()[-6:-3])
            elif "Atom3" in line:
                count = 2

            # Dihedral
            # store the columns as negtive since there is unexpected "B" or E in front of column
            if "TIME" in line:
                count = 0
            if count == 3:
                torsions.append(line.split()[-8:-4])
            elif "Atom4" in line:
                count = 3

        hes_ele = np.array([None,None,None,None])
        for a,b in bonds:
            hes_ele = np.vstack((hes_ele,[a,b,None,None]))
        for a,b,c in angles:
            hes_ele = np.vstack((hes_ele,[a,b,c,None]))
        for a,b,c,d in torsions:
            hes_ele = np.vstack((hes_ele,[a,b,c,d]))
        np.save("calc/geo",hes_ele)
        return
        
    def geometry(self,log):
        # Run Trajectory (Required for cpptraj)
        with open("./calc/"+self.name_dyn, 'w') as f:
            f.write(self.dyn_script)
        sp.call("sander -O -i calc/{} -o calc/traj.out -p calc/prmtop -c calc/amber.{}.rst -x calc/amber.{}.nc".format(self.name_dyn,self.name,self.name),shell=True)
        # Generate All geometry
        int_script = "bonds\nangles\ndihedrals\n"
        with open('./calc/'+self.name_int, 'w') as f:
            f.write(int_script)
        sp.call("cpptraj -p calc/prmtop < calc/{} > calc/{}".format(self.name_int,self.name_geo),shell=True)
        self.extract(log)
        
        return
    def run(self,check_tokens=False):
        logger.log(5, 'RUNNING: {}'.format(self.filename))
        self._index_output_log = []
        com_opts = self.get_com_opts()
        current_directory = os.getcwd()
        os.chdir(self.directory)
        log = open(self.name_log,'w')
        sp.call("mkdir calc",shell=True, stderr=log, stdin=log, stdout=log)
        if com_opts['opt']:
            logger.log(1, '  MINIMIZE & ANALYZE: {}'.format(self.filename))
            # Run leap
            sp.call("tleap -f {}".format(self.filename),shell=True, stderr=log, stdin=log, stdout=log) # parm7 rst7 files made
            # Run Min
            self.min_script = self.min_script.replace("aa","700")
            self.min_script = self.min_script.replace("bb","5")
            with open("./calc/"+self.name_min, 'w') as f:
                f.write(self.min_script)
            sp.call("sander -O -i calc/{} -o calc/{} -p calc/prmtop -c calc/inpcrd -r calc/amber.{}.rst".format(self.name_min,self.name_ene,self.name),shell=True, stderr=log, stdin=log, stdout=log)
        elif com_opts['sp']:
            logger.log(1, '  CALCULATE: {}'.format(self.filename))
            # Run leap
            sp.call("tleap -f {}".format(self.filename),shell=True, stderr=log, stdin=log, stdout=log) # parm7 rst7 files made
            # Run Min
            self.min_script = self.min_script.replace("aa","0")
            self.min_script = self.min_script.replace("bb","0")
            with open("./calc/"+self.name_min, 'w') as f:
                f.write(self.min_script)
            sp.call("sander -O -i calc/{} -o calc/{} -p calc/prmtop -c calc/inpcrd -r calc/amber.{}.rst".format(self.name_min,self.name_ene,self.name),shell=True, stderr=log, stdin=log, stdout=log)
        # check if energy calculation failed
        restart = 1
        while(restart==1):
            with open("./calc/"+self.name_ene,'r') as f:
                fline = f.readlines()
                for line in fline:
                    if "restarting should resolve the error" in line:
                        sp.call("sander -O -i calc/{} -o calc/{} -p calc/prmtop -c calc/amber.{}.rst -r calc/amber.{}.rst".format(self.name_min,self.name_ene,self.name,self.name),shell=True, stderr=log, stdin=log, stdout=log)
                        restart = 1
                    else:
                        restart = 0

        if com_opts['geo']:
            self.geometry(log)
        if com_opts['freq']:
            self.hessian(log)
            # if geo file is already present 
            # may not have geo file if hessian is only ran
            if os.path.isfile('./calc/'+self.name_geo):
                self.geo_extract()
        os.chdir(current_directory)


class TinkerHess(File):
    def __init__(self, path):
        super(TinkerHess, self).__init__(path)
        self._hessian = None
        self.natoms = None
    @property
    def hessian(self):
        if self._hessian is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            hessian = np.zeros([self.natoms * 3, self.natoms * 3], dtype=float)
            logger.log(5, '  -- Creatting {} Hessian Matrix.'.format(
                hessian.shape))
            with open(self.path, 'r') as f:
                lines = f.read()
            words = lines.split()
            diag = True
            row_num = 0
            col_num = 0
            line = -1
            index = 0
            for i, word in enumerate(words):
                match = re.compile('\d+[.]\d+').search(word)
                # First group of values are all of the diagonal elements. So
                # This will grab them first and put them in the correct index
                # of the Hessian.
                if diag and match:
                    hessian[row_num, col_num] = word
                    row_num += 1
                    col_num += 1
                # After the first group of values the line will read
                # 'Off-diagonal'. This signifies when the next elements are
                # H_i, j for section i.
                if word == 'Off-diagonal':
                    diag = False
                    line += 1
                    index = line + 1
                    row_num = 0
                    col_num = 0
                if not diag and match:
                    hessian[line, col_num + index] = word
                    hessian[row_num + index, line] = word
                    row_num += 1
                    col_num += 1
            # Convert hessian units to use kJ/mol instead of kcal/mol.
            self._hessian = hessian / co.HARTREE_TO_KCALMOL \
                * co.HARTREE_TO_KJMOL
            logger.log(5, '  -- Finished Creating {} Hessian matrix.'.format(
                hessian.shape))
            return self._hessian

class TinkerLog(File):
    def __init__(self, path):
        super(TinkerLog, self).__init__(path)
        self._structures = None
        self.name = None
    @property
    def structures(self):
        if self._structures == None:
            logger.log(10, 'READING: {}'.format(self.filename))
            self._structures = []
            with open(self.path, 'r') as f:
                sections = {'sp':1, 'minimization':2, 'hessian':2}
                count_previous = 0
                calc_section = 'sp'
                for line in f:
                    count_current = sections[calc_section]
                    if count_current != count_previous:
                        # Due to TINKER printing to standard error and standard
                        # out the redirection of this printout to the *q2mm.log
                        # does not print the bonds and angles in the correct
                        # order or in a consistent order. Therefore we need to
                        # sort the bonds, angles, and torsions in a standard
                        # way before extending the data in calculate.
                        bonds = []
                        angles = []
                        torsions = []
                        current_structure = Structure()
                        self._structures.append(current_structure)
                        count_previous += 1
                    section = None
                    if "SINGLE POINT" in line:
                        calc_section = 'minimization'
                        for bond in bonds:
                            # Not sure if I have to sort the atom list but
                            # I'm doing it anyway.
                            bond.atom_nums.sort()
                        # Sorts the bonds by the first atom and then by
                        # the second.
                        bonds.sort(key=lambda x: (x.atom_nums[0],
                                                  x.atom_nums[1]))
                        for angle in angles:
                            if angle.atom_nums[0] > angle.atom_nums[2]:
                                angle.atom_nums = [angle.atom_nums[2],
                                                   angle.atom_nums[1],
                                                   angle.atom_nums[0]]
                        for torsion in torsions:
                            if torsion.atom_nums[1] > torsion.atom_nums[2]:
                                torsion.atom_nums = [torsion.atom_nums[3],
                                                     torsion.atom_nums[2],
                                                     torsion.atom_nums[1],
                                                     torsion.atom_nums[0]]
                        angles.sort(key=lambda x: (x.atom_nums[1],
                                                   x.atom_nums[0],
                                                   x.atom_nums[2]))
                        torsions.sort(key=lambda x: (x.atom_nums[1],
                                                     x.atom_nums[2],
                                                     x.atom_nums[0],
                                                     x.atom_nums[3]))
                        current_structure.bonds.extend(bonds)
                        current_structure.angles.extend(angles)
                        current_structure.torsions.extend(torsions)
                    if 'END OF OPTIMIZED SINGLE POINT' in line:
                        calc_section = 'hessian'
                    if 'Bond' in line:
                        bond = self.read_line_for_bond(line)
                        if bond is not None:
                            bonds.append(bond)
                    if 'Angle' in line:
                        angle = self.read_line_for_angle(line)
                        if angle is not None:
                            angles.append(angle)
                    if 'Torsion' in line:
                        torsion = self.read_line_for_torsion(line)
                        if torsion is not None:
                            torsions.append(torsion)
                    if 'Total Potential Energy' in line:
                        energy = self.read_line_for_energy(line)
                        if energy is not None:
                            current_structure.props['energy']=energy
                    if "END OF CALCULATION" in line and \
                        calc_section != 'hessian':
                        last_ind = len(self._structures) - 1
                        self._structures.remove(self._structures[last_ind])
            logger.log(5, '  -- Imported {} structure(s)'.format(
                len(self._structures)))
        return self._structures

    def read_line_for_bond(self, line):
        # All bond data starts with the string "Bond" and then the rest of the
        # interaction information.
        match = re.compile('Bond\s+(\d+)-(\w+)\s+(\d+)-(\w+)\s+'
            '({0})\s+({0})\s+({0})'.format(co.RE_FLOAT)).search(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(3)]]
            value = float(match.group(6))
            return Bond(atom_nums=atom_nums, value=value)
        else:
            return None

    def read_line_for_angle(self, line):
        # All angle data starts with the string "Angle" and then the rest of the
        # interaction information.
        match = re.compile('Angle\s+(\d+)-(\w+)\s+(\d+)-(\w+)\s+(\d+)-(\w+)\s+'
            '({0})\s+({0})\s+({0})'.format(co.RE_FLOAT)).search(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(3),
                match.group(5)]]
            value = float(match.group(8))
            return Angle(atom_nums=atom_nums, value=value)
        else:
            return None

    def read_line_for_torsion(self, line):
        # All torsion data starts with the string "torsion" and then the rest of
        # the interaction information.
        match = re.compile('Torsion\s+(\d+)-(\w+)\s+(\d+)-(\w+)\s+'
            '(\d+)-(\w+)\s+(\d+)-(\w+)\s+({0})\s+({0})'.format(
                co.RE_FLOAT)).search(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(3),
                match.group(5), match.group(7)]]
            value = float(match.group(9))
            return Torsion(atom_nums=atom_nums, value=value)
        else:
            return None

    def read_line_for_energy(self, line):
        # The TPE is in units of kcal/mol, so we have to convert them to kJ/mol
        # for consistency purposes.
        match = re.compile('Total Potential Energy :\s+({0})'.format(
            co.RE_FLOAT)).search(line)
        if match:
            energy = float(match.group(1))
            energy *= co.HARTREE_TO_KJMOL / co.HARTREE_TO_KCALMOL
            return energy
        else:
            return None

class TinkerXYZ(File):
    def __init__(self, path):
        super(TinkerXYZ, self).__init__(path)
        self._index_output_log = None
        self._structures = None
        self.commands = None
        self.name = os.path.splitext(self.filename)[0]
        # Key file is needed to set the settings for the calculation, including
        # the parameters needed to perform the calculation.
        self.name_key = self.name + '.q2mm.key'
        # The log file is a file that contains the information redirected from
        # the TINKER calculations that are performed with Q2MM. This is not a
        # file setup by TINKER. TINKER will only print to the front end except
        # for select files such as a newly minimized structure. In this case
        # the minimized structure will be saved to *.q2mm.xyz.
        self.name_log = self.name + '.q2mm.log'
        self.name_xyz = self.name + '.q2mm.xyz'
        self.name_hes = self.name + '.q2mm.hes'
        self.name_1st_hess = self.name + '.hes'
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            struct = Structure()
            self._structures = [struct]
            with open(self.filename, 'r') as f:
                for line in f:
                    line = line.split()
                    if len(line) == 2:
                        struct.props['total atoms'] = int(line[0])
                        struct.props['title'] = line[1]
                        logger.log(5, '  -- Read {} atoms.'.format(
                            struct.props['total atoms']))
                    if len(line) > 2:
                        indx, ele, x, y, z, at, bonded_atom = line[0], \
                            line[1], line[2], line[3], line[4], \
                            line[5], line[6:]
                        struct.atoms.append(Atom(index=indx,
                            element=ele,
                            x=float(x),
                            y=float(y),
                            z=float(z),
                            atom_type=at,
                            atom_type_name=at,
                            bonded_atom_indices=bonded_atom))
            return self._structures
    def get_com_opts(self):
        com_opts = {'freq': False,
                    'opt': False,
                    'sp': False,
                    'tors': False}
        if any(x in ['tb', 'ta', 'tt', 'te', 'tea'] for x in self.commands):
            com_opts['sp'] = True
        if any(x in ['tbo','tao','tto','teo','teao'] for x in self.commands):
            com_opts['opt'] = True
            com_opts['sp'] = True
        if any(x in ['th', 'tjeig', 'tgeig'] for x in self.commands):
            com_opts['freq'] = True
            com_opts['opt'] = True
            com_opts['sp'] = True
        if any(x in ['tt', 'tto'] for x in self.commands):
            com_opts['tors'] = True
        return com_opts
    def run(self,check_tokens=False):
        logger.log(5, 'RUNNING: {}'.format(self.filename))
        self._index_output_log = []
        com_opts = self.get_com_opts()
        current_directory = os.getcwd()
        os.chdir(self.directory)
        if os.path.isfile(self.name_key):
            os.remove(self.name_key)
        if os.path.isfile(self.name_log):
            os.remove(self.name_log)
        if os.path.isfile(self.name_xyz):
            os.remove(self.name_xyz)
        if os.path.isfile(self.name_hes):
            os.remove(self.name_heis)
        with open(self.name_key, 'w') as f:
            f.write('parameters mm3.prm\
                     \nnoversion')
        if com_opts['sp']:
            logger.log(1, '  ANALYZE: {}'.format(self.filename))
            with open(self.name_log, 'w') as f:
                sp.call(
                    'analyze {}.xyz -k {} D'.format(self.name,
                    self.name_key), shell=True, stderr=f, stdin=f, stdout=f)
                # Not sure if these print outs are important, but I add in these
                # to define what section of the *q2mm.log file you are in. This
                # is especially helpful when needed to grab all of the bond,
                # angle, and torsional data since these have to be ordered
                # everytime they are grabbed.
                f.write("\n=======================\
                         \n= END OF SINGLE POINT =\
                         \n=======================\n")
        if com_opts['opt']:
            logger.log(1, '  MINIMIZE & ANALYZE: {}'.format(self.filename))
            with open(self.name_log, 'a') as f:
                # The float value is the convergence criteria.
                sp.call(
                    'minimize {}.xyz -k {} 0.01 {}'.format(self.name,
                    self.name_key, self.name_xyz), shell=True, stderr=f,
                    stdin=f, stdout=f)
                sp.call(
                    'analyze {} -k {} D'.format(self.name_xyz,
                    self.name_key), shell=True, stderr=f, stdin=f, stdout=f)
                f.write("\n=================================\
                         \n= END OF OPTIMIZED SINGLE POINT =\
                         \n=================================\n")
        if com_opts['freq']:
            logger.log(1, '  TESTHESS: {}'.format(self.filename))
            with open(self.name_log, 'a') as f:
                # Tinker will not take a file output argument if the there isn't
                # currently a file.hes. For example, file.xyz will write to
                # file.hes if file.hes doesn't already exist otherwise TINKER
                # will ask the user for a new filename.
                if os.path.isfile(self.name + '.hes'):
                    sp.call(
                        'testhess {} -k {} y n {}'.format(self.name,
                        self.name_key,self.name_hes), shell=True,
                        stderr=f, stdin=f, stdout=f)
                else:
                    sp.call(
                        'testhess {} -k {} y n'.format(self.name,
                        self.name_key), shell=True,
                        stderr=f, stdin=f, stdout=f)
                    os.rename(self.name + '.hes', self.name_hes)
                f.write("\n==================\
                         \n= END OF HESSIAN =\
                         \n==================\n")
        with open(self.name_log, 'a') as f:
            f.write("\n======================\
                     \n= END OF CALCULATION =\
                     \n======================\n")
        os.chdir(current_directory
)

class TinkerXYZ_FOR_GAUS(File):
    def __init__(self, path):
        super(TinkerXYZ_FOR_GAUS, self).__init__(path)
        self._index_output_log = None
        self._structures = None
        self.commands = None
        self.name = os.path.splitext(self.filename)[0]
        self.ref_xyz = self.name + ".xyz"
        self.name_key = "gaus." +self.name + '.q2mm.key'
        self.name_log = "gaus." +self.name + '.q2mm.log'
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            struct = Structure()
            self._structures = [struct]
            with open(self.filename, 'r') as f:
                for line in f:
                    line = line.split()
                    if len(line) == 2:
                        struct.props['total atoms'] = int(line[0])
                        struct.props['title'] = line[1]
                        logger.log(5, '  -- Read {} atoms. (GAUSSIAN2TINKER)'.format(
                            struct.props['total atoms']))
                    if len(line) > 2:
                        indx, ele, x, y, z, at, bonded_atom = line[0], \
                            line[1], line[2], line[3], line[4], \
                            line[5], line[6:]
                        struct.atoms.append(Atom(index=indx,
                            element=ele,
                            x=float(x),
                            y=float(y),
                            z=float(z),
                            atom_type=at,
                            atom_type_name=at,
                            bonded_atom_indices=bonded_atom))
            return self._structures
    def get_com_opts(self):
        com_opts = {'freq': False,
                    'opt': False,
                    'sp': False,
                    'tors': False}
        if any(x in ['gtb', 'gta', 'gtt'] for x in self.commands):
            com_opts['sp'] = True
        return com_opts
    def tinkerxyz_from_gaussian(self):
        ref_xyz = list()
        # Gaussian log to xyz coordinate
        with open(self.filename) as f:
            lines = f.readlines()
            xyz = 0
            comma = 0
            for i in range(len(lines)):
                line = lines[i]
                if xyz and comma and "," not in line:
                    #continue
                    break
                elif xyz and "," in line:
                    comma = 1
                    line = line.replace('\n','')
                    l = line.split(',')
                    ref_xyz.append(l[2:])
                elif "Charge =" in line:
                    xyz = 1
                
            
        self.filename = "gaus." + self.name + '.xyz' #replace xyz
        with open(self.ref_xyz,'r') as ref:
            with open(self.filename,'w') as f:
                lines = ref.readlines()
                for i in range(len(lines)):
                    line = lines[i]
                    l = line.split()
                    if i == 0:
                        f.write(line)
                    else:
                        if len(l) > 5:
                            l[2:5] = ref_xyz[i-1]
                            f.write("\t".join(l)+'\n')
        return
    def run(self,check_tokens=False):
        logger.log(5, 'RUNNING: {}'.format(self.filename))
        self.tinkerxyz_from_gaussian()
        self._index_output_log = []
        com_opts = self.get_com_opts()
        current_directory = os.getcwd()
        os.chdir(self.directory)
        if os.path.isfile(self.name_key):
            os.remove(self.name_key)
        if os.path.isfile(self.name_log):
            os.remove(self.name_log)
        with open(self.name_key, 'w') as f:
            f.write('parameters mm3.prm\
                     \nnoversion')
        if com_opts['sp']:
            logger.log(1, '  ANALYZE: {}'.format(self.filename))
            print("Gaussian Ref xyz file Running Tinker for Geo. Info.")
            with open(self.name_log, 'w') as f:
                sp.call(
                    'analyze gaus.{}.xyz -k {} D'.format(self.name,
                    self.name_key), shell=True, stderr=f, stdin=f, stdout=f)
                f.write("\n=======================\
                         \n= END OF SINGLE POINT =\
                         \n=======================\n")
        with open(self.name_log, 'a') as f:
            f.write("\n======================\
                     \n= END OF CALCULATION =\
                     \n======================\n")
        os.chdir(current_directory)

class GaussCom(File):
    """
    Used to write Gaussian command files to run quick calculations. The ony use
    thus far is for checking the RMS of the electrostatic potential using 
    partial charges obtained from a force field.
    """
    def __init__(self, path):
        super(GaussCom, self).__init__(path)
        self.name = os.path.splitext(self.filename)[0]
        self.commands = None
        self.old_chk = None
        self.name_chk = self.name + '.chk'
        self.name_log = self.name + '.log'
        self.atom_and_coords = []
        self.charge = None
        self.charge_list = None
        self.multiplicity = 1
        self.memory = 24
        self.procs = 24
        self.method = 'M06'

    def read_newzmat(self, old_check):
        self.old_chk = old_check
        if os.path.isfile('TEMP.com'):
            os.remove('TEMP.com')
        sp.call('newzmat -ichk -ocart {} {}'.format(self.old_chk, 'TEMP'), 
            shell=True)
        with open('TEMP.com', 'r') as f:
            for line in f:
                if re.search('[A-Z][a-z]?\s+{0}\s+{0}\s+{0}'.format(
                    co.RE_FLOAT),line):
                    ele, x, y, z = line.split()
                    self.atom_and_coords.append([ele,float(x),float(y),float(z)])
                if re.search('[+-]?\d,\d',line):
                    charge, multiplicity = line.split(',')
                    self.charge = charge
                    self.multiplicity = multiplicity    
                if re.search('^[#]',line):
                    split = line.split()
                    for keyword in split:
                        for method in co.gaussian_methods:
                            if method.upper() in keyword.upper():
                                self.method = method.upper()
        #os.remove('TEMP.com')

    def write_com(self):
        """
        Gaussian has a utility, newzmat, which could write the com file, but
        in my experience it is rather slow especially for larger chkpoint files.
        It is easier just writting the file ourselves. -Tony
        """
        elements_in_structure = []
        with open(self.filename, 'w') as f:
            if self.old_chk:
                f.write('%oldchk={}\n'.format(self.old_chk))
            f.write('%chk={}\n'.format(self.name_chk))
            f.write('%mem={}GB\n'.format(self.memory))
            f.write('%nprocshared={}\n'.format(self.procs))
            f.write('#N Guess=TCheck SCRF=Check GenChk pop=(chelpg,readradii) '
                    'IOp(6/20=30133) chkbasis {}\n\n'.format(self.method))
            f.write('SP with fixed charges\n\n')
            f.write('{} {}'.format(self.charge,self.multiplicity))
            if self.charge_list:
                try:
                    if len(self.charge_list) != len(self.atom_and_coords):
                        raise ValueError("Length of the charge list and atom "
                                         "list are not the same: {}".format(
                                         self.filename))
                except:
                    print("Just print something")
                for atom, charge in zip(self.atom_and_coords,self.charge_list):
                    if atom[0] not in elements_in_structure:
                        elements_in_structure.append(atom[0])
                    f.write('   '.join((' {}--{:8.6f}'.format(atom[0],charge),
                                        '{:14.10f}'.format(atom[1]),
                                        '{:14.10f}'.format(atom[2]),
                                        '{:14.10f}'.format(atom[3]),
                                        '\n'
                                        )))
            f.write('\n')
            for atom in elements_in_structure:
                f.write('{0} {1}\n'.format(atom,co.CHELPG_RADII[atom]))
            f.write('\n')

    def run_gaussian(self):
        """
        This runs the gaussian job. I have this labeled differently than our 
        typical "run" functions because I don't want to use this function until
        after we have calculated and collected partial charge data.
        """
        logger.log(5, 'RUNNING: {}'.format(self.filename))
        self._index_output_log = []
        current_directory = os.getcwd()
        os.chdir(self.directory)
        if os.path.isfile(self.name_log):
            os.remove(self.name_log)
        if os.path.isfile(self.name_chk):
            os.remove(self.name_chk)
        sp.call('g09 {}'.format(self.filename), shell=True)
        os.chdir(current_directory)

class GaussFormChk(File):
    """
    Used to retrieve data from Gaussian formatted checkpoint files.
    """
    def __init__(self, path):
        super(GaussFormChk, self).__init__(path)
        self.atoms = []
        # Not sure these should really be called the eigenvalues.
        self.evals = None
        self.low_tri = None
        self._hess = None
    @property
    def hess(self):
        if self._hess is None:
            self.read_self()
        return self._hess
    def read_self(self):
        logger.log(5, 'READING: {}'.format(self.filename))
        stuff = re.search(
            'Atomic numbers\s+I\s+N=\s+(?P<num_atoms>\d+)'
            '\n\s+(?P<anums>.*?)'
            'Nuclear charges.*?Current cartesian coordinates.*?\n(?P<coords>.*?)'
            'Force Field'
            '.*?Real atomic weights.*?\n(?P<masses>.*?)'
            'Atom fragment info.*?Cartesian Gradient.*?\n(?P<evals>.*?)'
            'Cartesian Force Constants.*?\n(?P<hess>.*?)'
            'Dipole Moment',
            open(self.path, 'r').read(), flags=re.DOTALL)
        anums = [int(x) for x in stuff.group('anums').split()]
        masses = [float(x) for x in stuff.group('masses').split()]
        coords = [float(x) for x in stuff.group('coords').split()]
        coords = [coords[i:i+3] for i in range(0, len(coords), 3)]
        for anum, mass, coord in zip(anums, masses, coords):
            self.atoms.append(
                Atom(
                    atomic_num = anum,
                    coords = coord,
                    exact_mass = mass)
                )
        logger.log(5, '  -- Read {} atoms.'.format(len(self.atoms)))
        self.evals = np.array(
            [float(x) for x in stuff.group('evals').split()], dtype=float)
        logger.log(5, '  -- Read {} eigenvectors.'.format(len(self.evals)))
        self.low_tri = np.array(
            [float(x) for x in stuff.group('hess').split()], dtype=float)
        one_dim = len(anums) * 3
        self._hess = np.empty([one_dim, one_dim], dtype=float)
        self._hess[np.tril_indices_from(self._hess)] = self.low_tri
        self._hess += np.tril(self._hess, -1).T
        # Convert to MacroModel units.
        self._hess *= co.HESSIAN_CONVERSION
        logger.log(5, '  -- Read {} Hessian.'.format(self._hess.shape))

class GaussLog(File):
    """
    Used to retrieve data from Gaussian log files.

    If you are extracting frequencies/Hessian data from this file, use
    the keyword NoSymmetry when running the Gaussian calculation.
    """
    def __init__(self, path):
        super(GaussLog, self).__init__(path)
        self._evals = None
        self._evecs = None
        self._structures = None
        self._esp_rms = None
    @property
    def evecs(self):
        if self._evecs is None:
            self.read_out()
        return self._evecs
    @property
    def evals(self):
        if self._evals is None:
            self.read_out()
        return self._evals
    @property
    def structures(self):
        if self._structures is None:
            #self.read_out()
            self.read_archive()
        return self._structures
    @property
    def esp_rms(self):
        if self._esp_rms is None:
            self._esp_rms = -1
            self.read_out()
        return self._esp_rms
    def read_out(self):
        """
        Read force constant and eigenvector data from a frequency
        calculation.
        """
        logger.log(5, 'READING: {}'.format(self.filename))
        self._evals = []
        self._evecs = []
        self._structures = []
        force_constants = []
        evecs = []
        with open(self.path, 'r') as f:
            # The keyword "harmonic" shows up before the section we're
            # interested in. It can show up multiple times depending on the
            # options in the Gaussian .com file.
            past_first_harm = False
            # High precision mode, turned on by including "freq=hpmodes" in the
            # Gaussian .com file.
            hpmodes = False
            file_iterator = iter(f)
            # This while loop breaks when the end of the file is reached, or
            # if the high quality modes have been read already.
            while True:
                try:
                    line = next(file_iterator)
                except:
                    # End of file.
                    break
                if 'Charges from ESP fit' in line:
                    pattern = re.compile('RMS=\s+({0})'.format(co.RE_FLOAT))
                    match = pattern.search(line)
                    self._esp_rms = float(match.group(1))
                # Gathering some geometric information.
                elif 'orientation:' in line:
                    self._structures.append(Structure())
                    next(file_iterator)
                    next(file_iterator)
                    next(file_iterator)
                    next(file_iterator)
                    line = next(file_iterator)
                    while not '---' in line:
                        cols = line.split()
                        self._structures[-1].atoms.append(
                            Atom(atomic_num=int(cols[1]),
                                 x=float(cols[3]),
                                 y=float(cols[4]),
                                 z=float(cols[5])))
                        line = next(file_iterator)
                    logger.log(5, '  -- Found {} atoms.'.format(
                            len(self._structures[-1].atoms)))
                elif 'Harmonic' in line:
                    # The high quality eigenvectors come before the low quality
                    # ones. If you see "Harmonic" again, it means you're at the
                    # low quality ones now, so break.
                    if past_first_harm:
                        break
                    else:
                        past_first_harm = True
                elif 'Frequencies' in line:
                    # We're going to keep reusing these.
                    # We accumulate sets of eigevectors and eigenvalues, add
                    # them to self._evecs and self._evals, and then reuse this
                    # for the next set.
                    del(force_constants[:])
                    del(evecs[:])
                    # Values inside line look like:
                    #     "Frequencies --- xxxx.xxxx xxxx.xxxx"
                    # That's why we remove the 1st two columns. This is
                    # consistent with and without "hpmodes".
                    # For "hpmodes" option, there are 5 of these frequencies.
                    # Without "hpmodes", there are 3.
                    # Thus the eigenvectors and eigenvalues will come in sets of
                    # either 5 or 3.
                    cols = line.split()
                    for frequency in map(float, cols[2:]):
                        # Has 1. or -1. depending on the sign of the frequency.
                        if frequency < 0.:
                            force_constants.append(-1.)
                        else:
                            force_constants.append(1.)
                        # For now this is empty, but we will add to it soon.
                        evecs.append([])

                    # Moving on to the reduced masses.
                    line = next(file_iterator)
                    cols = line.split()
                    # Again, trim the "Reduced masses ---".
                    # It's "Red. masses --" for without "hpmodes".
                    for i, mass in enumerate(map(float, cols[3:])):
                        # +/- 1 / reduced mass
                        force_constants[i] = force_constants[i] / mass

                    # Now we are on the line with the force constants.
                    line = next(file_iterator)
                    cols = line.split()
                    # Trim "Force constants ---". It's "Frc consts --" without
                    # "hpmodes".
                    for i, force_constant in enumerate(map(float, cols[3:])):
                        # co.AU_TO_MDYNA = 15.569141
                        force_constants[i] *= force_constant / co.AU_TO_MDYNA

                    # Force constants were calculated above as follows:
                    #    a = +/- 1 depending on the sign of the frequency
                    #    b = a / reduced mass (obtained from the Gaussian log)
                    #    c = b * force constant / conversion factor (force
                    #         (constant obtained from Gaussian log) (conversion
                    #         factor is inside constants module)

                    # Skip the IR intensities.
                    next(file_iterator)
                    # This is different depending on whether you use "hpmodes".
                    line = next(file_iterator)
                    # "Coord" seems to only appear when the "hpmodes" is used.
                    if 'Coord' in line:
                        hpmodes = True
                    # This is different depending on whether you use
                    # "freq=projected".
                    line = next(file_iterator)
                    # The "projected" keyword seems to add "IRC Coupling".
                    if 'IRC Coupling' in line:
                        line = next(file_iterator)
                    # We're on to the eigenvectors.
                    # Until the end of this section containing the eigenvectors,
                    # the number of columns remains constant. When that changes,
                    # we know we're to the next set of frequencies, force
                    # constants and eigenvectors.
                    cols = line.split()
                    cols_len = len(cols)

                    while len(cols) == cols_len:
                        # This will come after all the eigenvectors have been
                        # read. We can break out then.
                        if 'Harmonic' in line:
                            break
                        # If "hpmodes" is used, you have an extra column here
                        # that is simply an index.
                        if hpmodes:
                            cols = cols[1:]
                        # cols corresponds to line(s) (maybe only 1st line)
                        # under section "Coord Atom Element:" (at least for
                        # "hpmodes").

                        # Just the square root of the mass from co.MASSES.
                        # co.MASSES currently has the average mass.
                        # Gaussian may use the mass of the most abundant
                        # isotope. This may be a problem.
                        mass_sqrt = np.sqrt(list(co.MASSES.items())[int(cols[1]) - 1][1])

                        cols = cols[2:]
                        # This corresponds to the same line still, but without
                        # the atom elements.

                        # This loop expands the LoL, evecs, as so.
                        # Iteration 1:
                        # [[x], [x], [x], [x], [x]]
                        # Iteration 2:
                        # [[x, x], [x, x], [x, x], [x, x], [x, x]]
                        # ... etc. until the length of the sublist is equal to
                        # the number of atoms. Remember, for low precision
                        # eigenvectors it only adds in sets of 3, not 5.

                        # Elements of evecs are simply the data under
                        # "Coord Atom Element" multiplied by the square root
                        # of the weight.
                        for i in range(len(evecs)):
                            if hpmodes:
                                # evecs is a LoL. Length of sublist is
                                # equal to # of columns in section "Coord Atom
                                # Element" minus 3, for the 1st 3 columns
                                # (index, atom index, atomic number).
                                evecs[i].append(float(cols[i]) * mass_sqrt)
                            else:
                                # This is fow low precision eigenvectors. It's a
                                # funny way to go in sets of 3. Take a look at
                                # your low precision Gaussian log and it will
                                # make more sense.
                                for useless in range(3):
                                    x = float(cols.pop(0))
                                    evecs[i].append(x * mass_sqrt)
                        line = next(file_iterator)
                        cols = line.split()

                    # Here the overall number of eigenvalues and eigenvectors is
                    # increased by 5 (high precision) or 3 (low precision). The
                    # total number goes to 3N - 6 for non-linear and 3N - 5 for
                    # linear. Same goes for self._evecs.
                    for i in range(len(evecs)):
                        self._evals.append(force_constants[i])
                        self._evecs.append(evecs[i])
                    # We know we're done if this is in the line.
                    if 'Harmonic' in line:
                        break
        if self._evals and self._evecs:
            for evec in self._evecs:
                # Each evec is a single eigenvector.
                # Add up the sum of squares over an eigenvector.
                sum_of_squares = 0.
                # Appropriately named, element is an element of that single
                # eigenvector.
                for element in evec:
                    sum_of_squares += element * element
                # Now x is the inverse of the square root of the sum of squares
                # for an individual eigenvector.
                element = 1 / np.sqrt(sum_of_squares)
                for i in range(len(evec)):
                    evec[i] *= element
            self._evals = np.array(self._evals)
            self._evecs = np.array(self._evecs)
            logger.log(1, '>>> self._evals: {}'.format(self._evals))
            logger.log(1, '>>> self._evecs: {}'.format(self._evecs))
            logger.log(5, '  -- {} structures found.'.format(
                len(self.structures)))
    # May want to move some attributes assigned to the structure class onto
    # this filetype class.
    def read_archive(self):
        """
        Only reads last archive found in the Gaussian .log file.
        """
        logger.log(5, 'READING: {}'.format(self.filename))
        struct = Structure()
        self._structures = [struct]
        # Matches everything in between the start and end.
        # (?s)  - Flag for re.compile which says that . matches all.
        # \\\\  - One single \
        # Start - " 1\1\".
        # End   - Some number of \ followed by @. Not sure how many \ there
        #         are, so this matches as many as possible. Also, this could
        #         get separated by a line break (which would also include
        #         adding in a space since that's how Gaussian starts new lines
        #         in the archive).
        # We pull out the last one [-1] in case there are multiple archives
        # in a file.
#        print(self.path)
#        print(open(self.path,'r').read())
#        print(re.findall('(?s)(\s1\\\\1\\\\.*?[\\\\\n\s]+@)',open(self.path,'r').read()))
        try:
            arch = re.findall(
                '(?s)(\s1\\\\1\\\\.*?[\\\\\n\s]+@)',
                open(self.path, 'r').read())[-1]
            logger.log(5, '  -- Located last archive.')
        except IndexError:
            logger.warning("  -- Couldn't locate archive.")
            raise
        # Make it into one string.
        arch = arch.replace('\n ', '')
        # Separate it by Gaussian's section divider.
        arch = arch.split('\\\\')
        # Helps us iterate over sections of the archive.
        section_counter = 0
        # SECTION 0
        # General job information.
        arch_general = arch[section_counter]
        section_counter += 1
        stuff = re.search(
            '\s1\\\\1\\\\.*?\\\\.*?\\\\.*?\\\\.*?\\\\.*?\\\\(?P<user>.*?)'
            '\\\\(?P<date>.*?)'
            '\\\\.*?',
            arch_general)
        struct.props['user'] = stuff.group('user')
        struct.props['date'] = stuff.group('date')
        # SECTION 1
        # The commands you wrote.
        arch_commands = arch[section_counter]
        section_counter += 1
        # SECTION 2
        # The comment line.
        arch_comment = arch[section_counter]
        section_counter += 1
        # SECTION 3
        # Actually has charge, multiplicity and coords.
        arch_coords = arch[section_counter]
        section_counter +=1
        stuff = re.search(
            '(?P<charge>.*?)'
            ',(?P<multiplicity>.*?)'
            '\\\\(?P<atoms>.*)',
            arch_coords)
        struct.props['charge'] = stuff.group('charge')
        struct.props['multiplicity'] = stuff.group('multiplicity')
        # We want to do more fancy stuff with the atoms than simply add to
        # the properties dictionary.
        atoms = stuff.group('atoms')
        atoms = atoms.split('\\')
        # Z-matrix coordinates adds another section. We need to be aware of
        # this.
        probably_z_matrix = False
        for atom in atoms:
            stuff = atom.split(',')
            # An atom typically looks like this:
            #    C,0.1135,0.13135,0.63463
            if len(stuff) == 4:
                ele, x, y, z = stuff
            # But sometimes they look like this (notice the extra zero):
            #    C,0,0.1135,0.13135,0.63463
            # I'm not sure what that extra zero is for. Anyway, ignore
            # that extra whatever if it's there.
            elif len(stuff) == 5:
                ele, x, y, z = stuff[0], stuff[2], stuff[3], stuff[4]
            # And this would be really bad. Haven't seen anything else like
            # this yet.
            # 160613 - So, not sure when I wrote that comment, but something
            # like this definitely happens when using scans and z-matrices.
            # I'm going to ignore grabbing any atoms in this case.
            else:
                logger.warning(
                    'Not sure how to read coordinates from Gaussian acrhive!')
                probably_z_matrix = True
                section_counter += 1
                # Let's have it stop looping over atoms, but not fail anymore.
                break
                # raise Exception(
                #     'Not sure how to read coordinates from Gaussian archive!')
            struct.atoms.append(
                Atom(element=ele, x=float(x), y=float(y), z=float(z)))
        logger.log(20, '  -- Read {} atoms.'.format(len(struct.atoms)))
        # SECTION 4
        # All sorts of information here. This area looks like:
        #     prop1=value1\prop2=value2\prop3=value3
        arch_info = arch[section_counter]
        section_counter += 1
        arch_info = arch_info.split('\\')
        for thing in arch_info:
            prop_name, prop_value = thing.split('=')
            struct.props[prop_name] = prop_value
        # SECTION 5
        # The Hessian. Only exists if you did a frequency calculation.
        # Appears in lower triangular form.
        if not arch[section_counter] == '@':
            hess_tri = arch[section_counter]
            hess_tri = hess_tri.split(',')
            logger.log(
                5,
                '  -- Read {} Hessian elements in lower triangular '
                'form.'.format(len(hess_tri)))
            hess = np.zeros([len(atoms) * 3, len(atoms) * 3], dtype=float)
            logger.log(
                5, '  -- Created {} Hessian matrix.'.format(hess.shape))
            # Code for if it was in upper triangle (it's not).
            # hess[np.triu_indices_from(hess)] = hess_tri
            # hess += np.triu(hess, -1).T
            # Lower triangle code.
            hess[np.tril_indices_from(hess)] = hess_tri
            hess += np.tril(hess, -1).T
            hess *= co.HESSIAN_CONVERSION
            struct.hess = hess
            # SECTION 6
            # Not sure what this is.

        # stuff = re.search(
        #     '\s1\\\\1\\\\.*?\\\\.*?\\\\.*?\\\\.*?\\\\.*?\\\\(?P<user>.*?)'
        #     '\\\\(?P<date>.*?)'
        #     '\\\\.*?\\\\\\\\(?P<com>.*?)'
        #     '\\\\\\\\(?P<filename>.*?)'
        #     '\\\\\\\\(?P<charge>.*?)'
        #     ',(?P<multiplicity>.*?)'
        #     '\\\\(?P<atoms>.*?)'
        #     # This marks the end of what always shows up.
        #     '\\\\\\\\'
        #     # This stuff sometimes shows up.
        #     # And it breaks if it doesn't show up.
        #     '.*?HF=(?P<hf>.*?)'
        #     '\\\\.*?ZeroPoint=(?P<zp>.*?)'
        #     '\\\\.*?Thermal=(?P<thermal>.*?)'
        #     '\\\\.*?\\\\NImag=\d+\\\\\\\\(?P<hess>.*?)'
        #     '\\\\\\\\(?P<evals>.*?)'
        #     '\\\\\\\\\\\\',
        #     arch)
        # logger.log(5, '  -- Read archive.')
        # atoms = stuff.group('atoms')
        # atoms = atoms.split('\\')
        # for atom in atoms:
        #     ele, x, y, z = atom.split(',')
        #     struct.atoms.append(
        #         Atom(element=ele, x=float(x), y=float(y), z=float(z)))
        # logger.log(5, '  -- Read {} atoms.'.format(len(atoms)))
        # self._structures = [struct]
        # hess_tri = stuff.group('hess')
        # hess_tri = hess_tri.split(',')
        # logger.log(
        #     5,
        #     '  -- Read {} Hessian elements in lower triangular '
        #     'form.'.format(len(hess_tri)))
        # hess = np.zeros([len(atoms) * 3, len(atoms) * 3], dtype=float)
        # logger.log(
        #     5, '  -- Created {} Hessian matrix.'.format(hess.shape))
        # # Code for if it was in upper triangle, but it's not.
        # # hess[np.triu_indices_from(hess)] = hess_tri
        # # hess += np.triu(hess, -1).T
        # # Lower triangle code.
        # hess[np.tril_indices_from(hess)] = hess_tri
        # hess += np.tril(hess, -1).T
        # hess *= co.HESSIAN_CONVERSION
        # struct.hess = hess
        # # Code to extract energies.
        # # Still not sure exactly what energies we want to use.
        # struct.props['hf'] = float(stuff.group('hf'))
        # struct.props['zp'] = float(stuff.group('zp'))
        # struct.props['thermal'] = float(stuff.group('thermal'))
    def get_most_converged(self, structures=None):
        """
        Used with geometry optimizations that don't succeed. Sometimes
        intermediate geometries obtain better convergence than the
        final geometry. This function returns the class Structure for
        the most converged geometry, which can then be used to output
        the coordinates for the next optimization.
        """
        if structures is None:
            structures = self.structures
        structures_compared = 0
        best_structure = None
        best_yes_or_no = None
        fields = ['RMS Force', 'RMS Displacement', 'Maximum Force',
                  'Maximum Displacement']
        for i, structure in reversed(list(enumerate(structures))):
            yes_or_no = [value[2] for key, value in structure.props.items()
                         if key in fields]
            if not structure.atoms:
                logger.warning('  -- No atoms found in structure {}. '
                               'Skipping.'.format(i+1))
                continue
            if len(yes_or_no) == 4:
                structures_compared += 1
                if best_structure is None:
                    logger.log(10, '  -- Most converged structure: {}'.format(
                            i+1))
                    best_structure = structure
                    best_yes_or_no = yes_or_no
                elif yes_or_no.count('YES') > best_yes_or_no.count('YES'):
                    best_structure = structure
                    best_yes_or_no = yes_or_no
                elif yes_or_no.count('YES') == best_yes_or_no.count('YES'):
                    number_better = 0
                    for field in fields:
                        if structure.props[field][0] < \
                                best_structure.props[field][0]:
                            number_better += 1
                    if number_better > 2:
                        best_structure = structure
                        best_yes_or_no = yes_or_no
            elif len(yes_or_no) != 0:
                logger.warning(
                    '  -- Partial convergence criterion in structure: {}'.format(
                        self.path))
        logger.log(10, '  -- Compared {} out of {} structures.'.format(
                structures_compared, len(self.structures)))
        return best_structure
    def read_any_coords(self, coords_type='both'):
        logger.log(10, 'READING: {}'.format(self.filename))
        structures = []
        with open(self.path, 'r') as f:
            section_coords_input = False
            section_coords_standard = False
            section_convergence = False
            section_optimization = False
            for i, line in enumerate(f):
                    # Look for input coordinates.
                    if coords_type == 'input' or coords_type == 'both':
                        # Marks end of input coords for a given structure.
                        if section_coords_input and 'Distance matrix' in line:
                            section_coords_input = False
                            logger.log(5, '[L{}] End of input coordinates '
                                       '({} atoms).'.format(
                                    i+1, count_atom))
                        # Add atoms and coordinates to structure.
                        if section_coords_input:
                            match = re.match(
                                '\s+(\d+)\s+(\d+)\s+\d+\s+({0})\s+({0})\s+'
                                '({0})'.format(co.RE_FLOAT), line)
                            if match:
                                count_atom += 1
                                try:
                                    current_atom = current_structure.atoms[
                                        int(match.group(1))-1]
                                except IndexError:
                                    current_structure.atoms.append(Atom())
                                    current_atom = current_structure.atoms[-1]
                                if current_atom.atomic_num:
                                    assert current_atom.atomic_num == int(
                                        match.group(2)), \
                                        ("[L{}] Atomic numbers don't match "
                                         "(current != existing) "
                                         "({} != {}).".format(
                                                i+1, int(match.group(2)),
                                                current_atom.atomic_num))
                                else:
                                    current_atom.atomic_num = int(
                                        match.group(2))
                                current_atom.coords_type = 'input'
                                current_atom.x = float(match.group(3))
                                current_atom.y = float(match.group(4))
                                current_atom.z = float(match.group(5))
                        # Start of input coords for a given structure.
                        if not section_coords_input and \
                                'Input orientation:' in line:
                            current_structure = Structure()
                            structures.append(current_structure)
                            section_coords_input = True
                            count_atom = 0
                            logger.log(5, '[L{}] Start input coordinates '
                                       'section.'.format(i+1))
                    # Look for standard coordinates.
                    if coords_type == 'standard' or coords_type == 'both':
                        # End of coordinates for a given structure.
                        if section_coords_standard and \
                                ('Rotational constants' in line or
                                 'Leave Link' in line):
                            section_coords_standard = False
                            logger.log(5, '[L{}] End standard coordinates '
                                       'section ({} atoms).'.format(
                                    i+1, count_atom))
                        # Grab coordinates for each atom.
                        # Add atoms to the structure.
                        if section_coords_standard:
                            match = re.match('\s+(\d+)\s+(\d+)\s+\d+\s+({0})\s+'
                                             '({0})\s+({0})'.format(
                                    co.RE_FLOAT), line)
                            if match:
                                count_atom += 1
                                try:
                                    current_atom = current_structure.atoms[
                                        int(match.group(1))-1]
                                except IndexError:
                                    current_structure.atoms.append(Atom())
                                    current_atom = current_structure.atoms[-1]
                                if current_atom.atomic_num:
                                    assert current_atom.atomic_num == int(
                                        match.group(2)), \
                                        ("[L{}] Atomic numbers don't match "
                                         "(current != existing) "
                                         "({} != {}).".format(
                                                i+1, int(match.group(2)),
                                                current_atom.atomic_num))
                                else:
                                    current_atom.atomic_num = int(
                                        match.group(2))
                                current_atom.coords_type = 'standard'
                                current_atom.x = float(match.group(3))
                                current_atom.y = float(match.group(4))
                                current_atom.z = float(match.group(5))
                        # Start of standard coordinates.
                        if not section_coords_standard and \
                                'Standard orientation' in line:
                            current_structure = Structure()
                            structures.append(current_structure)
                            section_coords_standard = True
                            count_atom = 0
                            logger.log(5, '[L{}] Start standard coordinates '
                                       'section.'.format(i+1))
        return structures
    def read_optimization(self, coords_type='both'):
        """
        Finds structures from a Gaussian geometry optimization that
        are listed throughout the log file. Also finds data about
        their convergence.

        coords_type = "input" or "standard" or "both"
                      Using both may cause coordinates in one format
                      to be overwritten by whatever comes later in the
                      log file.
        """
        logger.log(10, 'READING: {}'.format(self.filename))
        structures = []
        with open(self.path, 'r') as f:
            section_coords_input = False
            section_coords_standard = False
            section_convergence = False
            section_optimization = False
            for i, line in enumerate(f):
                # Look for start of optimization section of log file and
                # set a flag that it has indeed started.
                if section_optimization and 'Optimization stopped.' in line:
                    section_optimization = False
                    logger.log(5, '[L{}] End optimization section.'.format(i+1))
                if not section_optimization and \
                        'Search for a local minimum.' in line:
                    section_optimization = True
                    logger.log(5, '[L{}] Start optimization section.'.format(
                            i+1))
                if section_optimization:
                    # Start of a structure.
                    if 'Step number' in line:
                        structures.append(Structure())
                        current_structure = structures[-1]
                        logger.log(5, '[L{}] Added structure '
                                   '(currently {}).'.format(
                                i+1, len(structures)))
                    # Look for convergence information related to a single
                    # structure.
                    if section_convergence and 'GradGradGrad' in line:
                        section_convergence = False
                        logger.log(5, '[L{}] End convergence section.'.format(
                                i+1))
                    if section_convergence:
                        match = re.match(
                            '\s(Maximum|RMS)\s+(Force|Displacement)\s+({0})\s+'
                            '({0})\s+(YES|NO)'.format(
                                co.RE_FLOAT), line)
                        if match:
                            current_structure.props['{} {}'.format(
                                    match.group(1), match.group(2))] = \
                                (float(match.group(3)),
                                 float(match.group(4)), match.group(5))
                    if 'Converged?' in line:
                        section_convergence = True
                        logger.log(5, '[L{}] Start convergence section.'.format(
                                i+1))
                    # Look for input coords.
                    if coords_type == 'input' or coords_type == 'both':
                        # End of input coords for a given structure.
                        if section_coords_input and 'Distance matrix' in line:
                            section_coords_input = False
                            logger.log(5, '[L{}] End input coordinates section '
                                       '({} atoms).'.format(
                                    i+1, count_atom))
                        # Add atoms and coords to structure.
                        if section_coords_input:
                            match = re.match(
                                '\s+(\d+)\s+(\d+)\s+\d+\s+({0})\s+({0})\s+'
                                '({0})'.format(
                                    co.RE_FLOAT), line)
                            if match:
                                count_atom += 1
                                try:
                                    current_atom = current_structure.atoms[
                                        int(match.group(1))-1]
                                except IndexError:
                                    current_structure.atoms.append(Atom())
                                    current_atom = current_structure.atoms[-1]
                                if current_atom.atomic_num:
                                    assert current_atom.atomic_num == \
                                        int(match.group(2)), \
                                        ("[L{}] Atomic numbers don't match "
                                         "(current != existing) "
                                         "({} != {}).".format(
                                                i+1, int(match.group(2)),
                                                current_atom.atomic_num))
                                else:
                                    current_atom.atomic_num = \
                                        int(match.group(2))
                                current_atom.coords_type = 'input'
                                current_atom.x = float(match.group(3))
                                current_atom.y = float(match.group(4))
                                current_atom.z = float(match.group(5))
                        # Start of input coords for a given structure.
                        if not section_coords_input and \
                                'Input orientation:' in line:
                            section_coords_input = True
                            count_atom = 0
                            logger.log(5, '[L{}] Start input coordinates '
                                       'section.'.format(i+1))
                    # Look for standard coords.
                    if coords_type == 'standard' or coords_type == 'both':
                        # End of coordinates for a given structure.
                        if section_coords_standard and \
                                ('Rotational constants' in line or
                                 'Leave Link' in line):
                            section_coords_standard = False
                            logger.log(5, '[L{}] End standard coordinates '
                                       'section ({} atoms).'.format(
                                    i+1, count_atom))
                        # Grab coords for each atom. Add atoms to the structure.
                        if section_coords_standard:
                            match = re.match('\s+(\d+)\s+(\d+)\s+\d+\s+({0})\s+'
                                             '({0})\s+({0})'.format(
                                    co.RE_FLOAT), line)
                            if match:
                                count_atom += 1
                                try:
                                    current_atom = current_structure.atoms[
                                        int(match.group(1))-1]
                                except IndexError:
                                    current_structure.atoms.append(Atom())
                                    current_atom = current_structure.atoms[-1]
                                if current_atom.atomic_num:
                                    assert current_atom.atomic_num == int(
                                        match.group(2)), \
                                        ("[L{}] Atomic numbers don't match "
                                         "(current != existing) "
                                         "({} != {}).".format(
                                            i+1, int(match.group(2)),
                                            current_atom.atomic_num))
                                else:
                                    current_atom.atomic_num = int(
                                        match.group(2))
                                current_atom.coords_type = 'standard'
                                current_atom.x = float(match.group(3))
                                current_atom.y = float(match.group(4))
                                current_atom.z = float(match.group(5))
                        # Start of standard coords.
                        if not section_coords_standard and \
                                'Standard orientation' in line:
                            section_coords_standard = True
                            count_atom = 0
                            logger.log(5, '[L{}] Start standard coordinates '
                                       'section.'.format(i+1))
        return structures

def conv_sch_str(sch_struct):
    """
    Converts a schrodinger.structure object to my own structure object.
    Sort of pointless. Probably remove soon.
    """
    my_struct = Structure()
    my_struct.props.update(sch_struct.property)
    for sch_atom in sch_struct.atom:
        my_atom = Atom()
        my_struct.atoms.append(my_atom)
        my_atom.atom_type = sch_atom.atom_type
        my_atom.atom_type_name = sch_atom.atom_type_name
        my_atom.atomic_num = sch_atom.atomic_number
        my_atom.bonded_atom_indices = \
            [x.index for x in sch_atom.bonded_atoms]
        my_atom.element = sch_atom.element
        my_atom.index = sch_atom.index
        my_atom.partial_charge = sch_atom.partial_charge
        my_atom.x, my_atom.y, my_atom.z = sch_atom.x, sch_atom.y, sch_atom.z
        my_atom.props.update(sch_atom.property)
    for sch_bond in sch_struct.bond:
        my_bond = Bond()
        my_struct.bonds.append(my_bond)
        my_bond.atom_nums = [sch_bond.atom1, sch_bond.atom2]
        my_bond.order = sch_bond.order
        my_bond.value = sch_bond.length
    return my_struct

class SchrodingerFile(File):
    """
    Parent class used for all Schrodinger files.
    """
    def conv_sch_str(self, sch_struct):
        """
        Converts a schrodinger.structure object to my own structure object.
        Sort of pointless. Probably remove soon.
        """
        my_struct = Structure()
        my_struct.props.update(sch_struct.property)
        for sch_atom in sch_struct.atom:
            my_atom = Atom()
            my_struct.atoms.append(my_atom)
            my_atom.atom_type = sch_atom.atom_type
            my_atom.atom_type_name = sch_atom.atom_type_name
            my_atom.atomic_num = sch_atom.atomic_number
            my_atom.bonded_atom_indices = \
                [x.index for x in sch_atom.bonded_atoms]
            my_atom.element = sch_atom.element
            my_atom.index = sch_atom.index
            my_atom.partial_charge = sch_atom.partial_charge
            my_atom.x, my_atom.y, my_atom.z = sch_atom.x, sch_atom.y, sch_atom.z
            my_atom.props.update(sch_atom.property)
        for sch_bond in sch_struct.bond:
            my_bond = Bond()
            my_struct.bonds.append(my_bond)
            my_bond.atom_nums = [sch_bond.atom1, sch_bond.atom2]
            my_bond.order = sch_bond.order
            my_bond.value = sch_bond.length
        return my_struct

class JaguarIn(SchrodingerFile):
    """
    Used to retrieve data from Jaguar .in files.
    """
    def __init__(self, path):
        super(JaguarIn, self).__init__(path)
        self._structures = None
        self._hessian = None
        self._empty_atoms = None
        self._lines = None
    @property
    def hessian(self):
        """
        Reads the Hessian from a Jaguar .in.

        Automatically removes Hessian elements corresponding to dummy atoms.
        """
        if self._hessian is None:
            num = len(self.structures[0].atoms) + len(self._empty_atoms)
            logger.log(5,
                       '  -- {} has {} atoms and {} dummy atoms.'.format(
                    self.filename,
                    len(self.structures[0].atoms),
                    len(self._empty_atoms)))
            assert num != 0, \
                'Zero atoms found when loading Hessian from {}!'.format(
                self.path)
            hessian = np.zeros([num * 3, num * 3], dtype=float)
            logger.log(5, '  -- Created {} Hessian matrix (including dummy '
                       'atoms).'.format(hessian.shape))
            with open(self.path, 'r') as f:
                section_hess = False
                for line in f:
                    if section_hess and line.startswith('&'):
                        section_hess = False
                        hessian += np.tril(hessian, -1).T
                    if section_hess:
                        cols = line.split()
                        if len(cols) == 1:
                            hess_col = int(cols[0])
                        elif len(cols) > 1:
                            hess_row = int(cols[0])
                            for i, hess_ele in enumerate(cols[1:]):
                                hessian[hess_row - 1, i + hess_col - 1] = \
                                    float(hess_ele)
                    if '&hess' in line:
                        section_hess = True
            for atom in self._empty_atoms:
                logger.log(1, '>>> _empty_atom {}: {}'.format(atom.index, atom))
            # Figure out the indices of the dummy atoms.
            dummy_indices = []
            for atom in self._empty_atoms:
                logger.log(1, '>>> atom.index: {}'.format(atom.index))
                index = (atom.index - 1) * 3
                dummy_indices.append(index)
                dummy_indices.append(index + 1)
                dummy_indices.append(index + 2)
            logger.log(1, '>>> dummy_indices: {}'.format(dummy_indices))
            # Delete these rows and columns.
            logger.log(1, '>>> hessian.shape: {}'.format(hessian.shape))
            logger.log(1, '>>> hessian:\n{}'.format(hessian))
            hessian = np.delete(hessian, dummy_indices, 0)
            hessian = np.delete(hessian, dummy_indices, 1)
            logger.log(1, '>>> hessian:\n{}'.format(hessian))
            logger.log(5, '  -- Created {} Hessian matrix (w/o dummy '
                       'atoms).'.format(hessian.shape))
            self._hessian = hessian * co.HESSIAN_CONVERSION
            logger.log(1, '>>> hessian.shape: {}'.format(hessian.shape))
        return self._hessian
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            sch_ob = jag_in.read(self.path)
            sch_struct = sch_ob.getStructure()
            structures = [self.conv_sch_str(sch_struct)]
            logger.log(5, '  -- Imported {} structure(s).'.format(
                    len(structures)))
            # This area is sketch. I added it so I could use Hessian data
            # generated from a Jaguar calculation that had a dummy atom.
            # No gaurantees this will always work.
            for i, structure in enumerate(structures):
                empty_atoms = []
                for atom in structure.atoms:
                    logger.log(1, '>>> atom {}: {}'.format(atom.index, atom))
                    if atom.element == '':
                        empty_atoms.append(atom)
                for atom in empty_atoms:
                    structure.atoms.remove(atom)
                if empty_atoms:
                    logger.log(5, 'Structure {}: {} empty atoms '
                               'removed.'.format(i + 1, len(empty_atoms)))
            self._empty_atoms = empty_atoms
            self._structures = structures
        return self._structures
    def gen_lines(self):
        """
        Attempts to figure out the lines of itself.

        Since it'd be difficult, the written version will be missing much
        of the data in the original. Maybe there's something in the
        Schrodinger API for that.

        However, I do want this to include the ability to write out an
        atomic section with the ESP data that we'd want.
        """
        lines = []
        mae_name = None
        lines.append('MAEFILE: {}'.format(mae_name))
        lines.append('&gen')
        lines.append('&')
        lines.append('&zmat')
        # Just use the 1st structure. I don't imagine a Jaguar input file
        # ever containing more than one structure.
        struct = self.structures[0]
        lines.extend(struct.format_coords(format='gauss'))
        lines.append('&')
        return lines

class JaguarOut(File):
    """
    Used to retrieve data from Schrodinger Jaguar .out files.
    """
    def __init__(self, path):
        super(JaguarOut, self).__init__(path)
        self._structures = None
        self._eigenvalues = None
        self._eigenvectors = None
        self._frequencies = None
        self._dummy_atom_eigenvector_indices = None
        # self._force_constants = None
    @property
    def structures(self):
        if self._structures is None:
            self.import_file()
        return self._structures
    @property
    def eigenvalues(self):
        if self._eigenvalues is None:
            self.import_file()
        return self._eigenvalues
    @property
    def eigenvectors(self):
        if self._eigenvectors is None:
            self.import_file()
        return self._eigenvectors
    @property
    def frequencies(self):
        if self._frequencies is None:
            self.import_file()
        return self._frequencies
    @property
    def dummy_atom_eigenvector_indices(self):
        if self._dummy_atom_eigenvector_indices is None:
            self.import_file()
        return self._dummy_atom_eigenvector_indices
    def import_file(self):
        logger.log(10, 'READING: {}'.format(self.filename))
        frequencies = []
        force_constants = []
        eigenvectors = []
        structures = []
        with open(self.path, 'r') as f:
            section_geometry = False
            section_eigenvalues = False
            section_eigenvectors = False
            for i, line in enumerate(f):
                if section_geometry:
                    cols = line.split()
                    if len(cols) == 0:
                        section_geometry = False
                        structures.append(current_structure)
                        continue
                    elif len(cols) == 1:
                        pass
                    else:
                        match = re.match(
                            '\s+([\d\w]+)\s+({0})\s+({0})\s+({0})'.format(
                                co.RE_FLOAT), line)
                        if match != None:
                            current_atom = Atom()
                            current_atom.element = match.group(1).translate(
                                None, digits)
                            current_atom.x = float(match.group(2))
                            current_atom.y = float(match.group(3))
                            current_atom.z = float(match.group(4))
                            current_structure.atoms.append(current_atom)
                            logger.log(0,
                                       '{0:<3}{1:>12.6f}{2:>12.6f}'
                                       '{3:>12.6f}'.format(
                                    current_atom.element, current_atom.x,
                                    current_atom.y, current_atom.z))
                if 'geometry:' in line:
                    section_geometry = True
                    current_structure = Structure()
                    logger.log(5, '[L{}] Located geometry.'.format(i + 1))
                if 'Number of imaginary frequencies' in line or \
                        'Writing vibrational' in line or \
                        'Thermochemical properties at' in line:
                    section_eigenvalues = False
                if section_eigenvectors is True:
                    cols = line.split()
                    if len(cols) == 0:
                        section_eigenvectors = False
                        eigenvectors.extend(temp_eigenvectors)
                        continue
                    else:
                        for i, x in enumerate(cols[2:]):
                            if not len(temp_eigenvectors) > i:
                                temp_eigenvectors.append([])
                            temp_eigenvectors[i].append(float(x))
                if section_eigenvalues is True and \
                        section_eigenvectors is False:
                    if 'frequencies' in line:
                        cols = line.split()
                        frequencies.extend(map(float, cols[1:]))
                    if 'force const' in line:
                        cols = line.split()
                        force_constants.extend(map(float, cols[2:]))
                        section_eigenvectors = True
                        temp_eigenvectors = [[]]
                if 'normal modes in' in line:
                    section_eigenvalues = True
        logger.log(1, '>>> len(frequencies): {}'.format(len(frequencies)))
        logger.log(1, '>>> frequencies:\n{}'.format(frequencies))
        # logger.log(1, '>>> frequencies:\n{}'.format(
        #         [x / co.FORCE_CONVERSION for x in frequencies]))
        # logger.log(1, '>>> frequencies:\n{}'.format(
        #         [x * 4.55633e-6 for x in frequencies]))
        # logger.log(1, '>>> frequencies:\n{}'.format(
        #         [x * 1.23981e-4 for x in frequencies]))
        # logger.log(1, '>>> frequencies:\n{}'.format(
        #         [x / 219474.6305 for x in frequencies]))
        eigenvalues = [- fc / co.FORCE_CONVERSION if f < 0 else
                         fc / co.FORCE_CONVERSION
                         for fc, f in zip(force_constants, frequencies)]
        logger.log(1, '>>> eigenvalues:\n{}'.format(eigenvalues))
        # Remove eigenvector components related to dummy atoms.
        # Find the index of the atoms that are dummies.
        dummy_atom_indices = []
        for i, atom in enumerate(structures[-1].atoms):
            if atom.is_dummy:
                dummy_atom_indices.append(i)
        logger.log(10, '  -- Located {} dummy atoms.'.format(len(dummy_atom_indices)))
        # Correlate those indices to the rows in the cartesian eigenvector.
        dummy_atom_eigenvector_indices = []
        for dummy_atom_index in dummy_atom_indices:
            start = dummy_atom_index * 3
            dummy_atom_eigenvector_indices.append(start)
            dummy_atom_eigenvector_indices.append(start + 1)
            dummy_atom_eigenvector_indices.append(start + 2)
        new_eigenvectors = []
        # Create new eigenvectors without the rows corresponding to the
        # dummy atoms.
        for eigenvector in eigenvectors:
            new_eigenvectors.append([])
            for i, eigenvector_row in enumerate(eigenvector):
                if i not in dummy_atom_eigenvector_indices:
                    new_eigenvectors[-1].append(eigenvector_row)
        # Replace old eigenvectors with new where dummy atoms aren't included.
        eigenvectors = np.array(new_eigenvectors)
        self._dummy_atom_eigenvector_indices = dummy_atom_eigenvector_indices
        self._structures = structures
        self._eigenvalues = np.array(eigenvalues)
        self._eigenvectors = np.array(eigenvectors)
        self._frequencies = np.array(frequencies)
        # self._force_constants = np.array(force_constants)
        logger.log(5, '  -- Read {} structures'.format(
                len(self.structures)))
        logger.log(5, '  -- Read {} frequencies.'.format(
                len(self.frequencies)))
        logger.log(5, '  -- Read {} eigenvalues.'.format(
                len(self.eigenvalues)))
        logger.log(5, '  -- Read {} eigenvectors.'.format(
                self.eigenvectors.shape))
        # num_atoms = len(structures[-1].atoms)
        # logger.log(5,
        #            '  -- ({}, {}) eigenvectors expected for linear '
        #            'molecule.'.format(
        #         num_atoms * 3 - 5, num_atoms * 3))
        # logger.log(5, '  -- ({}, {}) eigenvectors expected for nonlinear '
        #            'molecule.'.format(
        #         num_atoms * 3 - 6, num_atoms * 3))

class Mae(SchrodingerFile):
    """
    Used to retrieve and work with data from Schrodinger .mae files.
    """
    def __init__(self, path):
        super(Mae, self).__init__(path)
        self._index_output_mae = None
        self._index_output_mmo = None
        self._structures = None
        self.commands = None
        # Strings for keeping track of this file and output files.
        self.name = os.path.splitext(self.filename)[0]
        self.name_com = self.name + '.q2mm.com'
        self.name_log = self.name + '.q2mm.log'
        self.name_mae = self.name + '.q2mm.mae'
        self.name_mmo = self.name + '.q2mm.mmo'
        self.name_out = self.name + '.q2mm.out'
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            # It would be great if we could leave this as an iter.
            try:
                sch_structs = list(sch_str.StructureReader(self.path))
            except:
                logger.warning('Error reading {}.'.format(self.path))
                raise
            self._structures = [self.conv_sch_str(sch_struct)
                                for sch_struct in sch_structs]
            logger.log(5, '  -- Imported {} structure(s).'.format(
                    len(self._structures)))
        return self._structures
    def get_com_opts(self):
        """
        Takes the users arguments from calculate (ex. mb, me, etc.) and
        determines what has to be written to the .com file in order to
        generate the requested data using MacroModel.

        Returns
        -------
        dictionary of options used when writing a .com file
        """
        com_opts = {
            'freq': False,
            'opt': False,
            'opt_mmo': False,
            'sp': False,
            'sp_mmo': False,
            'strs': False,
            'tors': False}
        if len(self.structures) > 1:
            com_opts['strs'] = True
        if any(x in ['jb', 'ja', 'jt'] for x in self.commands):
            com_opts['sp_mmo'] = True
        if any(x in ['me', 'mea', 'mq', 'mqh', 'mqa', 'mgESP', 'mjESP'] 
            for x in self.commands):
            com_opts['sp'] = True
        # Command meig is depreciated.
        if any(x in ['mh', 'meig', 'mjeig', 'mgeig'] for x in self.commands):
            if com_opts['strs']:
                raise Exception(
                    "Can't obtain the Hessian from a Maestro file "
                    "containing multiple structures!\n"
                    "FILENAME: {}\n"
                    "COMMANDS:{}\n".format(
                        self.path, ' '.join(commands)))
            else:
                com_opts['freq'] = True
        if any(x in ['mb', 'ma', 'mt', 'meo', 'meao'] for x in self.commands):
            com_opts['opt'] = True
            com_opts['opt_mmo'] = True
        elif any(x in ['mb', 'ma', 'mt'] for x in self.commands):
            com_opts['opt'] = True
        if any(x in ['mt', 'jt'] for x in self.commands):
            com_opts['tors'] = True
        return com_opts
    def get_debg_opts(self, com_opts):
        """
        Determines what arguments are needed for the DEBG line used inside
        a MacroModel .com file.

        Returns
        -------
        list of integers
        """
        debg_opts = []
        debg_opts.append(57)
        # Leads to problems when an angle inside a torsion is ~ 0 or 180.
        # if com_opts['tors']:
        #     debg_opts.append(56)
        if com_opts['freq']:
            debg_opts.extend((210, 211))
        debg_opts.sort()
        debg_opts.insert(0, 'DEBG')
        while len(debg_opts) < 9:
            debg_opts.append(0)
        return debg_opts
    def write_com(self, sometext=None):
        """
        Writes the .com file with all the right arguments to generate
        the requested data.
        """
        # Setup new filename. User can add additional text.
        if sometext:
            pieces = self.name_com.split('.')
            pieces.insert(-1, sometext)
            self.name_com = '.'.join(pieces)
        # Even if the command file already exists, we still need to
        # determine these indices.
        self._index_output_mae = []
        self._index_output_mmo = []
        com_opts = self.get_com_opts()
        debg_opts = self.get_debg_opts(com_opts)
        com = '{}\n{}\n'.format(self.filename, self.name_mae)
        # Is this right? It seems to work, but looking back at this,
        # I'm not sure why we wouldn't always want to control using
        # MMOD. Also, that 2nd argument of MMOD only affects the color
        # of atoms. I don't think this needs to be included. At some
        # point, I am going to remove it and test everything to make
        # sure it's not essential.
        if debg_opts:
            com += co.COM_FORM.format(*debg_opts)
        else:
            com += co.COM_FORM.format('MMOD', 0, 1, 0, 0, 0, 0, 0, 0)
        # May want to turn on/off arg2 (continuum solvent).
        #com += co.COM_FORM.format('FFLD', 2, 0, 0, 0, 36.7, 0, 0, 0)
        com += co.COM_FORM.format('FFLD', 2, 0, 0, 0, 0, 0, 0, 0)
        # Also may want to turn on/off cutoffs using BDCO.
        ## We have noticed there are some oddities for electrostatic 
        ## interactions. In some cases "Residue-based cutoffs" are used which
        ## result inconsistent energy contributions. Inorder to keep this
        ## consistent between calculations in Q2MM and conformational seaching
        ## EXNB is used to set all vdW and electrostatic cutoffs to 99.0 
        ## ensuring all interactions are gathered. The seventh column is the
        ## cutoff for hydrogen bonds, and 4 is the default value.
        com += co.COM_FORM.format('EXNB', 0, 0, 0, 0, 99., 99., 4., 0)
        if com_opts['strs']:
            com += co.COM_FORM.format('BGIN', 0, 0, 0, 0, 0, 0, 0, 0)
        # Look into differences.
        com += co.COM_FORM.format('READ', -1, 0, 0, 0, 0, 0, 0, 0)
        if com_opts['sp'] or com_opts['sp_mmo'] or com_opts['freq']:
            com += co.COM_FORM.format('MINI', 9, 0, 0, 0, 0, 0, 0, 0)
            # self._index_output_mae.append('stupid_extra_structure')
            self._index_output_mae.append('pre')
        if com_opts['sp'] or com_opts['sp_mmo']:
            com += co.COM_FORM.format('ELST', 1, 0, 0, 0, 0, 0, 0, 0)
            self._index_output_mmo.append('pre')
            # Replaced by using a pointless MINI statement. For whatever
            # reason, that causes the .mmo file to be written without
            # needing this WRIT statement.
            # com += co.COM_FORM.format('WRIT', 0, 0, 0, 0, 0, 0, 0, 0)
            # self._index_output_mae.append('pre')
        if com_opts['freq']:
            # Now the WRIT is handled above.
            # com += co.COM_FORM.format('MINI', 9, 0, 0, 0, 0, 0, 0, 0)
            # self._index_output_mae.append('stupid_extra_structure')
            # What does arg1 as 3 even do?
            com += co.COM_FORM.format('RRHO', 3, 0, 0, 0, 0, 0, 0, 0)
            self._index_output_mae.append('hess')
        if com_opts['opt']:
            # Commented line was used in code from Per-Ola/Elaine.
            # com += co.COM_FORM.format('MINI', 9, 0, 50, 0, 0, 0, 0, 0)
            # TNCG has more risk of not converging, and may print NaN instead
            # of coordinates and forces to output.
            # arg1: 1 = PRCG, 9 = TNCG
            com += co.COM_FORM.format('MINI', 1, 0, 500, 0, 0, 0, 0, 0)
            self._index_output_mae.append('opt')
        if com_opts['opt_mmo']:
            com += co.COM_FORM.format('ELST', 1, 0, 0, 0, 0, 0, 0, 0)
            self._index_output_mmo.append('opt')
        if com_opts['strs']:
            com += co.COM_FORM.format('END', 0, 0, 0, 0, 0, 0, 0, 0)
        # If the file already exists, don't rewrite it.
        path_com = os.path.join(self.directory, self.name_com)
        if sometext and os.path.exists(path_com):
            logger.log(5, '  -- {} already exists. Skipping write.'.format(
                    self.name_com))
        else:
            with open(os.path.join(self.directory, self.name_com), 'w') as f:
                f.write(com)
            logger.log(5, 'WROTE: {}'.format(
                    os.path.join(self.name_com)))

    def run(self, max_fails=5, max_timeout=None, timeout=10, check_tokens=True):
        """
        Runs MacroModel .com files. This has to be more complicated than a
        simple subprocess command due to problems with Schrodinger tokens.
        This script checks the available tokens, and if there's not enough,
        waits to run MacroModel until there are.

        Arguments
        ---------
        max_timeout : int
                      Maximum number of attempts to look for Schrodinger
                      license tokens before giving up.
        max_fails : int
                    Maximum number of times the job can fail.
        timeout : float
                  Time waited in between lookups of Schrodinger license
                  tokens.
        """
        #print("Run " + str(self.filename) + " with commands:" + str(self.commands))
        current_directory = os.getcwd()
        os.chdir(self.directory)
        current_timeout = 0
        current_fails = 0
        licenses_available = False
        if check_tokens is True:
            logger.log(5, "  -- Checking Schrodinger tokens.")
            while True:
                token_string = sp.check_output(
                    '$SCHRODINGER/utilities/licutil -available', shell=True)
                if (sys.version_info > (3, 0)):
                  token_string = token_string.decode("utf-8")
                if 'SUITE' not in token_string:
                    licenses_available = True
                    break
                suite_tokens = co.LIC_SUITE.search(token_string)
                macro_tokens = co.LIC_MACRO.search(token_string)
                #suite_tokens = re.search(co.LIC_SUITE, token_string)
                #macro_tokens = re.search(co.LIC_MACRO, token_string)
                if not suite_tokens or not macro_tokens:
                    raise Exception(
                        'The command "$SCHRODINGER/utilities/licutil '
                        '-available" is not working with the current '
                        'regex in calculate.py.\nOUTPUT:\n{}'.format(
                            token_string))
                suite_tokens = int(suite_tokens.group(1))
                macro_tokens = int(macro_tokens.group(1))
                if suite_tokens > co.MIN_SUITE_TOKENS and \
                        macro_tokens > co.MIN_MACRO_TOKENS:
                    licenses_available = True
                    break
                else:
                    if max_timeout is not None and \
                            current_timeout > max_timeout:
                        pretty_timeout(
                            current_timeout, suite_tokens,
                            macro_tokens, end=True, name_com=self.name_com)
                        raise Exception(
                            "Not enough tokens to run {}. Waited {} seconds "
                            "before giving up.".format(
                                self.name_com, current_timeout))
                    pretty_timeout(current_timeout, suite_tokens, macro_tokens,
                                   name_com=self.name_com)
                    current_timeout += timeout
                    time.sleep(timeout)
        else:
            licenses_available = True
        if licenses_available:
            while True:
                try:
                    logger.log(5, 'RUNNING: {}'.format(self.name_com))
                    sp.check_output(
                        '$SCHRODINGER/bmin -WAIT {}'.format(
                            os.path.splitext(self.name_com)[0]), shell=True)
                    break
                except sp.CalledProcessError:
                    logger.warning('Call to MacroModel failed and I have no '
                                   'idea why!')
                    current_fails += 1
                    if current_fails < max_fails:
                        time.sleep(timeout)
                        continue
                    else:
                        raise
        os.chdir(current_directory)

def pretty_timeout(current_timeout, macro_tokens, suite_tokens, end=False,
                   level=10, name_com=None):
    """
    Logs information about the wait for Schrodinger tokens.

    Arguments
    ---------
    current_timeout : int
                      Number of times waited for Schrodinger tokens.
    macro_tokens : int
                   Current number of available MacroModel tokens.
    suite_tokens : int
                   Current number of available Schrodinger Suite tokens.
    end : bool
          If True, adds a pretty ending border to all these logs.
    level : int
            Logging level of the pretty messages.
    """
    if current_timeout == 0:
        if name_com:
            logger.warning('  -- Waiting on tokens to run {}.'.format(
                    name_com))
        logger.log(level,
                   '--' + ' (s) '.center(8, '-') +
                   '--' + ' {} '.format(co.LABEL_SUITE).center(17, '-') +
                   '--' + ' {} '.format(co.LABEL_MACRO).center(17, '-') +
                   '--')
    logger.log(level, '  {:^8d}  {:^17d}  {:^17d}'.format(
            current_timeout, macro_tokens, suite_tokens))
    if end is True:
        logger.log(level, '-' * 50)

class MacroModelLog(File):
    """
    Used to retrieve data from MacroModel log files.
    """
    def __init__(self, path):
        super(MacroModelLog, self).__init__(path)
        self._hessian = None
    @property
    def hessian(self):
        if self._hessian is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            with open(self.path, 'r') as f:
                lines = f.read()
            num_atoms = int(re.search('Read\s+(\d+)\s+atoms.', lines).group(1))
            logger.log(5, '  -- Read {} atoms.'.format(num_atoms))

            hessian = np.zeros([num_atoms * 3, num_atoms * 3], dtype=float)
            logger.log(5, '  -- Creating {} Hessian matrix.'.format(hessian.shape))
            words = lines.split()
            section_hessian = False
            start_row = False
            start_col = False
            for i, word in enumerate(words):
                # 1. Start of Hessian section.
                if word == 'Mass-weighted':
                    section_hessian = True
                    continue
                # 5. End of Hessian. Add last row of Hessian and break.
                if word == 'Eigenvalues:':
                    for col_num, element in zip(col_nums, elements):
                        hessian[row_num - 1, col_num - 1] = element
                    section_hessian = False
                    break
                # 4. End of a Hessian row. Add to matrix and reset.
                if section_hessian and start_col and word == 'Element':
                    for col_num, element in zip(col_nums, elements):
                        hessian[row_num - 1, col_num - 1] = element
                    start_col = False
                    start_row = True
                    row_num = int(words[i + 1])
                    col_nums = []
                    elements = []
                    continue
                # 2. Start of a Hessian row.
                if section_hessian and word == 'Element':
                    row_num = int(words[i + 1])
                    col_nums = []
                    elements = []
                    start_row = True
                    continue
                # 3. Okay, made it through the row number. Now look for columns
                #    and elements.
                if section_hessian and start_row and word == ':':
                    start_row = False
                    start_col = True
                    continue
                if section_hessian and start_col and '.' not in word and \
                        word != 'NaN':
                    col_nums.append(int(word))
                    continue
                if section_hessian and start_col and '.' in word or \
                        word == 'NaN':
                    elements.append(float(word))
                    continue
            self._hessian = hessian
            logger.log(5, '  -- Creating {} Hessian matrix.'.format(hessian.shape))
        return self._hessian

class MacroModel(File):
    """
    Extracts data from MacroModel .mmo files.
    """
    def __init__(self, path):
        super(MacroModel, self).__init__(path)
        self._structures = None
    @property
    def structures(self):
        if self._structures is None:
            logger.log(10, 'READING: {}'.format(self.filename))
            self._structures = []
            with open(self.path, 'r') as f:
                count_current = 0
                count_input = 0
                count_structure = 0
                count_previous = 0
                section = None
                for line in f:
                # This would probably be better as a function in the structure
                # class but I wanted this as upstream as possible so I didn't
                # have to worry about other coding issues. The MMO file lists
                # the bonds, angles, and torsions in some order that I am unsure
                # of. It seems consistent with the same filename but with two
                # files with the exact same structure the ordering is off. This
                # reorders the lists before being added to the structure class.
                    if 'Input filename' in line:
                        count_input += 1
                    if 'Input Structure Name' in line:
                        count_structure += 1
                    count_previous = count_current
                    # Sometimes only one of the above ("Input filename" and
                    # "Input Structure Name") is used, sometimes both are used.
                    # count_current will make sure you catch both.
                    count_current = max(count_input, count_structure)
                    # If these don't match, then we reached the end of a
                    # structure.
                    if count_current != count_previous:
                        bonds = []
                        angles = []
                        torsions = []
                        current_structure = Structure()
                        self._structures.append(current_structure)
                    # For each structure we come across, look for sections that
                    # we are interested in: those pertaining to bonds, angles,
                    # and torsions. Of course more could be added. We set the
                    # section to None to mark the end of a section, and we leave
                    # it None for parts of the file we don't care about.
                    if 'BOND LENGTHS AND STRETCH ENERGIES' in line:
                        section = 'bond'
                    if 'ANGLES, BEND AND STRETCH BEND ENERGIES' in line:
                        section = 'angle'
                    if 'BEND-BEND ANGLES AND ENERGIES' in line:
                        section = None
                    if 'DIHEDRAL ANGLES AND TORSIONAL ENERGIES' in line:
                        section = 'torsion'
                    if 'DIHEDRAL ANGLES AND TORSIONAL CROSS-TERMS' in line:
                        section = None
                    if section == 'bond':
                        bond = self.read_line_for_bond(line)
                        if bond is not None:
                            #current_structure.bonds.append(bond)
                            bonds.append(bond)
                    if section == 'angle':
                        angle = self.read_line_for_angle(line)
                        if angle is not None:
                            #current_structure.angles.append(angle)
                            angles.append(angle)
                    if section == 'torsion':
                        torsion = self.read_line_for_torsion(line)
                        if torsion is not None:
                            #current_structure.torsions.append(torsion)
                            torsions.append(torsion)
                    if 'Connection Table' in line:
                        # Sort the bonds, angles, and torsions before the start
                        # of a new structure
                        if bonds:
                            bonds.sort(key = lambda x: (x.atom_nums[0], x.atom_nums[1]))
                            current_structure.bonds.extend(bonds)
                        if angles:
                            angles.sort(key = lambda x: (x.atom_nums[1], x.atom_nums[0], x.atom_nums[2]))
                            current_structure.angles.extend(angles)
                        if torsions:
                            torsions.sort(key = lambda x: (x.atom_nums[1], x.atom_nums[2], x.atom_nums[0], x.atom_nums[3]))
                            current_structure.torsions.extend(torsions)
            logger.log(5, '  -- Imported {} structure(s).'.format(
                    len(self._structures)))
        return self._structures
    def read_line_for_bond(self, line):
        match = co.RE_BOND.match(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(2)]]
            value = float(match.group(3))
            comment = match.group(4).strip()
            ff_row = int(match.group(5))
            return Bond(atom_nums=atom_nums, comment=comment, value=value,
                        ff_row=ff_row)
        else:
            return None
    def read_line_for_angle(self, line):
        match = co.RE_ANGLE.match(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(2),
                                  match.group(3)]]
            # Reorder the terminal atoms so that the lower index atom is first.
            if atom_nums[0] > atom_nums[2]:
                atom_nums = [atom_nums[2],atom_nums[1],atom_nums[0]]
            value = float(match.group(4))
            comment = match.group(5).strip()
            ff_row = int(match.group(6))
            return Angle(atom_nums=atom_nums, comment=comment, value=value,
                         ff_row=ff_row)
        else:
            return None
    def read_line_for_torsion(self, line):
        match = co.RE_TORSION.match(line)
        if match:
            atom_nums = [int(x) for x in [match.group(1), match.group(2),
                                  match.group(3), match.group(4)]]
            if atom_nums[1] > atom_nums[2]:
                atom_nums = [atom_nums[3],
                             atom_nums[2],
                             atom_nums[1],
                             atom_nums[0]]
            value = float(match.group(5))
            comment = match.group(6).strip()
            ff_row = int(match.group(7))
            return Torsion(atom_nums=atom_nums, comment=comment, value=value,
                           ff_row=ff_row)
        else:
            return None

def select_structures(structures, indices, label):
        """
        Returns a list of structures where the index matches the label. This
        is used with the structures in the class MacroModel (.mmo's) and Mae
        (.mae's of course).

        Basically, you're not sure what structures appear in these files if the
        files were generated using calculate.py and the .com files it writes.
        Fear not! calculate.py keeps track of that for you (using indices) and
        knows which structures to use.

        indices - A list of strings (labels).
        label   - A string. Possible strings include:
                      'opt', 'pre', 'hess' (.mae only), and
                      'stupid_extra_structure'
        """
        selected = []
        idx_iter = iter(indices)
        for str_num, struct in enumerate(structures):
            try:
                idx_curr = next(idx_iter)
            except StopIteration:
                idx_iter = iter(indices)
                idx_curr = next(idx_iter)
            if idx_curr == label:
                selected.append((str_num, struct))
        return selected

# This could use some documentation. Looks pretty though.
def geo_from_points(*args):
    x1 = args[0][0]
    y1 = args[0][1]
    z1 = args[0][2]
    x2 = args[1][0]
    y2 = args[1][1]
    z2 = args[1][2]
    if len(args) == 2:
        bond = math.sqrt((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)
        return float(bond)
    x3 = args[2][0]
    y3 = args[2][1]
    z3 = args[2][2]
    if len(args) == 3:
        dist_21 = math.sqrt((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)
        dist_23 = math.sqrt((x2 - x3)**2 + (y2 - y3)**2 + (z2 - z3)**2)
        dist_13 = math.sqrt((x1 - x3)**2 + (y1 - y3)**2 + (z1 - z3)**2)
        angle = math.acos((dist_21**2 + dist_23**2 - dist_13**2) /
                          (2*dist_21*dist_23))
        angle = math.degrees(angle)
        return float(angle)
    x4 = args[3][0]
    y4 = args[3][1]
    z4 = args[3][2]
    if len(args) == 4:
        vect_21 = [x2 - x1, y2 - y1, z2 - z1]
        vect_32 = [x3 - x2, y3 - y2, z3 - z2]
        vect_43 = [x4 - x3, y4 - y3, z4 - z3]
        x_ab = np.cross(vect_21,vect_32)
        x_bc = np.cross(vect_32,vect_43)
        norm_ab = x_ab/(math.sqrt(x_ab[0]**2 + x_ab[1]**2 + x_ab[2]**2))
        norm_bc = x_bc/(math.sqrt(x_bc[0]**2 + x_bc[1]**2 + x_bc[2]**2))
        mag_ab = math.sqrt(norm_ab[0]**2 + norm_ab[1]**2 + norm_ab[2]**2)
        mag_bc = math.sqrt(norm_bc[0]**2 + norm_bc[1]**2 + norm_bc[2]**2)
        angle = math.acos(np.dot(norm_ab, norm_bc)/(mag_ab * mag_bc))
        torsion = angle * (180/math.pi)
        return torsion

class Structure(object):
    """
    Data for a single structure/conformer/snapshot.
    """
    __slots__ = ['atoms', 'bonds', 'angles', 'torsions', 'hess', 'props']
    def __init__(self):
        self.atoms = []
        self.bonds = []
        self.angles = []
        self.torsions = []
        self.hess = None
        self.props = {}
    @property
    def coords(self):
        """
        Returns atomic coordinates as a list of lists.
        """
        return [atom.coords for atom in self.atoms]
    def format_coords(self, format='latex', indices_use_charge=None):
        """
        Returns a list of strings/lines to easily generate coordinates
        in various formats.

        latex  - Makes a LaTeX table.
        gauss  - Makes output that matches Gaussian's .com filse.
        jaguar - Just like Gaussian, but include the atom number after the
                 element name in the left column.
        """
        # Formatted for LaTeX.
        if format == 'latex':
            output = ['\\begin{tabular}{l S[table-format=3.6] '
                      'S[table-format=3.6] S[table-format=3.6]}']
            for i, atom in enumerate(self.atoms):
                if atom.element is None:
                    ele = co.MASSES.items()[atom.atomic_num - 1][0]
                else:
                    ele = atom.element
                output.append('{0}{1} & {2:3.6f} & {3:3.6f} & '
                              '{4:3.6f}\\\\'.format(
                        ele, i+1, atom.x, atom.y, atom.z))
            output.append('\\end{tabular}')
            return output
        # Formatted for Gaussian .com's.
        elif format == 'gauss':
            output = []
            for i, atom in enumerate(self.atoms):
                if atom.element is None:
                    ele = co.MASSES.items()[atom.atomic_num - 1][0]
                else:
                    ele = atom.element
                # Used only for a problem Eric experienced.
                # if ele == '': ele = 'Pd'
                if indices_use_charge:
                    if atom.index in indices_use_charge:
                        output.append(
                            ' {0:s}--{1:.5f}{2:>16.6f}{3:16.6f}'
                            '{4:16.6f}'.format(
                                ele, atom.partial_charge, atom.x,
                                atom.y, atom.z))
                    else:
                        output.append(' {0:<8s}{1:>16.6f}{2:>16.6f}{3:>16.6f}'.format(
                                ele, atom.x, atom.y, atom.z))
                else:
                    output.append(' {0:<8s}{1:>16.6f}{2:>16.6f}{3:>16.6f}'.format(
                            ele, atom.x, atom.y, atom.z))
            return output
        # Formatted for Jaguar.
        elif format == 'jaguar':
            output = []
            for i, atom in enumerate(self.atoms):
                if atom.element is None:
                    ele = co.MASSES.items()[atom.atomic_num - 1][0]
                else:
                    ele = atom.element
                # Used only for a problem Eric experienced.
                # if ele == '': ele = 'Pd'
                label = '{}{}'.format(ele, atom.index)
                output.append(' {0:<8s}{1:>16.6f}{2:>16.6f}{3:>16.6f}'.format(
                        label, atom.x, atom.y, atom.z))
            return output
    def select_stuff(self, typ, com_match=None):
        """
        A much simpler version of select_data. It would be nice if select_data
        was a wrapper around this function.
        """
        stuff = []
        for thing in getattr(self, typ):
            if (com_match and any(x in thing.comment for x in com_match)) or \
                    com_match is None:
                stuff.append(thing)
        return stuff
    def select_data(self, typ, com_match=None, **kwargs):
        """
        Selects bonds, angles, or torsions from the structure and returns them
        in the format used as data.

        typ       - 'bonds', 'angles', or 'torsions'.
        com_match - String or None. If None, just returns all of the selected
                    stuff (bonds, angles, or torsions). If a string, selects
                    only those that have this string in their comment.

                    In .mmo files, the comment corresponds to the substructures
                    name. This way, we only fit bonds, angles, and torsions that
                    directly depend on our parameters.
        """
        data = []
        logger.log(1, '>>> typ: {}'.format(typ))
        for thing in getattr(self, typ):
            if (com_match and any(x in thing.comment for x in com_match)) or \
                    com_match is None:
                datum = thing.as_data(**kwargs)
                # If it's a torsion we have problems.
                # Have to check whether an angle inside the torsion is near 0 or 180.
                if typ == 'torsions':
                    atom_nums = [datum.atm_1, datum.atm_2, datum.atm_3, datum.atm_4]
                    angle_atoms_1 = [atom_nums[0], atom_nums[1], atom_nums[2]]
                    angle_atoms_2 = [atom_nums[1], atom_nums[2], atom_nums[3]]
                    for angle in self.angles:
                        if set(angle.atom_nums) == set(angle_atoms_1):
                            angle_1 = angle.value
                            break
                    for angle in self.angles:
                        if set(angle.atom_nums) == set(angle_atoms_2):
                            angle_2 = angle.value
                            break
                    try:
                        logger.log(1, '>>> atom_nums: {}'.format(atom_nums))
                        logger.log(1, '>>> angle_1: {} / angle_2: {}'.format(
                                angle_1, angle_2))
                    except UnboundLocalError:
                        logger.error('>>> atom_nums: {}'.format(atom_nums))
                        logger.error(
                            '>>> angle_atoms_1: {}'.format(angle_atoms_1))
                        logger.error(
                            '>>> angle_atoms_2: {}'.format(angle_atoms_2))
                        if 'angle_1' not in locals():
                            logger.error("Can't identify angle_1!")
                        else:
                            logger.error(">>> angle_1: {}".format(angle_1))
                        if 'angle_2' not in locals():
                            logger.error("Can't identify angle_2!")
                        else:
                            logger.error(">>> angle_2: {}".format(angle_2))
                        logger.warning('WARNING: Using torsion anyway!')
                        data.append(datum)
                    if -20. < angle_1 < 20. or 160. < angle_1 < 200. or \
                            -20. < angle_2 < 20. or 160. < angle_2 < 200.:
                        logger.log(
                            1, '>>> angle_1 or angle_2 is too close to 0 or 180!')
                        pass
                    else:
                        data.append(datum)
                    # atom_coords = [x.coords for x in atoms]
                    # tor_1 = geo_from_points(
                    #     atom_coords[0], atom_coords[1], atom_coords[2])
                    # tor_2 = geo_from_points(
                    #     atom_coords[1], atom_coords[2], atom_coords[3])
                    # logger.log(1, '>>> tor_1: {} / tor_2: {}'.format(
                    #     tor_1, tor_2))
                    # if -5. < tor_1 < 5. or 175. < tor_1 < 185. or \
                    #         -5. < tor_2 < 5. or 175. < tor_2 < 185.:
                    #     logger.log(
                    #         1,
                    #         '>>> tor_1 or tor_2 is too close to 0 or 180!')
                    #     pass
                    # else:
                    #     data.append(datum)
                else:
                    data.append(datum)
        assert data, "No data actually retrieved!"
        return data
    def get_aliph_hyds(self):
        """
        Returns the atom numbers of aliphatic hydrogens. These hydrogens
        are always assigned a partial charge of zero in MacroModel
        calculations.

        This should be subclassed into something is MM3* specific.
        """
        aliph_hyds = []
        for atom in self.atoms:
            if 40 < atom.atom_type < 49:
                for bonded_atom_index in atom.bonded_atom_indices:
                    bonded_atom = self.atoms[bonded_atom_index - 1]
                    if bonded_atom.atom_type == 3:
                        aliph_hyds.append(atom)
        logger.log(5, '  -- {} aliphatic hydrogen(s).'.format(len(aliph_hyds)))
        return aliph_hyds
    def get_hyds(self):
        """
        Returns the atom numbers of any default MacroModel type hydrogens.

        This should be subclassed into something is MM3* specific.
        """
        hyds = []
        for atom in self.atoms:
            if 40 < atom.atom_type < 49:
                for bonded_atom_index in atom.bonded_atom_indices:
                    hyds.append(atom)
        logger.log(5, '  -- {} hydrogen(s).'.format(len(hyds)))
        return hyds
    def get_dummy_atom_indices(self):
        """
        Returns a list of integers where each integer corresponds to an atom
        that is a dummy atom.

        Returns
        -------
        list of integers
        """
        dummies = []
        for atom in self.atoms:
            if atom.is_dummy:
                logger.log(
                    10,'  -- Identified {} as a dummy atom.'.format(atom))
                dummies.append(atom.index)
        return dummies

class Atom(object):
    """
    Data class for a single atom.

    Really, some of this atom type stuff should perhaps be in a MM3*
    specific atom class.
    """
    __slots__ = ['atom_type', 'atom_type_name', 'atomic_num', 'atomic_mass',
                 'bonded_atom_indices', 'coords_type', '_element',
                 '_exact_mass', 'index', 'partial_charge', 'x', 'y', 'z',
                 'props']
    def __init__(self, atom_type=None, atom_type_name=None, atomic_num=None,
                 atomic_mass=None, bonded_atom_indices=None, coords=None,
                 coords_type=None, element=None, exact_mass=None, index=None,
                 partial_charge=None, x=None, y=None, z=None):
        self.atom_type = atom_type
        self.atom_type_name = atom_type_name
        self.atomic_num = atomic_num
        self.atomic_mass = atomic_mass
        self.bonded_atom_indices = bonded_atom_indices
        self.coords_type = coords_type
        self._element = element
        self._exact_mass = exact_mass
        self.index = index
        self.partial_charge = partial_charge
        self.x = x
        self.y = y
        self.z = z
        if coords:
            self.x = coords[0]
            self.y = coords[1]
            self.z = coords[2]
        self.props = {}
    def __repr__(self):
            return '{}[{},{},{}]'.format(
                self.atom_type_name, self.x, self.y, self.z)
    @property
    def coords(self):
        return [self.x, self.y, self.z]
    @coords.setter
    def coords(self, value):
        try:
            self.x = value[0]
            self.y = value[1]
            self.z = value[2]
        except TypeError:
            pass
    @property
    def element(self):
        if self._element is None:
            self._element = co.MASSES.items()[self.atomic_num - 1][0]
        return self._element
    @element.setter
    def element(self, value):
        self._element = value
    @property
    def exact_mass(self):
        if self._exact_mass is None:
            self._exact_mass = co.MASSES[self.element]
        return self._exact_mass
    @exact_mass.setter
    def exact_mass(self, value):
        self._exact_mass = value
    # I have no idea if these atom types are actually correct.
    # Really, the user should specify custom atom types, such as dummies, in a
    # configuration file somewhere.
    @property
    def is_dummy(self):
        """
        Return True if self is a dummy atom, else return False.

        Returns
        -------
        bool
        """
        # I think 61 is the default dummy atom type in a Schrodinger atom.typ
        # file.
        # Okay, so maybe it's not. Anyway, Tony added an atom type 205 for
        # dummies. It'd be really great if we all used the same atom.typ file
        # someday.
        # Could add in a check for the atom_type number. I removed it.
        if self.atom_type_name == 'Du' or \
                self.element == 'X' or \
                self.atomic_num == -2:
            return True
        else:
            return False

class Bond(object):
    """
    Data class for a single bond.
    """
    __slots__ = ['atom_nums', 'comment', 'order', 'value', 'ff_row']
    def __init__(self, atom_nums=None, comment=None, order=None, value=None,
                 ff_row=None):
        self.atom_nums = atom_nums
        self.comment = comment
        self.order = order
        self.value = value
        self.ff_row = ff_row
    def __repr__(self):
        return '{}[{}]({})'.format(
            self.__class__.__name__, '-'.join(
                map(str, self.atom_nums)), self.value)
    def as_data(self, **kwargs):
        # Sort of silly to have all this stuff about angles and
        # torsions in here, but they both inherit from this class.
        # I suppose it'd make more sense to create a structural
        # element class that these all inherit from.
        # Warning that I recently changed these labels, and that
        # may have consequences.
        if self.__class__.__name__.lower() == 'bond':
            typ = 'b'
        elif self.__class__.__name__.lower() == 'angle':
            typ = 'a'
        elif self.__class__.__name__.lower() == 'torsion':
            typ = 't'
        datum = datatypes.Datum(val=self.value, typ=typ,ff_row=self.ff_row)
        for i, atom_num in enumerate(self.atom_nums):
            setattr(datum, 'atm_{}'.format(i+1), atom_num)
        for k, v in kwargs.items():
            setattr(datum, k, v)
        return datum

class Angle(Bond):
    """
    Data class for a single angle.
    """
    def __init__(self, atom_nums=None, comment=None, order=None, value=None,
                 ff_row=None):
        super(Angle, self).__init__(atom_nums, comment, order, value, ff_row)

class Torsion(Bond):
    """
    Data class for a single torsion.
    """
    def __init__(self, atom_nums=None, comment=None, order=None, value=None,
                 ff_row=None):
        super(Torsion, self).__init__(atom_nums, comment, order, value, ff_row)

def return_filetypes_parser():
    """
    Returns an argument parser for filetypes module.
    """
    parser = argparse.ArgumentParser(
        description=__doc__, formatter_class=RawTextHelpFormatter)
    parser.add_argument(
        '-i', '--input', type=str,
        help='Input filename.')
    parser.add_argument(
        '-o', '--output', type=str,
        help='Output filename.')
    parser.add_argument(
        '-p', '--print', action='store_true',
        help='Print coordinates for each structure.')
    parser.add_argument(
        '-n', '--num', type=int,
        help='Number of structures to display.')
    return parser

def detect_filetype(filename):
    path = os.path.abspath(filename)
    ext = os.path.splitext(path)[1]
    if ext == '.mae' or ext =='.maegz':
        file_ob = Mae(path)
    elif ext == '.log':
        file_ob = GaussLog(path)
        file_ob.read_out()
        # try:
        #     file_ob.read_archive()
        # except IndexError:
        #     pass
    elif ext == '.in':
        file_ob = JaguarIn(path)
    elif ext == '.out':
        file_ob = JaguarOut(path)
    else:
        raise Exception('Filetype not recognized.')
    return file_ob

def main(args):
    parser = return_filetypes_parser()
    opts = parser.parse_args(args)
    file_ob = detect_filetype(opts.input)
    if opts.print:
        if hasattr(file_ob, 'structures'):
            for i, structure in enumerate(file_ob.structures):
                print(' ' + ' STRUCTURE {} '.format(i + 1).center(56, '-'))
                output = structure.format_coords(format='gauss')
                for line in output:
                    print(line)
                if opts.num and i+1 == opts.num:
                    break
        if hasattr(file_ob, 'evals') and file_ob.evals:
            print('EIGENVALUES:')
            print(file_ob.evals)
        if hasattr(file_ob, 'evecs') and file_ob.evecs:
            print('EIGENVECTORS:')
            print(file_ob.evecs)
    if opts.output:
        file_ob.write(opts.output)

if __name__ == '__main__':
    import argparse
    import sys

    logging.config.dictConfig(co.LOG_SETTINGS)
    main(sys.argv[1:])