1import math, sys 2 3# Element info taken mostly from Rasmol but the effective covalent 4# radii are from NWChem with those probably coming from Hondo. 5 6cpk = {0: "light grey", 7 1: "sky blue", 8 2: "red", 9 3: "yellow", 10 4: "white", 11 5: "pink", 12 6: "gold", 13 7: "blue", 14 8: "orange", 15 9: "dark grey", 16 10: "brown", 17 11: "purple", 18 12: "deep pink", 19 13: "green", 20 14: "fire brick", 21 15: "mid green"} 22 23# Colors from X11/rgb.txt 24colors = { 25 "light grey": (211,211,211), 26 "sky blue": (135,206,235), 27 "red": (255,0,0), 28 "yellow": (255,255,0), 29 "white": (255,255,255), 30 "pink": (255,105,180), 31 "gold": (255,215,0), 32 "blue": (0,0,255), 33 "orange": (255,165,0), 34 "dark grey": (169,169,169), 35 "brown": (165,42,42), 36 "purple": (160,32,240), 37 "deep pink": (205,16,118), 38 "green": (0,255,0), 39 "fire brick": (178,34,34), 40 "mid green": (60,179,113)} 41 42for color in colors.keys(): 43 r,g,b = colors[color] 44 colors[color]=r/255.0,g/255.0,b/255.0 45 46class Element: 47 def __init__(self,symbol,covrad,vdwrad,col,name,atno): 48 self.symbol = symbol 49 self.covrad = covrad 50 self.vdwrad = vdwrad 51 self.col = colors[cpk[col]] 52 self.name = name 53 self.atno = atno 54 55 def __str__(self): 56 return "%3d %2s %16s %.1f %.1f (%.2f %.2f %.2f) " % ( 57 self.atno,self.symbol,self.name,self.covrad,self.vdwrad, 58 self.col[0],self.col[1],self.col[2]) 59 60elements = range(104) 61elements[ 0] = Element('Bq' ,0.00, 0.0, 0, "Dummy" , 0 ) 62elements[ 1] = Element('H' ,0.30, 2.75, 4, "HYDROGEN" , 1 ) 63elements[ 2] = Element('He',1.22, 5.50, 5, "HELIUM" , 2 ) 64elements[ 3] = Element('Li',1.23, 3.05, 14, "LITHIUM" , 3 ) 65elements[ 4] = Element('Be',0.89, 1.57, 12, "BERYLLIUM" , 4 ) 66elements[ 5] = Element('B' ,0.88, 3.87, 13, "BORON" , 5 ) 67elements[ 6] = Element('C' ,0.77, 3.87, 9, "CARBON" , 6 ) # was color 0 68elements[ 7] = Element('N' ,0.70, 3.50, 1, "NITROGEN" , 7 ) 69elements[ 8] = Element('O' ,0.66, 3.37, 2, "OXYGEN" , 8 ) 70elements[ 9] = Element('F' ,0.58, 3.25, 6, "FLUORINE" , 9 ) 71elements[ 10] = Element('Ne',1.60, 5.05, 12, "NEON" , 10 ) 72elements[ 11] = Element('Na',1.66, 5.50, 7, "SODIUM" , 11 ) 73elements[ 12] = Element('Mg',1.36, 3.75, 15, "MAGNESIUM" , 12 ) 74elements[ 13] = Element('Al',1.25, 3.75, 9, "ALUMINIUM" , 13 ) 75elements[ 14] = Element('Si',1.17, 5.50, 6, "SILICON" , 14 ) 76elements[ 15] = Element('P' ,1.10, 4.70, 8, "PHOSPHORUS" , 15 ) 77elements[ 16] = Element('S' ,1.04, 4.52, 3, "SULPHUR" , 16 ) 78elements[ 17] = Element('Cl',0.99, 4.37, 13, "CHLORINE" , 17 ) 79elements[ 18] = Element('Ar',1.91, 6.92, 12, "ARGON" , 18 ) 80elements[ 19] = Element('K' ,2.03, 5.97, 12, "POTASSIUM" , 19 ) 81elements[ 20] = Element('Ca',1.74, 4.87, 9, "CALCIUM" , 20 ) 82elements[ 21] = Element('Sc',1.44, 3.30, 12, "SCANDIUM" , 21 ) 83elements[ 22] = Element('Ti',1.32, 4.87, 9, "TITANIUM" , 22 ) 84elements[ 23] = Element('V' ,1.22, 2.65, 12, "VANADIUM" , 23 ) 85elements[ 24] = Element('Cr',1.19, 2.82, 9, "CHROMIUM" , 24 ) 86elements[ 25] = Element('Mn',1.17, 2.97, 9, "MANGANESE" , 25 ) 87elements[ 26] = Element('Fe',1.165,4.87, 8, "IRON" , 26 ) 88elements[ 27] = Element('Co',1.16, 2.82, 12, "COBALT" , 27 ) 89elements[ 28] = Element('Ni',1.15, 3.10, 10, "NICKEL" , 28 ) 90elements[ 29] = Element('Cu',1.17, 2.87, 10, "COPPER" , 29 ) 91elements[ 30] = Element('Zn',1.25, 2.87, 10, "ZINC" , 30 ) 92elements[ 31] = Element('Ga',1.25, 3.87, 12, "GALLIUM" , 31 ) 93elements[ 32] = Element('Ge',1.22, 9.99, 12, "GERMANIUM" , 32 ) 94elements[ 33] = Element('As',1.21, 2.07, 12, "ARSENIC" , 33 ) 95elements[ 34] = Element('Se',1.17, 2.25, 12, "SELENIUM" , 34 ) 96elements[ 35] = Element('Br',1.14, 4.37, 10, "BROMINE" , 35 ) 97elements[ 36] = Element('Kr',1.98, 4.75, 12, "KRYPTON" , 36 ) 98elements[ 37] = Element('Rb',2.22, 6.62, 12, "RUBIDIUM" , 37 ) 99elements[ 38] = Element('Sr',1.92, 5.05, 12, "STRONTIUM" , 38 ) 100elements[ 39] = Element('Y' ,1.62, 4.02, 12, "YTTRIUM" , 39 ) 101elements[ 40] = Element('Zr',1.45, 3.55, 12, "ZIRCONIUM" , 40 ) 102elements[ 41] = Element('Nb',1.34, 3.32, 12, "NIOBIUM" , 41 ) 103elements[ 42] = Element('Mo',1.29, 4.37, 12, "MOLYBDENUM" , 42 ) 104elements[ 43] = Element('Tc',1.27, 4.50, 12, "TECHNETIUM" , 43 ) 105elements[ 44] = Element('Ru',1.24, 3.00, 12, "RUTHENIUM" , 44 ) 106elements[ 45] = Element('Rh',1.25, 3.05, 12, "RHODIUM" , 45 ) 107elements[ 46] = Element('Pd',1.28, 3.60, 12, "PALLADIUM" , 46 ) 108elements[ 47] = Element('Ag',1.34, 3.87, 9, "SILVER" , 47 ) 109elements[ 48] = Element('Cd',1.41, 4.37, 12, "CADMIUM" , 48 ) 110elements[ 49] = Element('In',1.50, 3.62, 12, "INDIUM" , 49 ) 111elements[ 50] = Element('Sn',1.40, 4.17, 12, "TIN" , 50 ) 112elements[ 51] = Element('Sb',1.41, 2.80, 12, "ANTIMONY" , 51 ) 113elements[ 52] = Element('Te',1.37, 3.15, 12, "TELLURIUM" , 52 ) 114elements[ 53] = Element('I' ,1.33, 4.37, 11, "IODINE" , 53 ) 115elements[ 54] = Element('Xe',2.09, 5.25, 12, "XENON" , 54 ) 116elements[ 55] = Element('Cs',2.35, 7.52, 12, "CAESIUM" , 55 ) 117elements[ 56] = Element('Ba',1.98, 6.02, 8, "BARIUM" , 56 ) 118elements[ 57] = Element('La',1.69, 4.57, 12, "LANTHANUM" , 57 ) 119elements[ 58] = Element('Ce',1.65, 4.65, 12, "CERIUM" , 58 ) 120elements[ 59] = Element('Pr',1.65, 4.05, 12, "PRASEODYMIUM", 59 ) 121elements[ 60] = Element('Nd',1.64, 4.47, 12, "NEODYMIUM" , 60 ) 122elements[ 61] = Element('Pm',1.65, 4.40, 12, "PROMETHIUM" , 61 ) 123elements[ 62] = Element('Sm',1.66, 4.35, 12, "SAMARIUM" , 62 ) 124elements[ 63] = Element('Eu',1.65, 4.90, 12, "EUROPIUM" , 63 ) 125elements[ 64] = Element('Gd',1.61, 4.22, 12, "GADOLINIUM" , 64 ) 126elements[ 65] = Element('Tb',1.59, 4.15, 12, "TERBIUM" , 65 ) 127elements[ 66] = Element('Dy',1.59, 4.07, 12, "DYSPROSIUM" , 66 ) 128elements[ 67] = Element('Ho',1.58, 4.02, 12, "HOLMIUM" , 67 ) 129elements[ 68] = Element('Er',1.57, 3.97, 12, "ERBIUM" , 68 ) 130elements[ 69] = Element('Tm',1.56, 3.92, 12, "THULIUM" , 69 ) 131elements[ 70] = Element('Yb',1.56, 3.85, 12, "YTTERBIUM" , 70 ) 132elements[ 71] = Element('Lu',1.56, 3.82, 12, "LUTETIUM" , 71 ) 133elements[ 72] = Element('Hf',1.44, 3.50, 12, "HAFNIUM" , 72 ) 134elements[ 73] = Element('Ta',1.34, 3.05, 12, "TANTALUM" , 73 ) 135elements[ 74] = Element('W' ,1.30, 3.15, 12, "TUNGSTEN" , 74 ) 136elements[ 75] = Element('Re',1.28, 3.25, 12, "RHENIUM" , 75 ) 137elements[ 76] = Element('Os',1.26, 3.95, 12, "OSMIUM" , 76 ) 138elements[ 77] = Element('Ir',1.26, 3.05, 12, "IRIDIUM" , 77 ) 139elements[ 78] = Element('Pt',1.29, 3.87, 12, "PLATINUM" , 78 ) 140elements[ 79] = Element('Au',1.34, 3.62, 6, "GOLD" , 79 ) 141elements[ 80] = Element('Hg',1.44, 4.95, 12, "MERCURY" , 80 ) 142elements[ 81] = Element('Tl',1.55, 4.27, 12, "THALLIUM" , 81 ) 143elements[ 82] = Element('Pb',1.54, 5.40, 12, "LEAD" , 82 ) 144elements[ 83] = Element('Bi',1.52, 4.32, 12, "BISMUTH" , 83 ) 145elements[ 84] = Element('Po',1.53, 3.02, 12, "POLONIUM" , 84 ) 146elements[ 85] = Element('At',1.50, 2.80, 12, "ASTATINE" , 85 ) 147elements[ 86] = Element('Rn',2.20, 5.75, 12, "RADON" , 86 ) 148elements[ 87] = Element('Fr',3.24, 8.10, 12, "FRANCIUM" , 87 ) 149elements[ 88] = Element('Ra',2.68, 6.42, 12, "RADIUM" , 88 ) 150elements[ 89] = Element('Ac',2.25, 5.30, 12, "ACTINIUM" , 89 ) 151elements[ 90] = Element('Th',2.16, 4.60, 12, "THORIUM" , 90 ) 152elements[ 91] = Element('Pa',1.93, 4.00, 12, "PROTACTINIUM", 91 ) 153elements[ 92] = Element('U' ,3.00, 4.37, 12, "URANIUM" , 92 ) 154elements[ 93] = Element('Np',1.57, 4.27, 12, "NEPTUNIUM" , 93 ) 155elements[ 94] = Element('Pu',1.81, 4.17, 12, "PLUTONIUM" , 94 ) 156elements[ 95] = Element('Am',2.21, 4.15, 12, "AMERICIUM" , 95 ) 157elements[ 96] = Element('Cm',1.43, 4.12, 12, "CURIUM" , 96 ) 158elements[ 97] = Element('Bk',1.42, 4.10, 12, "BERKELIUM" , 97 ) 159elements[ 98] = Element('Cf',1.40, 4.07, 12, "CALIFORNIUM" , 98 ) 160elements[ 99] = Element('Es',1.39, 4.05, 12, "EINSTEINIUM" , 99 ) 161elements[100] = Element('Fm',1.38, 4.02, 12, "FERMIUM" ,100 ) 162elements[101] = Element('Md',1.37, 4.00, 12, "MENDELEVIUM" ,101 ) 163elements[102] = Element('No',1.36, 3.97, 12, "NOBELIUM" ,102 ) 164elements[103] = Element('Lr',1.34, 3.95, 12, "LAWRENCIUM" ,103 ) 165 166default_radii = range(104) 167default_colors = range(104) 168for i in range(1,len(elements)): 169 default_radii[i] = elements[i].covrad 170 default_colors[i] = elements[i].col 171 172def scale_coords(coords, fac): 173 natom = len(coords) 174 for atom in range(natom): 175 for i in range(3): 176 coords[atom][i] = coords[atom][i]*fac 177 178def dist(i,j): 179 return math.sqrt((i[0]-j[0])**2+ 180 (i[1]-j[1])**2+ 181 (i[2]-j[2])**2) 182 183def make_bonds(coords, z, radii): 184 bonds = [] 185 for i in range(len(coords)): 186 for j in range(i): 187 rij = dist(coords[i],coords[j]) 188 if rij < 1.1*(radii[z[i]]+radii[z[j]]): 189 #print "bond",i,j,rij 190 bonds.append([i,j]) 191 return bonds 192 193def moldx(coords,z,bonds=None,units="angstrom",fname="mol.dx", 194 radii=default_radii,colors=default_colors, 195 scale=None,shift=None): 196 ''' 197 198 Make an Opendx format file with a stick representation of the 199 molecule. Additional representations will be easy to add. 200 201 coords = coords[natom][3] input coords (unchanged on output) 202 203 z = z[natom] input atomic number of each atom (unchanged on output) 204 205 bonds = optional list of precomputed connections. If None, these 206 . are computed from the covalent radii. 207 208 units = the input units. If the input units are atomic, they will be 209 . converted to angstrom. The output units are always angstrom 210 . unless additional scaling is specified. 211 212 fname = the name of the output file 213 214 radii = covalent radii for forming bonds. Defaults will be used 215 . if not speficied. 216 217 colors = colors for the atoms (3-vector float in [0,1]). Defaults 218 . will be used if not specified. 219 220 scale = if scale (scalar float) is specified, the *output* 221 . coordinates are scaled by that factor. 222 223 shift = if shift (3-vector float) is specified, the *output* 224 . coordinates are shifted *after* scaling 225 . (e.g., x*scale+shift[0]) 226 227 ''' 228 229 au_to_angs = 0.529177249 230 if len(z) != len(coords): raise ValueError 231 232 if not bonds: 233 # Need angstroms for making bonds 234 if units in ["au","a.u.","atomic","bohr"]: 235 scale_coords(coords,au_to_angs) 236 237 bonds = make_bonds(coords,z,radii) 238 239 if units in ["au","a.u.","atomic","bohr"]: 240 scale_coords(coords,1.0/au_to_angs) 241 242 # We have the atomic coordinates, numbers, radii, and 243 # connectivity. The opendx stick representation is formed 244 # by splitting the bonds at the midpoint (weighted by 245 # the ratio of the radii) with appropriate coloring. 246 247 if not units in ["au","a.u.","atomic","bohr"]: 248 # Need atomic units for the output 249 scale_coords(coords,1.0/au_to_angs) 250 251 pos = [] 252 col = [] 253 con = [] 254 255 for i,j in bonds: 256 npos = len(pos) 257 rij = dist(coords[i],coords[j]) 258 si = radii[z[i]]/(radii[z[i]]+radii[z[j]]) 259 sj = 1.0-si 260 mid = [coords[i][0]*sj+coords[j][0]*si, 261 coords[i][1]*sj+coords[j][1]*si, 262 coords[i][2]*sj+coords[j][2]*si] 263 264 pos.append(coords[i]); col.append(colors[z[i]]) 265 pos.append(mid); col.append(colors[z[i]]) 266 267 pos.append(mid); col.append(colors[z[j]]) 268 pos.append(coords[j]); col.append(colors[z[j]]) 269 270 con.append([npos,npos+1]) 271 con.append([npos+2,npos+3]) 272 273 f = open(fname,'w+') 274 275 f.write("object 1 class array type float rank 1 shape 3 items %d data follows\n" % len(pos)) 276 for x,y,z in pos: 277 if scale: x, y, z = x*scale, y*scale, z*scale 278 if shift: x, y, z = x+shift[0], y+shift[1], z+shift[2] 279 f.write("%f %f %f\n" % (x,y,z)) 280 f.write('attribute "dep" string "positions"\n') 281 282 f.write("object 2 class array type int rank 1 shape 2 items %d data follows\n" % len(con)) 283 for i,j in con: f.write("%d %d \n" % (i,j)) 284 f.write('attribute "element type" string "lines"\n') 285 f.write('attribute "ref" string "positions"\n') 286 f.write('attribute "dep" string "connections"\n') 287 288 f.write("object 3 class array type float rank 0 items %d data follows\n" % len(pos)) 289 for p in pos: f.write("0 ") 290 f.write("\n") 291 f.write('attribute "dep" string "positions"\n') 292 293 f.write("object 4 class array type float rank 1 shape 3 items %d data follows\n" % len(col)) 294 for r,g,b in col: f.write("%f %f %f\n" % (r,g,b)) 295 f.write('attribute "dep" string "positions"\n') 296 297 f.write('object "irregular positions irregular connections" class field\n') 298 f.write('component "positions" value 1\n') 299 f.write('component "connections" value 2\n') 300 f.write('component "data" value 3\n') 301 f.write('component "colors" value 4\n') 302 f.write('attribute "name" string "irregular positions irregular connections"\n') 303 f.write('end\n') 304 305 f.close() 306 307 # Undo-Force a.u. while making the bonds (cannot do until after 308 # writing file since pos contains references to these values). 309 if not units in ["au","a.u.","atomic","bohr"]: 310 scale_coords(coords,au_to_angs) 311 312 313def symtoatn(sym): 314 sym = sym.lower() 315 for el in elements: 316 if el.symbol.lower() == sym: 317 return el.atno 318 raise "Uh?" 319 320if __name__ == "__main__": 321 # Read first molecule from the madness output file 322 f = sys.stdin 323 while 1: 324 line = f.readline().split() 325 if len(line) and line[0] == 'geometry': 326 coords = [] 327 z = [] 328 units = 'au' 329 while 1: 330 line = f.readline().split() 331 if len(line) == 0: continue 332 if line[0] == 'end': 333 break 334 elif line[0] == 'units': 335 units = line[1] 336 elif not (line[0] == 'units' or line[0] == 'eprec'): 337 z.append(symtoatn(line[0])) 338 coords.append([float(line[1]), float(line[2]), float(line[3])]) 339 moldx(coords, z, units=units, fname='molecule.dx') 340 sys.exit(0) 341 342 343