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