xyz2mol.py

#!/usr/bin/env python
# Arnfinn Hykkerud Steindal, Tromso Oct. 2010
# Converts a xyz-file to a mol-file for Dalton
# TODO: 1) convert from au to angstrom, if needed

import re
import sys
import shutil
import string
import os

import argparse as ap
from datetime import datetime

name2number = {'H': 1, 'He': 2, 'Li': 3, 'Be': 4, 'B': 5, 'C': 6, 'N': 7,
               'O': 8, 'F': 9, 'Ne': 10, 'Na': 11, 'Mg': 12, 'Al': 13,
               'Si': 14, 'P': 15, 'S': 16, 'Cl': 17, 'Ar': 18, 'K': 19,
               'Ca': 20, 'Sc': 21, 'Ti': 22, 'V': 23, 'Cr': 24, 'Mn': 25,
               'Fe': 26, 'Co': 27, 'Ni': 28, 'Cu': 29, 'Zn': 30, 'Ga': 31,
               'Ge': 32, 'As': 33, 'Se': 34, 'Br': 35, 'Kr': 36, 'Rb': 37,
               'Sr': 38, 'Y': 39, 'Zr': 40, 'Nb': 41, 'Mo': 42, 'Tc': 43,
               'Ru': 44, 'Rh': 45, 'Pd': 46, 'Ag': 47, 'Cd': 48, 'In': 49,
               'Sn': 50, 'Sb': 51, 'Te': 52, 'I': 53, 'Xe': 54, 'Cs': 55,
               'Ba': 56, 'La': 57, 'Ce': 58, 'Pr': 59, 'Nd': 60, 'Pm': 61,
               'Sm': 62, 'Eu': 63, 'Gd': 64, 'Tb': 65, 'Dy': 66, 'Ho': 67,
               'Er': 68, 'Tm': 69, 'Yb': 70, 'Lu': 71, 'Hf': 72, 'Ta': 73,
               'W': 74, 'Re': 75, 'Os': 76, 'Ir': 77, 'Pt': 78, 'Au': 79,
               'Hg': 80, 'Tl': 81, 'Pb': 82, 'Bi': 83, 'Po': 84, 'At': 85,
               'Rn': 86, 'Fr': 87, 'Ra': 88}

def xyz2mol(lines):
    all_atoms = []
    try:
        num_atoms = int(lines[0].split()[0])
    except:
        print "header: {0}".format(lines[0].split()[0])
        quit("Not correct format in the xyz-file header")
    for i in lines[2:]:
        words = i.split()
        atm_name = words[0]
        if len(atm_name) > 1:
            try:
                int(atm_name[1])
                atom = atm_name[0]
            except:
                atom = atm_name[0:2]
        else:
            atom = atm_name
        try:
            atom_charge = name2number[atom]
        except:
            quit("Element {0} not found in name2number dictionary.".format(atom))
        xcoord = words[1]
        ycoord = words[2]
        zcoord = words[3]
        all_atoms.append((atom_charge,atm_name,xcoord,ycoord,zcoord))
    sort_all = sorted(all_atoms, key=lambda atoms: atoms[0], reverse=True)
    return sort_all

def create_mollines(all_atoms,charge, basis):
    atomtypes = 0
    old_charge = 999
    tmp_mol = ""
    for i in all_atoms:
        if i[0] != old_charge:
            if old_charge != 999:
                tmp_mol += "     Charge={0}   Atoms={1}\n".format(old_charge,k)
                tmp_mol += group
            k = 0
            atomtypes += 1
            group = ""
            old_charge = i[0]
        lenx = max(12-len(i[2]),2)
        leny = max(12-len(i[3]),2)
        lenz = max(12-len(i[4]),2)
        group += i[1] + " "*lenx + i[2] + " "*leny + i[3] + " "*lenz + i[4] + "\n"
        k += 1
    tmp_mol += "     Charge={0}   Atoms={1}\n".format(old_charge,k)
    tmp_mol += group
#    atomtypes += 1
    mol_lines  = "BASIS\n{0}\n".format(basis)
    mol_lines += "Generated by Arnfinns xyz2mol.py\n"
    mol_lines += "date: {0}\n".format(datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
    mol_lines += "AtomTypes={0} Charge={1} Angstrom NoSymmetry\n".format(atomtypes,charge)
    mol_lines += tmp_mol
    return mol_lines


parser = ap.ArgumentParser()
parser.add_argument("-i", "--input",dest="filename", help="the file to read")
parser.add_argument("-b", "--basis",dest="basis", default='aug-cc-pVDZ', help='''the basisset. [default: %(default)s]''')
parser.add_argument("-q", "--charge",dest="charge", default = "0", help='''the molecular charge''')

args = parser.parse_args()

filename=args.filename
try:
    charge = int(args.charge)
except:
    quit("ERROR: {0} not a valid charge!".format(args.charge))

# read xyz file
xyz = open(filename,'r')
lines = xyz.readlines()
xyz.close()

# extract data from file
all_atoms = xyz2mol(lines)

# total number of protons and electrons
tot_protons = 0
for i in all_atoms:
    tot_protons += i[0]
tot_electrons = tot_protons - charge
print "Total number of protons: {0}".format(tot_protons)
print "Total number of electrons: {0}".format(tot_electrons)
if tot_electrons%2 != 0:
    tmp = ""
    if charge > 0:
        tmp = "+"
    print "WARNING!"
    print "Odd number of electrons: {0}! Correct total charge ({1}{2})?".format(tot_electrons,tmp,charge)

# Write mol file
a=len(filename)
mol=open(filename[0:a-4]+'.mol','w')
mol.write(create_mollines(all_atoms,charge,args.basis))
mol.close()