orca2dalton.py

#!/usr/bin/env python

'''
    Script to convert ORCA basis set file to Dalton library format.

    Reads:
        stdin

    Writes:
        stdout

    Example usage:
        $ ./orca2dalton.py < ano-pV5Z.bas > ano-pV5Z.dalton

    Author:
        Radovan Bast <lastname at kth.se>
'''

import sys
from chemistry import atomic_number, l_name_number, l_number_name


def zero_array(n):
    a = []
    for i in range(n):
        a.append(0.0)
    return a


def zero_matrix(k, l):
    m = []
    for i in range(k):
        m.append(zero_array(l))
    return m


def add_exp_to_list(l, e):
    l_new = l[:]
    e_found = False
    for x in l_new:
        if abs(x - e) < 1.0e-10:
            e_found = True
    if not e_found:
        l_new.append(e)
    return l_new


def find_element(l, e):
    for i, x in enumerate(l):
        if abs(x - e) < 1.0e-10:
            return i


def read_chunk(chunk):

    # first we get max L value
    l_max = 0
    for line in chunk:
        if len(line.split()) == 2:
            l_str = line.split()[0].lower()
            if l_str not in l_name_number.keys():
                sys.stderr.write('L quantum number "%s" not mapped\n' % l_str)
                sys.exit(1)
            l = l_name_number[l_str]
            if l > l_max: l_max = l

    exponents = []
    num_blocks = []
    for l in range(l_max + 1):
        exponents.append([])
        num_blocks.append(0)

    # get number of blocks and exponents
    for line in chunk:
        if len(line.split()) == 2:
            l_str = line.split()[0].lower()
            l = l_name_number[line.split()[0].lower()]
            num_blocks[l] += 1
        if len(line.split()) == 3:
            e = float(line.split()[1])
            exponents[l] = add_exp_to_list(exponents[l], e)

    # now we know all dimensions
    coefficients = []
    num_blocks_local = []
    for l in range(l_max + 1):
        coefficients.append(zero_matrix(len(exponents[l]), num_blocks[l]))
        num_blocks_local.append(0)

    # now we read in actual numbers
    for line in chunk:
        if len(line.split()) == 2:
            l_str = line.split()[0].lower()
            l = l_name_number[line.split()[0].lower()]
            num_blocks_local[l] += 1
        if len(line.split()) == 3:
            e = float(line.split()[1])
            c = float(line.split()[2])
            i = find_element(exponents[l], e)
            j = num_blocks_local[l] - 1
            coefficients[l][i][j] = c

    return exponents, coefficients


def main():

    lines = []
    # read lines from stdin
    for line in sys.stdin.read().splitlines():
        # remove trailing whitespace
        lines.append(line.rstrip())

    # title is the first line
    title = lines[0]

    atomic_number_lowercase = {}
    for atom in atomic_number.keys():
        atomic_number_lowercase[atom.lower()] = atomic_number[atom]

    # extract comments and atom text chunks
    comments = []
    atom_chunk = {}
    for line in lines[1:]: # ignore first line
        if len(line) > 0:
            if line[0] == '!':
                comments.append(line[1:])
            elif line == 'STOP':
                # we ignore this
                pass
            elif line == '':
                # we ignore this
                pass
            elif len(line.split()) == 1:
                if line.lower() in atomic_number_lowercase:
                    # we start a new chunk
                    current_atom = line.lower()
                    atom_chunk[current_atom] = []
                else:
                    sys.stderr.write('element "%s" not recognized - perhaps a typo?\n' % line)
                    sys.exit(1)
            else:
                # we append to the chunk
                atom_chunk[current_atom].append(line)

    # print header
    print('$ %s' % title)
    print('$')
    print('$ [automatically generated by orca2dalton.py]')
    print('$')
    if len(comments) > 0:
        print('$ comments copied from original file:')
        print('$ -----------------------------------')
        for line in comments:
            if line != '':
                print('$ %s' % line)
        print('$ -----------------------------------')
        print('$')

    # print all elements
    for atom in atomic_number:
        if atom.lower() in atom_chunk.keys():
            print('a %i' % atomic_number_lowercase[atom.lower()])
            print('$ %s' % atom)
            exponents, coefficients = read_chunk(atom_chunk[atom.lower()])
            for l in range(len(exponents)):
                print('$ %s-functions' % l_number_name[l])
                num_exponents = len(exponents[l])
                num_blocks = len(coefficients[l][0])
                print('%5i%5i%5i' % (num_exponents, num_blocks, 0))
                for i in range(len(exponents[l])):
                    s = '%15f' % exponents[l][i]
                    k = 0
                    for j in range(len(coefficients[l][i])):
                        if k == 6:
                            s += '\n   '
                            k = 0
                        s += '%12.8f' % coefficients[l][i][j]
                        k += 1
                    print(s)


if __name__ == '__main__':
    main()