1#!/usr/bin/env python
2#
3#   Polarizable Embedding (PE) library
4#   Copyright (C) 2013, 2014 The PE library developers. See the CONTRIBUTORS file
5#                      in the top-level directory of this distribution.
6#
7#   This tool is part of the PE library.
8#
9#   The PE library is free software: you can redistribute it and/or modify
10#   it under the terms of the GNU Lesser General Public License as
11#   published by the Free Software Foundation, either version 3 of the
12#   License, or (at your option) any later version.
13#
14#   The PE library is distributed in the hope that it will be useful,
15#   but WITHOUT ANY WARRANTY; without even the implied warranty of
16#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17#   GNU Lesser General Public License for more details.
18#
19#   You should have received a copy of the GNU Lesser General Public License
20#   along with the PE library. If not, see <http://www.gnu.org/licenses/>.
21#
22#   Contact information:
23#
24#   Jogvan Magnus Haugaard Olsen
25#   E-mail: foeroyingur@gmail.com
26#
27
28import os
29import sys
30import argparse as ap
31
32
33parser = ap.ArgumentParser(description='Potential Converter: convert potential files to new format',
34                           usage='%(prog)s [options]',
35                           fromfile_prefix_chars='@')
36
37parser.add_argument('--version', action='version', version='Potential Converter 1.0')
38
39parser.add_argument('-i', dest='inputfile', metavar='INPUT_FILE',
40                    help='''Specify the name of the input file.''')
41
42parser.add_argument('-o', dest='outputfile', metavar='OUTPUT_FILE',
43                    help='''Specify the name of the output file.''')
44
45parser.add_argument('--input-type', dest='inputtype', metavar='TYPE',
46                    choices=['QMMM'], default='QMMM',
47                    help='''Specify the format of the input file. Valid choices
48                            are %(choices)s. [default: %(default)s]''')
49
50args = parser.parse_args()
51
52pyver = sys.version_info
53if pyver[0] < 2 or pyver[0] == 2 and pyver[1] < 7:
54    exit('ERROR: Python >= 2.7 required.')
55
56if not args.inputfile:
57    exit('ERROR: no input file specified.')
58elif not os.path.isfile(args.inputfile):
59    exit('ERROR: input file not found.')
60elif not args.outputfile:
61    exit('ERROR: no outputfile specified.')
62elif os.path.isfile(args.outputfile):
63    exit('ERROR: output file already exists.')
64
65inputdir, inputfile = os.path.split(args.inputfile)
66outputdir, outputfile = os.path.split(args.outputfile)
67
68charge2elem = { 0:  'X',  1:  'H',  2: 'He',  3: 'Li',  4: 'Be',  5:  'B',
69                6:  'C',  7:  'N',  8:  'O',  9:  'F', 10: 'Ne', 11: 'Na',
70               12: 'Mg', 13: 'Al', 14: 'Si', 15:  'P', 16:  'S', 17: 'Cl',
71               18: 'Ar', 19:  'K', 20: 'Ca', 21: 'Sc', 22: 'Ti', 23:  'V',
72               24: 'Cr', 25: 'Mn', 26: 'Fe', 27: 'Co', 28: 'Ni', 29: 'Cu',
73               30: 'Zn', 31: 'Ga', 32: 'Ge', 33: 'As', 34: 'Se', 35: 'Br',
74               36: 'Kr', 37: 'Rb', 38: 'Sr', 39:  'Y', 40: 'Zr', 41: 'Nb',
75               42: 'Mo', 43: 'Tc', 44: 'Ru', 45: 'Rh', 46: 'Pd', 47: 'Ag',
76               48: 'Cd', 49: 'In', 50: 'Sn', 51: 'Sb', 52: 'Te', 53:  'I',
77               54: 'Xe', 55: 'Cs', 56: 'Ba', 57: 'La', 58: 'Ce', 59: 'Pr',
78               60: 'Nd', 61: 'Pm', 62: 'Sm', 63: 'Eu', 64: 'Gd', 65: 'Tb',
79               66: 'Dy', 67: 'Ho', 68: 'Er', 69: 'Tm', 70: 'Yb', 71: 'Lu',
80               72: 'Hf', 73: 'Ta', 74:  'W', 75: 'Re', 76: 'Os', 77: 'Ir',
81               78: 'Pt', 79: 'Au', 80: 'Hg', 81: 'Tl', 82: 'Pb', 83: 'Bi',
82               84: 'Po', 85: 'At', 86: 'Rn', 87: 'Fr', 88: 'Ra'}
83
84au2aa = 0.5291772108
85
86if args.inputtype == 'QMMM':
87    fin = open('{}'.format(args.inputfile), 'r')
88    line = fin.readline().split()
89    unit = line[0]
90    line = fin.readline().split()
91    nsites = int(line[0])
92    mulorder = int(line[1])
93    polorder = int(line[2])
94    lexlst = int(line[3])
95    if len(line) == 5:
96        lelems = 1
97    elif len(line) == 4:
98        lelems = 0
99    elems = []
100    coords = []
101    M0s = []
102    M1s = []
103    M2s = []
104    M3s = []
105    M4s = []
106    M5s = []
107    alphas = []
108    exclists = []
109    if mulorder == 5:
110        pad = 56
111    elif mulorder == 4:
112        pad = 35
113    elif mulorder == 3:
114        pad = 20
115    elif mulorder == 2:
116        pad = 10
117    elif mulorder == 1:
118        pad = 4
119    elif mulorder == 0:
120        pad = 1
121    for line in fin:
122        line = line.split()
123        if not line:
124            continue
125        exclists.append(map(int, line[0:lexlst]))
126        coords.append(map(float, line[lexlst+lelems:lexlst+lelems+3]))
127        if lelems == 1:
128            elems.append(int(line[lexlst]))
129        elif lelems == 0:
130            elems.append(0)
131        if mulorder >= 0:
132            M0s.append(float(line[lexlst+lelems+3]))
133        if mulorder >= 1:
134            M1s.append(map(float, line[lexlst+lelems+4:lexlst+lelems+7]))
135        if mulorder >= 2:
136            M2s.append(map(float, line[lexlst+lelems+7:lexlst+lelems+13]))
137        if mulorder >= 3:
138            M3s.append(map(float, line[lexlst+lelems+13:lexlst+lelems+23]))
139        if mulorder >= 4:
140            M4s.append(map(float, line[lexlst+lelems+23:lexlst+lelems+38]))
141        if mulorder >= 5:
142            M5s.append(map(float, line[lexlst+lelems+38:lexlst+lelems+59]))
143        if polorder == 1:
144            isoalpha = float(line[lexlst+lelems+3+pad])
145            alphas.append([isoalpha, 0.0, 0.0, isoalpha, 0.0, isoalpha])
146        elif polorder == 2:
147            start = lexlst+lelems+pad+3
148            end = lexlst+lelems+pad+9
149            alphas.append(map(float, line[start:end]))
150    fin.close()
151
152    if exclists[0][0] == exclists[1][0]:
153        exit('ERROR: conversion does not work for this kind of exclusion list')
154
155    for i, excl in enumerate(exclists):
156        if excl[0] != i+1:
157            exit('ERROR: conversion does not work for this kind of exclusion list')
158
159    if unit == 'AU':
160        coords = [[elem * au2aa for elem in coord] for coord in coords]
161    ndec = len(str(nsites))
162    body = '@COORDINATES\n'
163    body += '{}\n'.format(nsites)
164    body += 'AA\n'
165    for i, coord in enumerate(coords):
166        body += '{0} {1[0]:14.8f} {1[1]:14.8f} {1[2]:14.8f}\n'.format(charge2elem[elems[i]], coord)
167    if mulorder >= 0:
168        body += '@MULTIPOLES\n'
169    if mulorder >= 0:
170        body += 'ORDER 0\n'
171        body += '{0}\n'.format(nsites)
172        for i, M0 in enumerate(M0s):
173            body += '{0:{1}d} {2:14.8f}\n'.format(i+1, ndec, M0)
174    if mulorder >= 1:
175        body += 'ORDER 1\n'
176        body += '{0}\n'.format(nsites)
177        for i, M1 in enumerate(M1s):
178            body += '{0:{1}d}'.format(i+1, ndec)
179            for j in range(3):
180                body += ' {0:14.8f}'.format(M1[j])
181            body += '\n'
182    if mulorder >= 2:
183        body += 'ORDER 2\n'
184        body += '{0}\n'.format(nsites)
185        for i, M2 in enumerate(M2s):
186            body += '{0:{1}d}'.format(i+1, ndec)
187            for j in range(6):
188                body += ' {0:14.8f}'.format(M2[j])
189            body += '\n'
190    if mulorder >= 3:
191        body += 'ORDER 3\n'
192        body += '{0}\n'.format(nsites)
193        for i, M3 in enumerate(M3s):
194            body += '{0:{1}d}'.format(i+1, ndec)
195            for j in range(10):
196                body += ' {0:14.8f}'.format(M3[j])
197            body += '\n'
198    if mulorder >= 4:
199        body += 'ORDER 4\n'
200        body += '{0}\n'.format(nsites)
201        for i, M4 in enumerate(M4s):
202            body += '{0:{1}d}'.format(i+1, ndec)
203            for j in range(15):
204                body += ' {0:14.8f}'.format(M4[j])
205            body += '\n'
206    if mulorder >= 5:
207        body += 'ORDER 5\n'
208        body += '{0}\n'.format(nsites)
209        for i, M5 in enumerate(M5s):
210            body += '{0:{1}d}'.format(i+1, ndec)
211            for j in range(21):
212                body += ' {0:14.8f}'.format(M5[j])
213            body += '\n'
214    if polorder >= 1:
215        body += '@POLARIZABILITIES\n'
216        body += 'ORDER 1 1\n'
217        body += '{0}\n'.format(nsites)
218        for i, alpha in enumerate(alphas):
219            body += '{0:{1}d}'.format(i+1, ndec)
220            for j in range(6):
221                body += ' {0:14.8f}'.format(alpha[j])
222            body += '\n'
223        if exclists:
224            body += 'EXCLISTS\n'
225            body += '{0} {1}\n'.format(nsites, lexlst)
226            for i, exclist in enumerate(exclists):
227                for ex in exclist:
228                    body += ' {0:{1}}'.format(ex, ndec)
229                body += '\n'
230        body += '\n'
231    fout = open('{}'.format(args.outputfile), 'w')
232    fout.write(body)
233    fout.close()
234
235else:
236    exit('ERROR: potential is not supported')
237    fin = open('{}'.format(args.inputfile), 'r')
238    oldpot = fin.read()
239    fin.close()
240