xyz2molcas2.py

#!/usr/bin/env python

import os
import argparse as ap

def get_atomnum(ele):
    list={
        "H":1,
        "C":6,
        "N":7,
        "O":8
        }
    try:
        return list[ele]
    except:
        print "Element "+ele+" not in list. Please add it!"


def error1(charge,num_protons,file):
    print "File: "+file+": Wrong charge specified for a closed shell molecule"
    print "Charge "+str(charge)+" given with total of "+str(num_protons)+" protons."
def error2(charge):
    print "Setting charge to " +str(charge-1)
    

def get_occu(file,charge):
    num_protons = 0
    xyz = open(file, 'r')
    for i in xyz.readlines():
        try:
            ele = i[0]
        except:
            print "empty line? "+i
        if ele in ["H","C","N","O"]:
            num_protons = num_protons + get_atomnum(ele)

    test = str(num_protons/2.0).split('.')
    if int(test[1])==5:
        if charge==0:
            error1(charge,num_protons,file)
            error2(charge)
            num_occ = num_protons/2 + 1
        else:
            num_occ = (num_protons-charge)/2
    elif int(test[1])==0:
        if charge==1 or charge==-1:
            error1(charge,num_protons,file)
            print "Please check your molecule."
            quit()
        elif charge==0:
            num_occ = num_protons/2

    return num_occ


parser = ap.ArgumentParser(description='LoProp job submitter 0.1')

parser.add_argument('--xyz', dest='xyzfile', metavar='XYZ_FILE',
                    help='''the name of the xyz input file that contains the
                            molecule.''')

parser.add_argument('--method', dest='method', metavar='METHOD',
                     default='DFT', choices=['DFT', 'HF'],
                     help='''Choose QM method. Valid choices are
                             %(choices)s. If you use DFT then you can also
                             specify the xc-functional using the
                             --xcfun option. [default: %(default)s]''')

parser.add_argument('--basis', dest='basis', metavar='BASISSET',
                     default='A-6-31PGP',
                     help='''Specify basis set. [default: %(default)s]''')

parser.add_argument('--xcfun', dest='xcfun', default='B3LYP',
                     metavar='XCFUN',
                     help='''Specify xc-functional. [default: %(default)s]''')

parser.add_argument('--bond-midpoints', dest='bonds', action='store_true',
                    default=False,
                    help='''Also calculate parameters on bond midpoints.
                            [default: %(default)s]''')

parser.add_argument('--charge', dest='charge', metavar='CHARGE',
                     default=0, type=int,
                     help='''Specify charge. [default: %(default)s]''')

parser.add_argument('--multi', dest='multi', metavar='MULTIPOLES',
                     default=2, type=int,
                     help='''Specify the max multipole moment. [default: %(default)s]''')

parser.add_argument('--occu', dest='occu', action='store_false',
                     default=True,
                     help='''Use occupated orbitals (OCCU) instead of charge. [default: %(default)s]''')


args = parser.parse_args()

if not args.xyzfile:
    exit('You must specify a .xyz input file.')

if args.xyzfile[-4:] == '.xyz':
    xyzfile = args.xyzfile[:-4]
else:
    xyzfile = args.xyzfile

if not os.path.isfile('{0}.xyz'.format(xyzfile)):
    exit('{0}.xyz not found.'.format(xyzfile))

if args.method == 'DFT':
    method = 'KSDFT\n' + args.xcfun
elif args.method == 'HF':
    method = ''
else:
    exit('Invalid QM method for force field calculation.')

if args.bonds:
    bmids = ''
else:
    bmids = 'BONDS\n0.0\n'

if args.occu:
    num_occu=get_occu(xyzfile+".xyz",args.charge)
    occu_line="OCCU\n"+str(num_occu)+"\n"
else:
    occu_line="CHARGE\n" + str(args.charge) + "\n"


input = ("&GATEWAY &END\n" +
         "TITLE\n" +
         "Generated by subloprop 0.1\n" +
         "COORD\n" +
         "{0}.xyz\n".format(xyzfile) +
         "BASIS\n" +
         "{0}\n".format(args.basis) +
         "GROUP\n" +
         "C1\n" +
         "END OF INPUT\n" +
         "&SEWARD &END\n" +
         "MULT\n" +
         "2\n" +
         "CHOL\n" +
         "END OF INPUT\n" +
         "&SCF &END\n" +
         "{0}".format(occu_line) +
#         'CHARGE\n' +
#         '{0}\n'.format(args.charge) +
         'CHOL\n' +
         '{0}\n'.format(method) +
         'END OF INPUT\n' +
         '&LOPROP &END\n' +
         'MPPROP\n' +
         '{0}\n'.format(args.multi) +
         '{0}'.format(bmids) +
         'END OF INPUT\n')
finp = open('{0}.inp'.format(xyzfile), 'w')
finp.write(input)
finp.close()