1# 2# 3# lmpsdata.py 4# 5# For reading, writing and manipulating lammps data files 6# For calculation of certain properties using lammps data files 7# For creating VMD input text files using lammps data files 8# All x,y,z calculations assume the information includes image flags 9 10class Lmpsdata: 11 def __init__(self,file,atomtype): 12 """initiates lammps data structures""" 13 self.atomtype=atomtype 14 self.keywords=[] 15 self.atoms=[] 16 self.angles=[] 17 self.bonds=[] 18 self.dihedrals=[] 19 self.dipoles=[] 20 self.impropers=[] 21 self.masses=[] 22 self.shapes=[] 23 self.velocities=[] 24 self.anglecoef=[] 25 self.bondcoef=[] 26 self.dihedralcoef=[] 27 self.impropercoef=[] 28 self.paircoef=[] 29 self.read(file) 30 31 def read(self,file): 32 """Reads in lammps data file 33 Skips the first line of the file (Comment line) 34 First reads the header portion of the file (sectflag=1) 35 blank lines are skipped 36 header keywords delineate an assignment 37 if no header keyword is found, body portion begins 38 Second reads the body portion of the file (sectflag=2) 39 first line of a section has only a keyword 40 next line is skipped 41 remaining lines contain values 42 a blank line signifies the end of that section 43 if no value is listed on a line than a keyword must be used 44 File is read until the end 45 The keywords read in are stored in self.keywords""" 46 if file=='': 47 print 'no file is given. Will have to build keywords and class structures manually' 48 return 49 sectflag=1 50 f=open(file,'r') 51 f.readline() 52 for line in f: 53 row=line.split() 54 if sectflag==1: 55 if len(row)==0: 56 #skip line the line is blank 57 pass 58 elif len(row)==1: 59 # Set Sectflag to 2, assume line is a keyword 60 sectflag=2 61 checkkey=1 #ensures keyword will be checked in body portion 62 keyword=row 63 elif len(row)==2: 64 if row[1]=='atoms': 65 self.atomnum=row[0] 66 self.keywords.append('atoms') 67 elif row[1]=='bonds': 68 self.bondnum=row[0] 69 self.keywords.append('bonds') 70 elif row[1]=='angles': 71 self.anglenum=row[0] 72 self.keywords.append('angles') 73 elif row[1]=='dihedrals': 74 self.dihedralnum=row[0] 75 self.keywords.append('dihedrals') 76 elif row[1]=='impropers': 77 self.impropernum=row[0] 78 self.keywords.append('impropers') 79 else: 80 # Set Sectflag to 2, assume line is a keyword 81 sectflag=2 82 checkkey=1 #ensures keyword will be checked in body portion 83 keyword=row 84 elif len(row)==3: 85 if row[1]=='atom' and row[2]=='types': 86 self.atomtypenum=row[0] 87 self.keywords.append('atom types') 88 elif row[1]=='bond' and row[2]=='types': 89 self.bondtypenum=row[0] 90 self.keywords.append('bond types') 91 elif row[1]=='angle' and row[2]=='types': 92 self.angletypenum=row[0] 93 self.keywords.append('angle types') 94 elif row[1]=='dihedral' and row[2]=='types': 95 self.dihedraltypenum=row[0] 96 self.keywords.append('dihedral types') 97 elif row[1]=='improper' and row[2]=='types': 98 self.impropertypenum=row[0] 99 self.keywords.append('improper types') 100 else: 101 # Set Sectflag to 2, assume line is a keyword 102 sectflag=2 103 checkkey=1 #ensures keyword will be checked in body portion 104 keyword=row 105 elif len(row)==4: 106 if row[2]=='xlo' and row[3]=='xhi': 107 self.xdimension=[row[0], row[1]] 108 self.keywords.append('xlo xhi') 109 elif row[2]=='ylo' and row[3]=='yhi': 110 self.ydimension=[row[0], row[1]] 111 self.keywords.append('ylo yhi') 112 elif row[2]=='zlo' and row[3]=='zhi': 113 self.zdimension=[row[0], row[1]] 114 self.keywords.append('zlo zhi') 115 else: 116 # Set Sectflag to 2, assume line is a keyword 117 sectflag=2 118 checkkey=1 #ensures keyword will be checked in body portion 119 keyword=row 120 elif len(row)==5: 121 if row[1]=='extra' and row[2]=='bond' and row[3]=='per' and row[4]=='atom': 122 self.extrabonds=row[0] 123 self.keywords.append('extra bond per atom') 124 else: 125 # Set Sectflag to 2, assume line is a keyword 126 sectflag=2 127 checkkey=1 #ensures keyword will be checked in body portion 128 keyword=row 129 elif len(row)==6: 130 if row[3]=='xy' and row[4]=='xz' and row[5]=='yz': 131 self.tilt=[row[0], row[1], row[2]] 132 self.keywords.append('xy xz yz') 133 else: 134 # Set Sectflag to 2, assume line is a keyword 135 sectflag=2 136 checkkey=1 #ensures keyword will be checked in body portion 137 keyword=row 138 else: 139 # set sectflag to 2, assume line is a keyword 140 sectflag=2 141 checkkey=1 #ensures keyword will be checked in body portion 142 keyword=row 143 elif sectflag==2: 144 if checkkey==1: 145 if len(keyword)==1: 146 if keyword[0]=='Atoms' or keyword[0]=='Velocities' or keyword[0]=='Masses' or\ 147 keyword[0]=='Shapes' or keyword[0]=='Dipoles' or keyword[0]=='Bonds' or\ 148 keyword[0]=='Angles' or keyword[0]=='Dihedrals' or keyword[0]=='Impropers': 149 bodyflag=1 150 blanknum=0 151 self.keywords.append(keyword[0]) 152 checkkey=0 153 else: 154 bodyflag=0 155 checkkey=0 156 elif len(keyword)==2: 157 if row[1]=='Coeffs' and (row[0]=='Pair' or row[0]=='Bond' or row[0]=='Angle' or\ 158 row[0]=='Dihedral' or row[0]=='Improper'):#class 2 force field keywords not included 159 bodyflag=1 160 blanknum=0 161 self.keywords.append('{0} {1}'.format(keyword[0],keyword[1])) 162 checkkey=0 163 else: 164 bodyflag=0 165 checkkey=0 166 else: 167 #egnore line and set bodyflag=0 168 bodyflag=0 169 checkkey=0 170 if bodyflag==0: #bodyflag 0 means no body keyword has been found 171 if len(row)==1: 172 if row[0]=='Atoms' or row[0]=='Velocities' or row[0]=='Masses' or\ 173 row[0]=='Shapes' or row[0]=='Dipoles' or row[0]=='Bonds' or\ 174 row[0]=='Angles' or row[0]=='Dihedrals' or row[0]=='Impropers': 175 bodyflag=1 176 blanknum=0 177 keyword=row 178 self.keywords.append(keyword[0]) 179 else: 180 #egnore line 181 pass 182 elif len(row)==2: 183 if row[1]=='Coeffs' and (row[0]=='Pair' or row[0]=='Bond' or row[0]=='Angle' or\ 184 row[0]=='Dihedral' or row[0]=='Improper'):#class 2 force field keywords not included 185 bodyflag=1 186 blanknum=0 187 keyword=row 188 self.keywords.append('{0} {1}'.format(keyword[0],keyword[1])) 189 else: 190 #egnore line 191 pass 192 else: 193 #egnore line 194 pass 195 elif bodyflag==1: #currently assumes 1 or more blank lines are between body data keywords 196 if len(row)==0: 197 blanknum+=1 198 if blanknum>1: 199 bodyflag=0 200 elif len(keyword)==1: 201 if keyword[0]=='Atoms': 202 try: 203 int(row[0]) 204 except ValueError: 205 keyword=row 206 checkkey=1 207 else: 208 self.atoms.append(row) 209 elif keyword[0]=='Velocities': 210 try: 211 int(row[0]) 212 except ValueError: 213 keyword=row 214 checkkey=1 215 else: 216 self.velocities.append(row) 217 elif keyword[0]=='Masses': 218 try: 219 int(row[0]) 220 except ValueError: 221 keyword=row 222 checkkey=1 223 else: 224 self.masses.append(row) 225 elif keyword[0]=='Shapes': 226 try: 227 int(row[0]) 228 except ValueError: 229 keyword=row 230 checkkey=1 231 else: 232 self.shapes.append(row) 233 elif keyword[0]=='Dipoles': 234 try: 235 int(row[0]) 236 except ValueError: 237 keyword=row 238 checkkey=1 239 else: 240 self.dipoles.append(row) 241 elif keyword[0]=='Bonds': 242 try: 243 int(row[0]) 244 except ValueError: 245 keyword=row 246 checkkey=1 247 else: 248 self.bonds.append(row) 249 elif keyword[0]=='Angles': 250 try: 251 int(row[0]) 252 except ValueError: 253 keyword=row 254 checkkey=1 255 else: 256 self.angles.append(row) 257 elif keyword[0]=='Dihedrals': 258 try: 259 int(row[0]) 260 except ValueError: 261 keyword=row 262 checkkey=1 263 else: 264 self.dihedrals.append(row) 265 elif keyword[0]=='Impropers': 266 try: 267 int(row[0]) 268 except ValueError: 269 keyword=row 270 checkkey=1 271 else: 272 self.impropers.append(row) 273 else: 274 #egnore line and change bodyflag to 0 275 bodyflag=0 276 elif len(keyword)==2: 277 if keyword[0]=='Pair' and keyword[1]=='Coeffs': 278 try: 279 int(row[0]) 280 except ValueError: 281 keyword=row 282 checkkey=1 283 else: 284 self.paircoef.append(row) 285 elif keyword[0]=='Bond' and keyword[1]=='Coeffs': 286 try: 287 int(row[0]) 288 except ValueError: 289 keyword=row 290 checkkey=1 291 else: 292 self.bondcoef.append(row) 293 elif keyword[0]=='Angle' and keyword[1]=='Coeffs': 294 try: 295 int(row[0]) 296 except ValueError: 297 keyword=row 298 checkkey=1 299 else: 300 self.anglecoef.append(row) 301 elif keyword[0]=='Dihedral' and keyword[1]=='Coeffs': 302 try: 303 int(row[0]) 304 except ValueError: 305 keyword=row 306 checkkey=1 307 else: 308 self.dihedralcoef.append(row) 309 elif keyword[0]=='Improper' and keyword[1]=='Coeffs': 310 try: 311 int(row[0]) 312 except ValueError: 313 keyword=row 314 checkkey=1 315 else: 316 self.impropercoef.append(row) 317 else: 318 #egnore line and change bodyflag to 0 319 bodyflag=0 320 else: 321 #egnore line and change bodyflag to 0 322 bodyflag=0 323 f.close() 324 325 def write(self,file,modflag): 326 """Write lammps data files using the lammps keywords and lammpsdata structures 327 writes first line of the file as a Comment line 328 if no modifications to any of the lammpsdata structures (modflag=0) 329 Use Keywords to write lammpsdata structures directly 330 if modifications to any of the lammpsdata structures (modflag=1) 331 Key Lammpsdata structures like atom numbers, coefficient numbers 332 need to be modified to match the other modified lammpsdata structures 333 This section will still use the keywords to write lammpsdata structures. 334 For all modflags, the code will write data to the file until all of the 335 keyword's data structures have been finished writing. 336 The keywords are stored in self.keywords""" 337 f=open(file,'w') 338 f.write('polymer data file\n') 339 for row in self.keywords: 340 if row=='atoms': 341 if modflag==0: 342 f.write('{0} {1}\n'.format(self.atomnum,row)) 343 elif modflag==1: 344 f.write('{0} {1}\n'.format(len(self.atoms),row)) 345 elif row=='bonds': 346 if modflag==0: 347 f.write('{0} {1}\n'.format(self.bondnum,row)) 348 elif modflag==1: 349 f.write('{0} {1}\n'.format(len(self.bonds),row)) 350 elif row=='angles': 351 if modflag==0: 352 f.write('{0} {1}\n'.format(self.anglenum,row)) 353 elif modflag==1: 354 f.write('{0} {1}\n'.format(len(self.angles),row)) 355 elif row=='dihedrals': 356 if modflag==0: 357 f.write('{0} {1}\n'.format(self.dihedralnum,row)) 358 elif modflag==1: 359 f.write('{0} {1}\n'.format(len(self.dihedrals),row)) 360 elif row=='impropers': 361 if modflag==0: 362 f.write('{0} {1}\n'.format(self.impropernum,row)) 363 elif modflag==1: 364 f.write('{0} {1}\n'.format(len(self.impropers),row)) 365 elif row=='atom types': 366 if modflag==0: 367 f.write('{0} {1}\n'.format(self.atomtypenum,row)) 368 elif modflag==1: 369 f.write('{0} {1}\n'.format(len(self.masses),row)) 370 elif row=='bond types': 371 if modflag==0: 372 f.write('{0} {1}\n'.format(self.bondtypenum,row)) 373 elif modflag==1: 374 f.write('{0} {1}\n'.format(len(self.bondcoef),row)) 375 elif row=='angle types': 376 if modflag==0: 377 f.write('{0} {1}\n'.format(self.angletypenum,row)) 378 elif modflag==1: 379 f.write('{0} {1}\n'.format(len(self.anglecoef),row)) 380 elif row=='dihedral types': 381 if modflag==0: 382 f.write('{0} {1}\n'.format(self.dihedraltypenum,row)) 383 elif modflag==1: 384 f.write('{0} {1}\n'.format(len(self.dihedralcoef),row)) 385 elif row=='improper types': 386 if modflag==0: 387 f.write('{0} {1}\n'.format(self.impropertypenum,row)) 388 elif modflag==1: 389 f.write('{0} {1}\n'.format(len(self.impropercoef),row)) 390 elif row=='xlo xhi': 391 f.write('{0} {1} {2}\n'.format(self.xdimension[0],self.xdimension[1],row)) 392 elif row=='ylo yhi': 393 f.write('{0} {1} {2}\n'.format(self.ydimension[0],self.ydimension[1],row)) 394 elif row=='zlo zhi': 395 f.write('{0} {1} {2}\n'.format(self.zdimension[0],self.zdimension[1],row)) 396 elif row=='extra bond per atom': 397 f.write('{0} {1}\n'.format(self.extrabonds,row)) 398 elif row=='xy xz yz': 399 f.write('{0} {1} {2} {3}\n'.format(self.tilt[0],self.tilt[1],self.tilt[2],row)) 400 elif row=='Atoms': 401 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 402 f.write('\n') #new line between body keyword and body data 403 for line in self.atoms: 404 f.write('\n') #creates a new line for adding body data 405 for item in line: 406 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 407 f.write('\n') #allows space to be added between the end of body data and a new keyword 408 elif row=='Velocities': 409 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 410 f.write('\n') #new line between body keyword and body data 411 for line in self.velocities: 412 f.write('\n') #creates a new line for adding body data 413 for item in line: 414 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 415 f.write('\n') #allows space to be added between the end of body data and a new keyword 416 elif row=='Masses': 417 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 418 f.write('\n') #new line between body keyword and body data 419 for line in self.masses: 420 f.write('\n') #creates a new line for adding body data 421 for item in line: 422 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 423 f.write('\n') #allows space to be added between the end of body data and a new keyword 424 elif row=='Shapes': 425 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 426 f.write('\n') #new line between body keyword and body data 427 for line in self.shapes: 428 f.write('\n') #creates a new line for adding body data 429 for item in line: 430 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 431 f.write('\n') #allows space to be added between the end of body data and a new keyword 432 elif row=='Dipoles': 433 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 434 f.write('\n') #new line between body keyword and body data 435 for line in self.dipoles: 436 f.write('\n') #creates a new line for adding body data 437 for item in line: 438 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 439 f.write('\n') #allows space to be added between the end of body data and a new keyword 440 elif row=='Bonds': 441 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 442 f.write('\n') #new line between body keyword and body data 443 for line in self.bonds: 444 f.write('\n') #creates a new line for adding body data 445 for item in line: 446 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 447 f.write('\n') #allows space to be added between the end of body data and a new keyword 448 elif row=='Angles': 449 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 450 f.write('\n') #new line between body keyword and body data 451 for line in self.angles: 452 f.write('\n') #creates a new line for adding body data 453 for item in line: 454 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 455 f.write('\n') #allows space to be added between the end of body data and a new keyword 456 elif row=='Dihedrals': 457 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 458 f.write('\n') #new line between body keyword and body data 459 for line in self.dihedrals: 460 f.write('\n') #creates a new line for adding body data 461 for item in line: 462 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 463 elif row=='Impropers': 464 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 465 f.write('\n') #new line between body keyword and body data 466 for line in self.impropers: 467 f.write('\n') #creates a new line for adding body data 468 for item in line: 469 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 470 f.write('\n') #allows space to be added between the end of body data and a new keyword 471 elif row=='Pair Coeffs': 472 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 473 f.write('\n') #new line between body keyword and body data 474 for line in self.paircoef: 475 f.write('\n') #creates a new line for adding body data 476 for item in line: 477 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 478 f.write('\n') #allows space to be added between the end of body data and a new keyword 479 elif row=='Bond Coeffs': 480 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 481 f.write('\n') #new line between body keyword and body data 482 for line in self.bondcoef: 483 f.write('\n') #creates a new line for adding body data 484 for item in line: 485 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 486 f.write('\n') #allows space to be added between the end of body data and a new keyword 487 elif row=='Angle Coeffs': 488 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 489 f.write('\n') #new line between body keyword and body data 490 for line in self.anglecoef: 491 f.write('\n') #creates a new line for adding body data 492 for item in line: 493 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 494 f.write('\n') #allows space to be added between the end of body data and a new keyword 495 elif row=='Dihedral Coeffs': 496 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 497 f.write('\n') #new line between body keyword and body data 498 for line in self.dihedralcoef: 499 f.write('\n') #creates a new line for adding body data 500 for item in line: 501 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 502 f.write('\n') #allows space to be added between the end of body data and a new keyword 503 elif row=='Improper Coeffs': 504 f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords 505 f.write('\n') #new line between body keyword and body data 506 for line in self.impropercoef: 507 f.write('\n') #creates a new line for adding body data 508 for item in line: 509 f.write(' {0}'.format(item)) #adds in each piece of body data with space imbetween 510 f.write('\n') #allows space to be added between the end of body data and a new keyword 511 else: 512 pass 513 f.close() 514 515 def atomorder(self): 516 """Takes self.atoms and organizes the atom id from least to greatest. 517 If the atom ids are already ordered this algorithm will do nothing.""" 518 current=range(len(self.atoms[0])) # initialize current [assumes self.atoms coloumn #'s does not change] 519 for i in range(1,len(self.atoms)): #when i=0, self.atoms will not change; therefore its skipped 520 for k in range(len(self.atoms[i])): 521 current[k]=self.atoms[i][k] 522 for j in range(i-1,-1,-1): 523 for k in range(len(self.atoms[j])): 524 self.atoms[j+1][k]=self.atoms[j][k] 525 if int(current[0]) > int(self.atoms[j][0]): 526 for k in range(len(current)): 527 self.atoms[j+1][k]=current[k] 528 break 529 elif j==0: 530 for k in range(len(current)): 531 self.atoms[j][k]=current[k] 532 #dont need a break here because this is the last j value in the for loop 533 534 def addatoms(self, atoms, retlist=False): 535 """Appends atoms to self.atoms. Assumes atoms are written in the correct atomtype format. 536 Change added atom numbers in self.atoms so self.atoms will be in increasing sequential order. 537 If retlist is True return a list of the modified atom numbers otherwise exit.""" 538 initpos=len(self.atoms) #Store index of first appended atoms 539 for item in atoms: 540 self.atoms.append(range(len(item))) #initializing spots for the new atoms to go 541 numberchange=[] #initiate number change where the changed atom numbers will be stored 542 count=0 543 for i in range(initpos,len(self.atoms)): 544 for j in range(len(self.atoms[i])): 545 if j==0: 546 self.atoms[i][0]=str(i+1) 547 numberchange.append(str(i+1)) 548 else: 549 self.atoms[i][j]=atoms[count][j] 550 count+=1 551 if retlist: return numberchange 552 else: return #need to redo this algorithm so their is no referencing going on. 553 554 def adddata(self,data,keyword): 555 """Adds data to a keyword's existing data structure 556 All body keywords are viable except Atoms which has its own method 557 All header keywords will do nothing in this algrorithm because they have their own algorithm 558 changes the body id number so id numbers will be in sequential order""" 559 if keyword=='Velocities': 560 pos=len(self.velocities) 561 for item in data: 562 self.velocities.append(item) #adds data to existing data structure 563 self.velocities[pos][0]=str(pos+1) #changing body id number 564 pos+=1 565 elif keyword=='Masses': 566 pos=len(self.masses) 567 for item in data: 568 self.masses.append(item) #adds data to existing data structure 569 self.masses[pos][0]=str(pos+1) #changing body id number 570 pos+=1 571 elif keyword=='Shapes': 572 pos=len(self.shapes) 573 for item in data: 574 self.shapes.append(item) #adds data to existing data structure 575 self.shapes[pos][0]=str(pos+1) #changing body id number 576 pos+=1 577 elif keyword=='Dipoles': 578 pos=len(self.dipoles) 579 for item in data: 580 self.dipoles.append(item) #adds data to existing data structure 581 self.dipoles[pos][0]=str(pos+1) #changing body id number 582 pos+=1 583 elif keyword=='Bonds': 584 pos=len(self.bonds) 585 for item in data: 586 self.bonds.append(item) #adds data to existing data structure 587 self.bonds[pos][0]=str(pos+1) 588 pos+=1 589 elif keyword=='Angles': 590 pos=len(self.angles) 591 for item in data: 592 self.angles.append(item) #adds data to existing data structure 593 self.angles[pos][0]=str(pos+1) #changing body id number 594 pos+=1 595 elif keyword=='Dihedrals': 596 pos=len(self.dihedrals) 597 for item in data: 598 self.dihedrals.append(item) #adds data to existing data structure 599 self.dihedrals[pos][0]=str(pos+1) #changing body id number 600 pos+=1 601 elif keyword=='Impropers': 602 pos=len(self.impropers) 603 for item in data: 604 self.impropers.append(item) #adds data to existing data structure 605 self.impropers[pos][0]=str(pos+1) #changing body id number 606 pos+=1 607 elif keyword=='Pair Coeffs': 608 pos=len(self.paircoef) 609 for item in data: 610 self.paircoef.append(item) #adds data to existing data structure 611 self.paircoef[pos][0]=str(pos+1) #changing body id number 612 pos+=1 613 elif keyword=='Bond Coeffs': 614 pos=len(self.bondcoef) 615 for item in data: 616 self.bondcoef.append(item) #adds data to existing data structure 617 self.bondcoef[pos][0]=str(pos+1) #changing body id number 618 pos+=1 619 elif keyword=='Angle Coeffs': 620 pos=len(self.anglecoef) 621 for item in data: 622 self.anglecoef.append(item) #adds data to existing data structure 623 self.anglecoef[pos][0]=str(pos+1) #changing body id number 624 pos+=1 625 elif keyword=='Dihedral Coeffs': 626 pos=len(self.dihedralcoef) 627 for item in data: 628 self.dihedralcoef.append(item) #adds data to existing data structure 629 self.dihedralcoef[pos][0]=str(pos+1) #changing body id number 630 pos+=1 631 elif keyword=='Improper Coeffs': 632 pos=len(self.impropercoef) 633 for item in data: 634 self.impropercoef.append(item) #adds data to existing data structure 635 self.impropercoef[pos][0]=str(pos+1) #changing body id number 636 pos+=1 637 638# def modifiydata(data,keyword): Will not implement 639# """for modifying header data""" 640 641 def changeatomnum(self,data,originalatoms,atomchanges,vflag=False): 642 """Takes data which contains atom ids 643 and alters those ids from the original atom id system in originalatoms 644 to the new atom id system in atomchanges.""" 645 #set up boolean array to match the size of data and with all True values 646# print atomchanges 647 array=booleanarray(len(data),len(data[0]),True) 648# print originalatoms 649 if vflag: # for changing atomnumbers for velocities 650 for i in range(len(originalatoms)): #len of originalatoms should match len of atomchanges 651 if originalatoms[i][0]==atomchanges[i]: continue 652 for j in range(len(data)): 653 for k in range(1): 654 if data[j][k]==originalatoms[i][0]: #checks if the atom id in data 655 #matches the atom id in origianalatoms 656 if array.getelement(j,k): #if the boolean array is true 657 #set the boolean array to False and 658 #change the atom id in data to atomchanges 659 array.setelement(j,k,False) 660 data[j][k]=atomchanges[i] 661 break 662 #the change of boolean array to False insures 663 #the atom id in data will only be changed once. 664 else: # for changing atom numbers for everything else 665 for i in range(len(originalatoms)): #len of originalatoms should match len of atomchanges 666 if originalatoms[i][0]==atomchanges[i]: continue 667 for j in range(len(data)): 668 for k in range(2,len(data[j])): 669 if data[j][k]==originalatoms[i][0]: #checks if the atom id in data 670 #matches the atom id in origianalatoms 671 if array.getelement(j,k): #if the boolean array is true 672 #set the boolean array to False and 673 #change the atom id in data to atomchanges 674 array.setelement(j,k,False) 675 data[j][k]=atomchanges[i] 676 break 677 #the change of boolean array to False insures 678 #the atom id in data will only be changed once. 679# print data 680 return data 681 682 def deletebodydata(self,keyword): 683 """Changes a keyword's class data structure to an empty structure""" 684 if keyword=='Velocities': 685 self.velocities=[] 686 elif keyword=='Masses': 687 self.masses=[] 688 elif keyword=='Shapes': 689 self.shapes=[] 690 elif keyword=='Dipoles': 691 self.dipoles=[] 692 elif keyword=='Bonds': 693 self.bonds=[] 694 elif keyword=='Angles': 695 self.angles=[] 696 elif keyword=='Dihedrals': 697 self.dihedrals=[] 698 elif keyword=='Impropers': 699 self.impropers=[] 700 elif keyword=='Pair Coeffs': 701 self.paircoef=[] 702 elif keyword=='Bond Coeffs': 703 self.bondcoef=[] 704 elif keyword=='Angle Coeffs': 705 self.anglecoef=[] 706 elif keyword=='Dihedral Coeffs': 707 self.dihedralcoef=[] 708 elif keyword=='Improper Coeffs': 709 self.impropercoef=[] 710 elif keyword=='Atoms': 711 self.atoms=[] 712 713 def extractmolecules(self,molecule): 714 """Takes the variable molecule and 715 extracts the individual molecules' data back into lmpsdata. 716 This extraction takes place through a 4 step process. 717 Step 1: Use a molecule's keywords to alter the lmpsdata data structures to empty. 718 To accomplish this procedure use the lmpsdata method deletebodydata. 719 Step 2: Add the molecules' atoms to lmpsdata's atoms using the method addatoms. 720 Return a list of atom id changes for each molecule. 721 Step 3: Utilize each molecules list of atom id changes to change their data's atom id numbers. 722 Uses changeatomnum and returns the altered molecule's data. 723 Step 4: Add the altered molecules' data to lmpsdata's data using the method adddata""" 724 #Use molecule index 0's keyword to change the equivalent lmpsdata structures to empty 725 print 'extracting the molecules back to data' 726 for keyword in molecule[0].keywords: # step 1 727 self.deletebodydata(keyword) 728 atomchanges=[] #initializing step 2 729 for i in range(len(molecule)): # step 2 730 atomchanges.append(self.addatoms(molecule[i].atoms,True)) 731 for i in range(len(molecule)): # step 3 732 for keyword in molecule[i].keywords: 733 if keyword=='Angles': 734 molecule[i].angles=self.changeatomnum(molecule[i].angles,molecule[i].atoms,atomchanges[i]) 735 if keyword=='Bonds': 736 molecule[i].bonds=self.changeatomnum(molecule[i].bonds,molecule[i].atoms,atomchanges[i]) 737 elif keyword=='Dihedrals': 738 molecule[i].dihedrals=self.changeatomnum(molecule[i].dihedrals,molecule[i].atoms,atomchanges[i]) 739 elif keyword=='Velocities': 740 molecule[i].velocities=self.changeatomnum(molecule[i].velocities,molecule[i].atoms,atomchanges[i],True) 741 elif keyword=='Impropers': 742 molecule[i].impropers=self.changeatomnum(molecule[i].impropers,molecule[i].atoms,atomchanges[i]) 743 for i in range(len(molecule)): # step 4 744 for keyword in molecule[i].keywords: 745 if keyword=='Angles': 746 self.adddata(molecule[i].angles,keyword) 747 elif keyword=='Bonds': 748 self.adddata(molecule[i].bonds,keyword) 749 elif keyword=='Dihedrals': 750 self.adddata(molecule[i].dihedrals,keyword) 751 elif keyword=='Velocities': 752 self.adddata(molecule[i].velocities,keyword) 753 elif keyword=='Impropers': 754 self.adddata(molecule[i].impropers,keyword) 755 756 def density(self,ringsize,init,final,file=' '): 757 """Create spherical shells from the initial radius to the final radius 758 The spherical shell's thickness is defined by ringsize 759 Calculate the mass of particles within each spherical shell 760 Calculate the volume of each spherical shell 761 Then calculate the density in each spherical shell 762 Current code is written with the assumption the spheres are centered around the origin 763 If a file is listed place the results in the variable file otherwise return the results 764 The distance calculations are written assuming image flags are in the data file""" 765 rho=[] 766 r=[init] #initial radius to examine 767 radius=init+ringsize 768 while radius<final: #finding the rest of the radii to examine 769 r.append(radius) 770 radius+=ringsize 771 for radius in r: 772 volume=4.0/3.0*3.1416*((radius+ringsize)**3 - radius**3) 773 mass=0.0 774 for row in self.atoms: 775 if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\ 776 self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\ 777 self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\ 778 self.atomtype=='peri': 779 l=len(row)-1 780 dist=float(row[l-3])**2+float(row[l-4])**2+float(row[l-5])**2 781 if dist<(radius+ringsize)**2 and dist>=radius**2: 782 if self.atomtype=='atomic' or self.atomtype=='charge' or\ 783 self.atomtype=='colloid' or self.atomtype=='electron' or\ 784 self.atomtype=='granular' or self.atomtype=='peri': 785 type=int(row[1])-1 786 else: 787 type=int(row[2])-1 788 mass+=float(self.masses[type][1]) 789 #assumes self.masses is written in atomtype order 790 elif self.atomtype=='dipole': 791 l=len(row)-1 792 dist=float(row[l-6])**2+float(row[l-7])**2+float(row[l-8])**2 793 if dist<(radius+ringsize)**2 and dist>=radius**2: 794 type=int(row[1])-1 795 mass+=float(self.masses[type][1]) 796 #assumes self.masses is written in atomtype order 797 elif self.atomtype=='ellipsoid': 798 l=len(row)-1 799 dist=float(row[l-7])**2+float(row[l-8])**2+float(row[l-9])**2 800 if dist<(radius+ringsize)**2 and dist>=radius**2: 801 type=int(row[1]) 802 mass+=float(self.masses[type][1]) 803 #assumes self.masses is written in atomtype order 804 elif self.atomtype=='hybrid': 805 dist=float(row[2])**2+float(row[3])**2+float(row[4])**2 806 if dist<(radius+ringsize)**2 and dist>=radius**2: 807 type=int(row[1]) 808 mass+=float(self.masses[type][1]) 809 #assumes self.masses is written in atomtype order 810 rho.append(mass/volume) 811 if file==' ': 812 return r,rho 813 else: 814 f=open(file,'w') 815 for i in range(len(r)): 816 f.write('{0} {1}\n'.format(r[i],rho[i])) 817 f.close() 818 return 819 820 def createxyz(self,file, routine='mass', values=None): 821 """Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user. 822 The mass version is assessed by setting the routine to 'mass' which is the default method. 823 The other version is assessed by setting the routine to 'atomtype'. 824 The other version takes values which is a list containing the value the user wants to assign those atomtypes to. 825 The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1. 826 This makes it easier to find the atomtype and assign the value from the list 827 All atom data is assumed to have image flags in the data.""" 828 f=open(file,'w') 829 f.write('{0}\n'.format(len(self.atoms))) 830 f.write('atoms\n') 831 if routine=='mass': 832 for line in self.atoms: 833 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 834 self.atomtype=='molecular': 835 type=int(line[2])-1 #3rd position 836 else: 837 type=int(line[1])-1 #2nd position 838 mass=self.masses[type][1] 839 if mass=='12.0107': elementnum=6 840 elif mass=='15.9994': elementnum=8 841 elif mass=='26.9815':elementnum=13 842 else: 843 print 'no matching mass value. The method will exit' 844 return 845 l=len(line)-1 846 if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\ 847 self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\ 848 self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\ 849 self.atomtype=='peri': 850 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 851 elif self.atomtype=='dipole': 852 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6])) 853 elif self.atomtype=='ellipsoid': 854 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7])) 855 elif self.atomtype=='hybrid': 856 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4])) 857 f.close() 858 elif routine=='atomtype': 859 for line in self.atoms: 860 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 861 self.atomtype=='molecular': 862 type=int(line[2])-1 #3rd position 863 else: 864 type=int(line[1])-1 #2nd position 865 elementnum=values[type] 866 l=len(line)-1 867 if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\ 868 self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\ 869 self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\ 870 self.atomtype=='peri': 871 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 872 elif self.atomtype=='dipole': 873 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6])) 874 elif self.atomtype=='ellipsoid': 875 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7])) 876 elif self.atomtype=='hybrid': 877 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4])) 878 f.close() 879 880 881class booleanarray: 882 """A class that stores boolean values in a list of lists.""" 883 def __init__(self,rownum, colnum, initval): 884 """ initialize a list of lists (array) with 885 rownum correspondinig to the number of lists in the list and 886 colnum corresponding to the number of elements in the list's list. 887 initval is the value the list of lists will be initialized with. 888 initval should be a boolean value.""" 889 # initializing self.array 890 self.array=[] 891 for i in range(rownum): 892 self.array.append(range(colnum)) 893 # setting elements of self.array to initval 894 for i in range(rownum): 895 for j in range(colnum): 896 self.setelement(i,j,initval) 897 898 def setelement(self,rownum, colnum, value): 899 """Assigns value to the list of lists (array) element at rownum and colnum.""" 900 self.array[rownum][colnum]=value 901 902 def getelement(self,rownum,colnum): 903 """Returns element in the list of lists (array) at rownum and colnum.""" 904 return self.array[rownum][colnum] 905 906def atomdistance(a,b,atomtype): 907 """Returns the distance between atom a and b. 908 Atomtype is the style the atom data is written in. 909 All atom data is assumed to have image flags in the data. 910 Atom a and atom b are assumed to be the same atomtype.""" 911 from math import sqrt 912 if atomtype=='angle' or atomtype=='atomic' or atomtype=='bond' or\ 913 atomtype=='charge' or atomtype=='colloid' or atomtype=='electron' or\ 914 atomtype=='full' or atomtype=='granular' or atomtype=='molecular' or\ 915 atomtype=='peri': 916 l=len(a)-1 917 dist=sqrt((float(a[l-3])-float(b[l-3]))**2\ 918 +(float(a[l-4])-float(b[l-4]))**2\ 919 +(float(a[l-5])-float(b[l-5]))**2) 920 elif atomtype=='dipole': 921 l=len(a)-1 922 dist=sqrt((float(a[l-6])-float(b[l-6]))**2\ 923 +(float(a[l-7])-float(b[l-7]))**2\ 924 +(float(a[l-8])-float(b[l-8]))**2) 925 elif atomtype=='ellipsoid': 926 l=len(a)-1 927 dist=sqrt((float(a[l-7])-float(b[l-7]))**2\ 928 +(float(a[l-8])-float(b[l-8]))**2\ 929 +(float(a[l-9])-float(b[l-9]))**2) 930 elif atomtype=='hybrid': 931 dist=sqrt((float(a[2])-float(b[2]))**2\ 932 +(float(a[3])-float(b[3]))**2\ 933 +(float(a[4])-float(b[4]))**2) 934 return dist 935 936def distance(a,coord,atomtype): 937 """Returns the distance between atom a and coord. 938 Atomtype is the style the atom data is written in. 939 All atom data is assumed to have image flags in the data.""" 940 from math import sqrt #need to alter this slightly 941 if atomtype=='angle' or atomtype=='atomic' or atomtype=='bond' or\ 942 atomtype=='charge' or atomtype=='colloid' or atomtype=='electron' or\ 943 atomtype=='full' or atomtype=='granular' or atomtype=='molecular' or\ 944 atomtype=='peri': 945 l=len(a)-1 946 dist=sqrt((float(a[l-3])-coord[2])**2\ 947 +(float(a[l-4])-coord[1])**2\ 948 +(float(a[l-5])-coord[0])**2) 949 elif atomtype=='dipole': 950 l=len(a)-1 951 dist=sqrt((float(a[l-6])-coord[2])**2\ 952 +(float(a[l-7])-coord[1])**2\ 953 +(float(a[l-8])-coord[0])**2) 954 elif atomtype=='ellipsoid': 955 l=len(a)-1 956 dist=sqrt((float(a[l-7])-coord[2])**2\ 957 +(float(a[l-8])-coord[1])**2\ 958 +(float(a[l-9])-coord[0])**2) 959 elif atomtype=='hybrid': 960 dist=sqrt((float(a[2])-coord[0])**2\ 961 +(float(a[3])-coord[1])**2\ 962 +(float(a[4])-coord[0])**2) 963 return dist 964 965 966class particlesurface: 967 def __init__(self,particle,cutoff,atomid,atomtype,shape='sphere'): 968 """Builds a particle surface with a specific shape from a particle 969 The atoms chosen from the surface will have the specific atomid 970 atomid will be given in terms of an integer and not a string.""" 971 self.particle=particle.atoms 972 self.atomtype=atomtype 973 self.cutoff=cutoff 974 self.surface=[] 975 self.particlelocator=[] #stores surface particle locations with respect to the particle 976 self.deletedatoms=[] 977 if shape=='sphere': self.createspheresurf(atomid) 978 979 def createspheresurf(self,atomid): 980 """Assumes sphere is centered around 0,0,0. 981 Finds maximum distance particle with atomid. 982 Than all atoms within cutoff distance of max distance 983 will be included into self.surface 984 Also will build a list of where the surface particles are 985 in relationship to the particle. 986 Will create a booleanarray to keep track of any changes made to the surface.""" 987 particledist=[] 988 for row in self.particle: #Stores all of the distances of the particle in particledist 989 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 990 self.atomtype=='molecular': 991 type=int(row[2]) #3rd position 992 else: 993 type=int(row[1]) #2nd position 994 if type==atomid: #only calculates the distance if the atom has the correct atomid number 995 particledist.append(distance(row,[0,0,0],self.atomtype)) 996 else: 997 particledist.append(0.0) 998 self.maxdist=max(particledist)# finds the max dist. 999 for i in range(len(particledist)): 1000 if particledist[i]>=(self.maxdist-self.cutoff): 1001 self.particlelocator.append(i) 1002 self.surface.append(self.particle[i]) 1003 self.bool=booleanarray(len(self.surface),1,True) 1004 1005 def getsurfatom(self,rownum): 1006 return self.surface[rownum] 1007 1008 def getbool(self,rownum): 1009 return self.bool.getelement(rownum,0) 1010 1011 def setbool(self,rownum,value): 1012 self.bool.setelement(rownum,0,value) 1013 1014 def removesurfatom(self,rownum): 1015 """note this does not actually remove the surface atom. 1016 Instead this adds a value to the list deletedatoms. 1017 This list is later used in extractparticle to actually delete those atoms""" 1018 self.deletedatoms.append(self.particlelocator[rownum]) 1019 1020 def extractparticle(self): 1021 """"deletesatoms from particle from the list of deletedatoms 1022 and than returns the particle.""" 1023 self.particle=self.deleteatoms(self.particle,self.deletedatoms) 1024 return self.particle 1025 1026 def deleteatoms(self,structure,rows): 1027 """delete atoms from particle and shifts the structure down""" 1028 new=[] 1029 #multiple copying of b to the rows being replaced. 1030 if rows==[]: 1031 for line in structure: 1032 new.append(line) #if no rows need replacing copy structure 1033 return new 1034 for i in range(rows[0]): 1035 new.append(structure[i]) #copy structure to new until the first replaced row 1036 count=0 1037 for i in range(rows[0],len(structure)-len(rows)):# to replace rows and shift undeleted rows over 1038 for j in range(i+1+count,len(structure)): 1039 for val in rows: 1040 if val==j: 1041 count+=1 1042 break 1043 if val==j:continue 1044 else: 1045 new.append(structure[j]) 1046 break 1047 return new 1048 1049 def addatom(self,atomtype,charge=None,moleculenum=None): 1050 """Adds an atom to the particle surface between maxdist and the cutoff distance. 1051 This atom is stored in self.particle and self.surface. 1052 In the current coding the particle is centered at (0,0,0). 1053 This method has a required input of atomtype and optional inputs of charge and moleculenum. 1054 This method currently does not support correctly self.atomtype of dipole, electron, ellipsoid, granular, peri or hybrid 1055 Will use the image flags 0,0,0. 1056 Note: The atomtype here is different from the atomtype attribute as part of this class 1057 Note: If the added atom will be required for bonding later than you will need to extract the particle data 1058 and rebuild the surface, because this method doesn't include updates to the required variables due to some coding issues""" 1059 import math, random 1060 # Checks to make sure self.atomtype is not set to an unsupported value 1061 if self.atomtype=='dipole' or self.atomtype=='electron' or\ 1062 self.atomtype=='ellipsoid' or self.atomtype=='peri' or self.atomtype=='hybrid': 1063 print self.atomtype, 'is not currently supported. But, you can add support by adding the needed functionality and inputs' 1064 return 1065 1066 print 'adding an atom to the surface' 1067 # Randomly assigns the position of a new atom in a spherical shell between self.cutoff and self.maxdist 1068 foundpoint=False 1069 while (foundpoint==False): 1070 x=random.uniform(-self.maxdist,self.maxdist) 1071 y=random.uniform(-self.maxdist,self.maxdist) 1072 try: 1073 z=random.uniform(math.sqrt(self.cutoff**2-x**2-y**2),math.sqrt(self.maxdist**2-x**2-y**2)) 1074 except ValueError: 1075 continue 1076 distance=math.sqrt(x**2+y**2+z**2) 1077 if distance<self.maxdist and distance>=self.maxdist-self.cutoff: 1078 foundpoint=True 1079 1080 # initialize new row in self.particle and self.surface 1081 self.particle.append([]) 1082 self.surface.append([]) 1083 1084 pposition=len(self.particle)-1 #particle position 1085 sposition=len(self.surface)-1 #surface position 1086 1087 # Adds the atom id number to the new row 1088 self.particle[pposition].append(str(pposition+1)) 1089 self.surface[sposition].append(str(pposition+1)) 1090 1091 # Adds the atomtype and the moleculenum if needed to the new row 1092 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 1093 self.atomtype=='molecular': 1094 self.particle[pposition].append(str(moleculenum)) 1095 self.surface[sposition].append(str(moleculenum)) 1096 self.particle[pposition].append(str(atomtype)) 1097 self.surface[sposition].append(str(atomtype)) 1098 else: 1099 self.particle[pposition].append(str(atomtype)) 1100 self.surface[sposition].append(str(atomtype)) 1101 1102 # Adds the charge if needed to the new row 1103 if self.atomtype=='charge' or self.atomtype=='full': 1104 self.particle[pposition].append(str(charge)) 1105 self.surface[sposition].append(str(charge)) 1106 1107 # Adds the atom's position to the new row 1108 self.particle[pposition].append(str(x)) 1109 self.surface[sposition].append(str(x)) 1110 self.particle[pposition].append(str(y)) 1111 self.surface[sposition].append(str(y)) 1112 self.particle[pposition].append(str(z)) 1113 self.surface[sposition].append(str(z)) 1114 1115 # Adds the atom's image flags to the new row 1116 self.particle[pposition].append('0') 1117 self.surface[sposition].append('0') 1118 self.particle[pposition].append('0') 1119 self.surface[sposition].append('0') 1120 self.particle[pposition].append('0') 1121 self.surface[sposition].append('0') 1122 1123 1124 def createxyz(self,file,data,routine='mass', values=None): 1125 """This shows the particle surface. To show the particle surface after bonding has occurred, 1126 you will need to extract the particle than reinsert the particle into the class and use createxyz. 1127 Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user. 1128 The mass version is assessed by setting the routine to 'mass' which is the default method. 1129 The other version is assessed by setting the routine to 'atomtype'. 1130 The other version takes values which is a list containing the value the user wants to assign those atomtypes to. 1131 The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1. 1132 This makes it easier to find the atomtype and assign the value from the list 1133 All atom data is assumed to have image flags in the data.""" 1134 f=open(file,'w') 1135 f.write('{0}\n'.format(len(self.surface))) 1136 f.write('atoms\n') 1137 if routine=='mass': 1138 for line in self.surface: 1139 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 1140 self.atomtype=='molecular': 1141 type=int(line[2])-1 #3rd position 1142 else: 1143 type=int(line[1])-1 #2nd position 1144 mass=data.masses[type][1] 1145 if mass=='12.0107': elementnum=6 1146 elif mass=='15.9994': elementnum=8 1147 elif mass=='26.981539':elementnum=13 1148 else: 1149 print 'no matching mass value. The method will exit' 1150 return 1151 l=len(line)-1 1152 if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\ 1153 self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\ 1154 self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\ 1155 self.atomtype=='peri': 1156 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 1157 elif self.atomtype=='dipole': 1158 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6])) 1159 elif self.atomtype=='ellipsoid': 1160 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7])) 1161 elif self.atomtype=='hybrid': 1162 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4])) 1163 f.close() 1164 elif routine=='atomtype': 1165 for line in self.surface: 1166 if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\ 1167 self.atomtype=='molecular': 1168 type=int(line[2])-1 #3rd position 1169 else: 1170 type=int(line[1])-1 #2nd position 1171 elementnum=values[type] 1172 l=len(line)-1 1173 if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\ 1174 self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\ 1175 self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\ 1176 self.atomtype=='peri': 1177 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 1178 elif self.atomtype=='dipole': 1179 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6])) 1180 elif self.atomtype=='ellipsoid': 1181 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7])) 1182 elif self.atomtype=='hybrid': 1183 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4])) 1184 f.close() 1185 1186def molecules(data,init,final,processors, method='all'): 1187 """ There are two ways to initialize the class Lmpsmolecule. 1188 Method controls which way is chosen. 1189 Acceptable values for Method are 'all', 'atom'. 1190 The default method is 'all'""" 1191 molecule=[] 1192 from multiprocessing import Pool 1193 p=Pool(processors) 1194 for i in range(init,final+1): 1195 molecule.append(p.apply_async(Lmpsmolecule,(i,data,method,))) 1196 for i in range(len(molecule)): 1197 molecule[i]=molecule[i].get() 1198 p.close() 1199 p.join() 1200 return molecule 1201 1202class Lmpsmolecule: #Technically should be a meta class but written as a separate class for easier coding. 1203 def __init__(self,moleculenum,data,method): 1204 """initiates lammps molecule structures 1205 and than extract the appropriate molecular structures from the base class data""" 1206# ***** Lammps Molecule Structures 1207 self.keywords=['Atoms']#include atoms here since this will be in every instance of this class 1208 self.atoms=[] #extract 1209 self.angles=[] #extract if keyword is there 1210 self.bonds=[] #extract if keyword is there 1211 self.dihedrals=[] #extract if keyword is there 1212 self.impropers=[] #extract if keyword is there 1213 self.velocities=[] #extract if keyword is there 1214# ***** Molecule number 1215 self.moleculenum=moleculenum 1216# ***** Extracting the moleculue's atom information from data 1217 self.extract(data,method) 1218# ***** If methods is 'all' then self.extract will extract all of the molecule's information from data 1219 1220 def extract(self,data,method): 1221 """uses base class data to extract the molecules atoms 1222 and any other molecule data from molecule moleculenum. This extraction is done in a two step process. 1223 Step 1: extract data.atoms with moleculenum and renumber self.atoms beginning from 1 1224 store extracted atom numbers from data.atoms in a list called temp 1225 Step 2: pass temp and data to method changeatomnum 1226 which extracts the rest of the Lammps Molecule structures from data 1227 and changes the atomnumbers in those structures to match the ones already in self.atoms""" 1228 #checking to make sure atomtype is a valid molecule 1229 #if not a valid molecule print error and exit method 1230 if data.atomtype=='full' or data.atomtype=='molecular': 1231# ***** Sets the molecule's atomtype 1232 self.atomtype=data.atomtype 1233 else: 1234 print "not a valid molecular structure" 1235 return 1236 #extract the molecule self.moleculenum from data.atom to self.atom 1237 atomnum=1 1238 temp=[] 1239 for i in range(len(data.atoms)): 1240 if int(data.atoms[i][1])==self.moleculenum: 1241 temp.append(data.atoms[i][0]) #store extracted atomnumbers into temp 1242 self.atoms.append(data.atoms[i]) #store extracted atom into self.atom 1243 self.atoms[atomnum-1][0]=str(atomnum) #change extracted atom to the correct atomnumber 1244 atomnum+=1 #update atomnumber 1245 #extract the rest of the Lammps Molecule structures from data using changeatomnum 1246 if method=='atom': 1247 return 1248 else: 1249 self.changeatomnum(temp,data) 1250 1251 def changeatomnum(self,temp,data=' '): 1252 """changes the atomnumbers in Lmpsmolecule structures angles, bonds, dihedrals, impropers, 1253 and velocities in a three step process. 1254 If data is defined extract the rest of the keywords used for Lmpsmolecule. 1255 Step 1A:If data is defined copy from data one of the above structures. 1256 If data is not defined skip this step. 1257 Step 1B:If data is defined remove the rows of copy which do not contain any values from temp 1258 Step 2: If data is defined extract from copy to an above Lmpsmolecule structure and remove the extracted rows 1259 If data is not defined skip this step. 1260 Step 3: Use temp and self.atoms to convert the atomnumbers in Lmpsmolecule structures 1261 to the correct values. 1262 temp and self.atoms line up, so temp[i] and self.atoms[i][0] 1263 correspond to the current Lmpsmolecule structures atom numbers 1264 and correct values 1265 will use booleanarray class to ensure that lmpsmolecule's structure values are altered only once.""" 1266 #Step 1 and Step 2 1267 if data!=' ': 1268 #extract molecular keywords from data.keywords 1269 for keywords in data.keywords: 1270 if keywords=='Angles' or keywords=='Bonds' or keywords=='Dihedrals' or\ 1271 keywords=='Impropers' or keywords=='Velocities': 1272 self.keywords.append(keywords) 1273 for keywords in self.keywords: 1274 if keywords!='Atoms': print 'extracting the data from', keywords 1275 if keywords=='Angles': 1276 copy=self.copy(data.angles) #Step 1A 1277 dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep 1278 for j in range(len(copy)): #Step 1B 1279 dnr.append(0) #adds 0 to all list elements of dnr 1280 for item in temp: #finds copied structure that has temp values 1281 for j in range(len(copy)): 1282 for i in range(2,len(copy[j])): 1283 if copy[j][i]==item: 1284 dnr[j]=1 #changes jth list element of dnr to 1 1285 break 1286 remove=[] 1287 for j in range(len(dnr)): # finds dnr values that are still 0 1288 if dnr[j]==0: 1289 remove.append(j) #and appends their index value to the list remove 1290 copy=self.deleterows(copy,remove) #removes all unneeded rows from copy 1291 structnum=1 1292 print 'the length of data is', len(copy) 1293 for item in temp: #Step 2 1294 found=[] 1295 for j in range(len(copy)): 1296 for i in range(2,len(copy[j])): 1297 if copy[j][i]==item: 1298 found.append(j) 1299 self.angles.append(copy[j]) 1300 l=len(self.angles)-1 1301 self.angles[l][0]=str(structnum) 1302 structnum+=1 1303 break 1304 # print 'the item is', item 1305 # print 'found is', found 1306 copy=self.deleterows(copy,found) 1307 elif keywords=='Bonds': 1308 copy=self.copy(data.bonds) #Step 1A 1309 dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep 1310 for j in range(len(copy)): #Step 1B 1311 dnr.append(0) #adds 0 to all list elements of dnr 1312 for item in temp: #finds copied structure that has temp values 1313 for j in range(len(copy)): 1314 for i in range(2,len(copy[j])): 1315 if copy[j][i]==item: 1316 dnr[j]=1 #changes jth list element of dnr to 1 1317 break 1318 remove=[] 1319 for j in range(len(dnr)): # finds dnr values that are still 0 1320 if dnr[j]==0: 1321 remove.append(j) #and appends their index value to the list remove 1322 copy=self.deleterows(copy,remove) #removes all unneeded rows from copy 1323 structnum=1 1324 print 'the length of data is', len(copy) 1325 for item in temp: #Step 2 1326 found=[] 1327 for j in range(len(copy)): 1328 for i in range(2,len(copy[j])): 1329 if copy[j][i]==item: 1330 found.append(j) 1331 self.bonds.append(copy[j]) 1332 l=len(self.bonds)-1 1333 self.bonds[l][0]=str(structnum) 1334 structnum+=1 1335 break 1336 # print 'the item is', item 1337 # print 'found is', found 1338 copy=self.deleterows(copy,found) 1339 elif keywords=='Dihedrals': 1340 copy=self.copy(data.dihedrals) #Step 1A 1341 dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep 1342 for j in range(len(copy)): #Step 1B 1343 dnr.append(0) #adds 0 to all list elements of dnr 1344 for item in temp: #finds copied structure that has temp values 1345 for j in range(len(copy)): 1346 for i in range(2,len(copy[j])): 1347 if copy[j][i]==item: 1348 dnr[j]=1 #changes jth list element of dnr to 1 1349 break 1350 remove=[] 1351 for j in range(len(dnr)): # finds dnr values that are still 0 1352 if dnr[j]==0: 1353 remove.append(j) #and appends their index value to the list remove 1354 copy=self.deleterows(copy,remove) #removes all unneeded rows from copy 1355 structnum=1 1356 print 'the length of data is', len(copy) 1357 structnum=1 1358 for item in temp: #Step 2 1359 found=[] 1360 for j in range(len(copy)): 1361 for i in range(2,len(copy[j])): 1362 if copy[j][i]==item: 1363 found.append(j) 1364 self.dihedrals.append(copy[j]) 1365 l=len(self.dihedrals)-1 1366 self.dihedrals[l][0]=str(structnum) 1367 structnum+=1 1368 break 1369 # print 'the item is', item 1370 # print 'found is', found 1371 copy=self.deleterows(copy,found) 1372 elif keywords=='Impropers': 1373 copy=self.copy(data.impropers) #Step 1B 1374 dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep 1375 for j in range(len(copy)): #Step 1B 1376 dnr.append(0) #adds 0 to all list elements of dnr 1377 for item in temp: #finds copied structure that has temp values 1378 for j in range(len(copy)): 1379 for i in range(2,len(copy[j])): 1380 if copy[j][i]==item: 1381 dnr[j]=1 #changes jth list element of dnr to 1 1382 break 1383 remove=[] 1384 for j in range(len(dnr)): # finds dnr values that are still 0 1385 if dnr[j]==0: 1386 remove.append(j) #and appends their index value to the list remove 1387 copy=self.deleterows(copy,remove) #removes all unneeded rows from copy 1388 structnum=1 1389 print 'the length of data is', len(copy) 1390 for item in temp: #Step 2 1391 found=[] 1392 for j in range(len(copy)): 1393 for i in range(2,len(copy[j])): 1394 if copy[j][i]==item: 1395 found.append(j) 1396 self.impropers.append(copy[j]) 1397 l=len(self.impropers)-1 1398 self.impropers[l][0]=str(structnum) 1399 structnum+=1 1400 break 1401 # print 'the item is', item 1402 # print 'found is', found 1403 copy=self.deleterows(copy,found) 1404 elif keywords=='Velocities': 1405 copy=self.copy(data.velocities) #Step 1 1406 dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep 1407 for j in range(len(copy)): #Step 1B 1408 dnr.append(0) #adds 0 to all list elements of dnr 1409 for item in temp: #finds copied structure that has temp values 1410 for j in range(len(copy)): 1411 for i in range(1): 1412 if copy[j][i]==item: 1413 dnr[j]=1 #changes jth list element of dnr to 1 1414 break 1415 remove=[] 1416 for j in range(len(dnr)): # finds dnr values that are still 0 1417 if dnr[j]==0: 1418 remove.append(j) #and appends their index value to the list remove 1419 copy=self.deleterows(copy,remove) #removes all unneeded rows from copy 1420 structnum=1 1421 print 'the length of data is', len(copy) 1422 for item in temp: #Step 2 1423 found=[] 1424 for j in range(len(copy)): 1425 for i in range(1): 1426 if copy[j][i]==item: 1427 found.append(j) 1428 self.velocities.append(copy[j]) 1429 l=len(self.velocities)-1 1430 self.velocities[l][0]=str(structnum) 1431 structnum+=1 1432 break 1433 # print 'the item is', item 1434 # print 'found is', found 1435 copy=self.deleterows(copy,found) 1436 1437 #Step 3 1438 for keywords in self.keywords: 1439 if keywords!='Atoms': print 'altering data structure values for', keywords 1440 if keywords=='Angles': 1441 copy=self.copy(self.angles) 1442 array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values 1443 for i in range(len(temp)): 1444 if temp[i]==self.atoms[i][0]: continue 1445 for j in range(len(copy)): 1446 for k in range(2,len(copy[j])): 1447 if copy[j][k]==temp[i]: 1448 if array.getelement(j,k): #if the boolean array is true 1449 self.angles[j][k]=self.atoms[i][0] 1450 array.setelement(j,k,False) 1451 break 1452 elif keywords=='Bonds': 1453 copy=self.copy(self.bonds) 1454 array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values 1455 for i in range(len(temp)): 1456 if temp[i]==self.atoms[i][0]: continue 1457 for j in range(len(copy)): 1458 for k in range(2,len(copy[j])): 1459 if copy[j][k]==temp[i]: 1460 if array.getelement(j,k): #if the boolean array is true 1461 self.bonds[j][k]=self.atoms[i][0] 1462 array.setelement(j,k,False) 1463 break 1464 elif keywords=='Dihedrals': 1465 copy=self.copy(self.dihedrals) 1466 array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values 1467 for i in range(len(temp)): 1468 if temp[i]==self.atoms[i][0]: continue 1469 for j in range(len(copy)): 1470 for k in range(2,len(copy[j])): 1471 if copy[j][k]==temp[i]: 1472 if array.getelement(j,k): #if the boolean array is true 1473 self.dihedrals[j][k]=self.atoms[i][0] 1474 array.setelement(j,k,False) 1475 break 1476 elif keywords=='Impropers': 1477 copy=self.copy(self.impropers) 1478 array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values 1479 for i in range(len(temp)): 1480 if temp[i]==self.atoms[i][0]: continue 1481 for j in range(len(copy)): 1482 for k in range(2,len(copy[j])): 1483 if copy[j][k]==temp[i]: 1484 if array.getelement(j,k): #if the boolean array is true 1485 self.impropers[j][k]=self.atoms[i][0] 1486 array.setelement(j,k,False) 1487 break 1488 elif keywords=='Velocities': 1489 copy=self.copy(self.velocities) 1490 array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values 1491 for i in range(len(temp)): 1492 if temp[i]==self.atoms[i][0]: continue 1493 for j in range(len(copy)): 1494 for k in range(1): 1495 if copy[j][k]==temp[i]: 1496 if array.getelement(j,k): #if the boolean array is true 1497 self.velocities[j][k]=self.atoms[i][0] 1498 array.setelement(j,k,False) 1499 break 1500 1501 def copy(self,structure): 1502 """copies structure to same and returns same.""" 1503 same=[] 1504 for row in structure: 1505 same.append(row) 1506 return same 1507 1508 def deleterows(self,structure,rows): # run through this {structure is list of strings and rows is list 1509 #of numbers} 1510 """delete rows in a structure and shifts the structure up 1511 rows must be in increasing order for this algorithm to work correctly""" 1512 new=[] 1513 #multiple copying of b to the rows being replaced. 1514 if rows==[]: 1515 for line in structure: 1516 new.append(line) #if no rows need replacing copy structure 1517 return new 1518 for i in range(rows[0]): 1519 new.append(structure[i]) #copy structure to new until the first replaced row 1520 count=0 1521 for i in range(rows[0],len(structure)-len(rows)):# to replace rows and shift undeleted rows over 1522 for j in range(i+1+count,len(structure)): 1523 for val in rows: 1524 if val==j: 1525 count+=1 1526 break 1527 if val==j:continue 1528 else: 1529 new.append(structure[j]) 1530 break 1531 return new 1532 1533 def modifyatom(self,atomnumber,column,newvalue): 1534 """modifies self.atom[atomnumber-1][column] to newvalue. 1535 *note: newvalue needs to be a string 1536 if column is 0 than changeatomnum method might need runnning for all atoms 1537 that have had their atomnumber(column=0) modified. 1538 changeatomnum method is not ran in this method when column=0""" 1539 self.atoms[atomnumber-1][column]=newvalue 1540 1541 1542 def deleteatoms(self,atomnumbers,atomid): 1543 """Algorithm to find all atoms bonded to atomnumbers in the direction of atomid. 1544 Atomid can be a list or a single integer. The single integer corresponds to atom's atomtype value. 1545 The list corresponds to the atom's atomid values. The only difference between both cases is in the top portion of code. 1546 When all bonded atoms have been found; ie: (the modified return values from findbonds yields an empty list); 1547 delete those bonded atoms and all molecule structures which contain those atoms. 1548 Convert the list of bonded atoms into rows and delete those rows from molecule.atoms""" 1549 print 'finding atoms to delete' 1550 bondedatoms=[] 1551 # 1st iteration 1552 nextatoms=self.findbonds(atomnumbers) 1553 try: # tests whether atomid is a list or a single integer 1554 atomid[0] 1555 except TypeError: 1556 #need to remove atoms from nextatoms which dont have the proper atom id 1557 testflag=True 1558 i=0 1559 while testflag: 1560 if i>len(nextatoms)-1: break #For the case where i becomes greater than the len of nextatoms break loop 1561 if int(self.atoms[nextatoms[i]-1][2])!=atomid:#uses the atom id for the row in atoms and than checks atom id 1562 del nextatoms[i] #delets the atom at i 1563 i-=1 #and than decreases i by 1 so next atom in the list will line up when i is increased 1564 i+=1 #increase i by 1 1565 else: 1566 #need to remove atoms from nextatoms which dont have the proper atom id 1567 testflag=True 1568 i=0 1569 while testflag: 1570 if i>len(nextatoms)-1: break #For the case where i becomes greater than the len of nextatoms break loop 1571 keep=False 1572 for id in atomid: 1573 if int(self.atoms[nextatoms[i]-1][0])==id: #checking if atomid is in next atom 1574 keep=True #keep this atom 1575 break 1576 if not keep: 1577 del nextatoms[i] #delets the atom at i 1578 i-=1 #and than decreases i by 1 so next atom in the list will line up when i is increased 1579 i+=1 #increase i by 1 1580 1581 #append next atoms into bondedatoms 1582 #copy next atoms into prevatoms 1583 prevatoms=[] 1584 for atom in nextatoms: 1585 bondedatoms.append(atom) 1586 prevatoms.append(atom) 1587 1588 1589 #2nd iteration 1590 if prevatoms==[]: 1591 print 'no bonds were found in first iteration that had atomid criteria' 1592 return 1593 1594 nextatoms=self.findbonds(prevatoms) 1595 #need to remove atoms from nextatoms which are in atomnumbers 1596 for atom in atomnumbers: 1597 for i in range(len(nextatoms)): 1598 if nextatoms[i]==atom: #checking if atom is in next atom 1599 del nextatoms[i] #delete the atom at i 1600 break 1601 if nextatoms==[]: break #all bonds from find bonds have already been added to bondedatoms 1602 1603 #append next atoms into bondedatoms 1604 #copy next atoms into prevatoms 1605 prevatoms=[] 1606 for atom in nextatoms: 1607 bondedatoms.append(atom) 1608 prevatoms.append(atom) 1609 1610 #iterative process for finding the rest of the atoms bonded to atomnumbers in the direction of atomid. 1611 while prevatoms!=[]: 1612 nextatoms=self.findbonds(prevatoms) 1613 #need to remove atoms from nextatoms which are in the prevatoms 1614 for atom in bondedatoms: 1615 for i in range(len(nextatoms)): 1616 if nextatoms[i]==atom: #checking if atom is in next atom 1617 del nextatoms[i] #delete the atom at i 1618 break 1619 if nextatoms==[]: break #all bonds from find bonds have already been added to bondedatoms 1620 1621 #append next atoms into bondedatoms 1622 #copy next atoms into prevatoms 1623 prevatoms=[] 1624 for atom in nextatoms: 1625 bondedatoms.append(atom) 1626 prevatoms.append(atom) 1627 1628 print 'the atoms to delete are', bondedatoms 1629 print 'deleting atoms from structures' 1630 #delete bonded atoms from the molecule's structures except the atom structure 1631 for i in range(1,len(self.keywords)): #goes through all keywords except the atom keyword 1632 print self.keywords[i] 1633 if self.keywords[i]=='Angles': 1634 rows=self.findatomnumbers(bondedatoms,self.angles,False) 1635 rows=self.listorder(rows) #to order the rows in increasing order 1636 self.angles=self.deleterows(self.angles,rows) #requires that rows be in increasing order 1637 elif self.keywords[i]=='Bonds': 1638 rows=self.findatomnumbers(bondedatoms,self.bonds,False) 1639 rows=self.listorder(rows) #to order the rows in increasing order 1640 self.bonds=self.deleterows(self.bonds,rows)#requires that rows be in increasing order 1641 elif self.keywords[i]=='Dihedrals': 1642 rows=self.findatomnumbers(bondedatoms,self.dihedrals,False) 1643 rows=self.listorder(rows) #to order the rows in increasing order 1644 self.dihedrals=self.deleterows(self.dihedrals,rows) #requires that rows be in increasing order 1645 elif self.keywords[i]=='Impropers': 1646 rows=self.findatomnumbers(bondedatoms,self.impropers,False) 1647 rows=self.listorder(rows) #to order the rows in increasing order 1648 self.impropers=self.deleterows(self.impropers,rows) #requires that rows be in increasing order 1649 elif self.keywords[i]=='Velocities': 1650 rows=self.findatomnumbers(bondedatoms,self.velocities,True) 1651 rows=self.listorder(rows) #to order the rows in increasing order 1652 self.velocities=self.deleterows(self.velocities,rows) #requires that rows be in increasing order 1653 1654 print 'Atoms' 1655 #convert bondedatoms from atom numbers to row numbers 1656 for i in range(len(bondedatoms)): 1657 bondedatoms[i]-=1 1658 bondedatoms=self.listorder(bondedatoms) #to order the row numbers in increasing order 1659 #delete bonded atoms (row numbers) from the atom structure 1660 self.atoms=self.deleterows(self.atoms,bondedatoms) #requires that row numbers be in increasing order 1661 1662 def findatomnumbers(self,atomnumbers,structure,vflag): #need to read through this algorithm 1663 """Algorithm to find atomnumbers in a molecule structure except. 1664 Atoms structure is not handled in here. 1665 Returns a list of the rows in which the atomnumbers are contained in the molecule structure""" 1666 rows=[] 1667 if vflag: #for handling velocity structure 1668 for atom in atomnumbers: 1669 for i in range(len(structure)): 1670 if int(structure[i][0])==atom: 1671 #duplicate rows for vflag=True are not possible. 1672 rows.append(i) 1673 break 1674 else: #for handling all other structures except atoms 1675 for atom in atomnumbers: 1676 for i in range(len(structure)): 1677 for j in range(2,len(structure[i])): 1678 if int(structure[i][j])==atom: 1679 #need to make sure duplicate value of rows are not being added 1680 if rows==[]: 1681 rows.append(i) #appends the row number 1682 else: 1683 #checking for duplicate values of rows 1684 duplicateflag=0 1685 for k in range(len(rows)): 1686 if rows[k]==i: 1687 duplicateflag=1 1688 break 1689 if duplicateflag==0: #if no duplicates adds bond number 1690 rows.append(i) #appends the row number 1691 break 1692 print 'the finished row is', rows 1693 return rows 1694 1695 def listorder(self,struct): 1696 """Takes struct and organizes the list from least to greatest. 1697 If the list is already ordered this algorithm will do nothing.""" 1698 if len(struct)==1: return struct #with the length at 1; there is only one element and therefore nothing to order 1699 for i in range(1,len(struct)): #when i=0, struct will not change; therefore its skipped 1700 copy=struct[i] 1701 for j in range(i-1,-1,-1): 1702 struct[j+1]=struct[j] 1703 if copy >struct[j]: 1704 struct[j+1]=copy 1705 break 1706 elif j==0: 1707 struct[j]=copy 1708 #dont need a break here because this is the last j value in the for loop 1709 print 'the organized row is', struct 1710 return struct 1711 1712 1713 1714 def findbonds(self,atomnumbers): 1715 """Algorithm to find all atoms bonded to atomnumbers. 1716 Returns a list of atomnumbers""" 1717 #finds the bonds in which the atomnumbers are located 1718 bondids=[] 1719 for i in range(len(atomnumbers)): 1720 for j in range(len(self.bonds)): 1721 for k in range(2,len(self.bonds[j])): 1722 if int(self.bonds[j][k])==atomnumbers[i]: 1723 if bondids==[]: 1724 bondids.append(int(self.bonds[j][0])) #appends the bond number 1725 else: 1726 #checking for duplicates of bondids 1727 duplicateflag=0 1728 for l in range(len(bondids)): 1729 if bondids[l]==int(self.bonds[j][0]): 1730 duplicateflag=1 1731 break 1732 if duplicateflag==0: #if no duplicates adds bond number 1733 bondids.append(int(self.bonds[j][0])) 1734 break 1735 1736 #Using bondids find the atoms bonded to atomnumbers 1737 bondedatoms=[] 1738 from math import fabs 1739 for id in bondids: 1740 for atoms in atomnumbers: 1741 found=False 1742 for i in range(2,len(self.bonds[id-1])): 1743 if int(self.bonds[id-1][i])==atoms: 1744 j=int(fabs(i-5)) #switches the index from the atomnumber location to the other location 1745 bondedatoms.append(int(self.bonds[id-1][j])) #appends the atomnuber at the other location 1746 found=True 1747 break 1748 if found==True: break 1749 return bondedatoms 1750 1751 def findparticlebondingpoints(self,particle,atomid,cutoffdistance,bondnumber): 1752 """Particle is a particlesurfaceobject, atomid is an int. 1753 Finds atoms with atomid in the molecule which are less than the cutoffdistance from atoms in the particle. 1754 The found atoms and particles locations are stored in a 3 dimensional list called possiblebonds 1755 The first index corresponds with the molecule's atom 1756 The second index correspons with the molecule/particle combination 1757 The third index corresponds with whether the value is the molecule or the particle 1758 Always bonds the two ends of the molecule that meet cutoff requirement. 1759 All other possible bonds are randomly chosen until the required bondnumbers are met. 1760 After every bond is chosen, the particle object's boolean list is updated, and possiblebonds is updated. 1761 The update to possiblebonds involves removing the row from which the bonded molecule's atom is located 1762 and also removing the particle atom and it's corresponding bonded atom from other rows of possiblebonds. 1763 The final bonds are all stored in self.bonding as a 2 dimensional list""" 1764 possiblebonds=[] 1765 for i in range(len(self.atoms)): 1766 if int(self.atoms[i][2])==atomid: 1767 row=[] 1768 for j in range(len(particle.surface)): 1769 #assumes molecule and particle have same atomtype if not true than atomdistance will give bad value 1770 if atomdistance(self.atoms[i],particle.surface[j],self.atomtype)<=cutoffdistance: 1771 if particle.getbool(j): row.append([i,j]) #if not bonded than can form a possible bond 1772 if row!=[]:possiblebonds.append(row) #need to correct this.... 1773 1774 #initiate section which assigns bonds into bonding information 1775 bonds=0 1776 self.bondinginformation=[] #initiate new class member 1777 1778 #Checks to see if no bonds are possible [2 possible cases] 1779 if possiblebonds==[]: 1780 print 'no possible bonds can be formed' 1781 return 1782 if bondnumber==0: 1783 print 'bondnumber is 0; so, no possible bonds can form' 1784 return 1785 1786 #section which assigns a bond to the first molecule atom which can be bonded. 1787 self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,0)) 1788 bonds+=1 1789 if bonds==bondnumber: return 1790 del possiblebonds[0] #deletes possible bonds to the first molecule which can be bonded 1791 l=len(self.bondinginformation)-1 1792 possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1]) #updates possiblebonds 1793 #by removing any bonds which contain the newly bonded particle. 1794 1795 if possiblebonds==[]:return 1796 #section which finds the last molecule atom which can be bonded 1797 l=len(possiblebonds)-1 1798 while possiblebonds[l]==[]:# to find the last molecule atom which can be bonded 1799 if possiblebonds==[]:return 1800 del possiblebonds[l] #since there are no more possible bonds in this location delete 1801 l-=1 #go to next possible spot where the last molecule atom could be bonded 1802 1803 #section which assigns a bond to the last molecule atom which can be bonded. 1804 self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,l)) 1805 bonds+=1 1806 if bonds==bondnumber: return 1807 del possiblebonds[l] #deletes possible bonds to the last molecule which can be bonded 1808 l=len(self.bondinginformation)-1 1809 possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1]) 1810 1811 if bondnumber-bonds>=len(possiblebonds): #the rest of the bonds are assigned in order 1812 while possiblebonds!=[]: 1813 if possiblebonds[0]==[]: #if row of possible bonds is empty then delete the row 1814 del possiblebonds[0] 1815 else: #else find a bond in the row of possible bonds 1816 self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,0)) 1817 #dont need to update bonds since possiblebonds will become empty at the same time or before bonds 1818 #is equal to bondnumber 1819 del possiblebonds[0] 1820 l=len(self.bondinginformation)-1 1821 possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1]) 1822 return 1823 else: #the rest of the bonds are assigned randomly 1824 from random import randint #use to randomly choose an index to bond. 1825 while bonds<bondnumber: 1826 if possiblebonds==[]:break #exits while loop when possiblebonds has become an empty set. 1827 #this ensures that in the case there are not egnough viable bonds from possiblebonds to 1828 #reach the bondnumber. Than, the while loop will not become infinite. 1829 l=len(possiblebonds)-1 1830 i=randint(0,l) 1831 if possiblebonds[i]==[]: #if row of possible bonds is empty then delete the row 1832 del possiblebonds[i] 1833 else: #else find a bond in the row of possible bonds 1834 self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,i)) 1835 bonds+=1 1836 del possiblebonds[i] 1837 l=len(self.bondinginformation)-1 1838 possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1]) 1839 return 1840 1841 def particlebondinginfo(self,particle,possiblebonds,index1): 1842 """Takes list of possiblebonds and assigns a bond from possiblebonds[index1]. 1843 Then assigns false to particle.setbool at the bonding location 1844 to ensure no more bonds can form with that particle. 1845 Returns the assigned bond.""" 1846 from random import randint #use to randomly choose 0 and 1 where 1 is keep. 1847 for i in range(len(possiblebonds[index1])): 1848 if i==len(possiblebonds[index1])-1: #automattically bonds under this condition 1849 break 1850 else: #bond has random chance of forming 1851 if randint(0,1)==1: #bond forms 1852 break 1853 else: #no bond forms 1854 continue 1855 particle.setbool(possiblebonds[index1][i][1],False)#makes sure no more bonds can form 1856 return possiblebonds[index1][i] 1857 1858 def updatepossiblebonds(self,possiblebonds,particlenum): 1859 """finds the particle number in the remaining possible bonds than deletes that bonding information. 1860 This algorithm assumes that particlenum can exist only once in each row of possiblebonds.""" 1861 for i in range(len(possiblebonds)): 1862 for j in range(len(possiblebonds[i])): 1863 if possiblebonds[i][j][1]==particlenum: 1864 del possiblebonds[i][j] 1865 break #go to next row of possiblebonds 1866 return possiblebonds 1867 1868 def bondtoparticle(self,particle,atomid,newid,newcharge): 1869 """Bonds molecule to particle. Particle is a particlesurface object. 1870 Moves the molecule's atoms bonded to the particle to the particle's position. 1871 Alters the atomtype of the molecule's atoms bonded to the particle to newid. 1872 Alters the charge of the molecule's atoms bonded to the particle to newcharge 1873 Because bonds have formed between the molecule's atoms and the particle's atoms, 1874 atoms with atomid on the molecule need to be removed otherwise the molecule's atoms will have to many bonds. 1875 self.deleteatoms will take care of deleting the atoms atomid and 1876 atoms down the polymer chain in the direction away from the molecule/particle bond. 1877 Now remove the bonded particle atoms from the particle surface becaused the bonded molecule atoms 1878 have replaced those particle atoms and their bonds with the rest of the particle. 1879 Newid must be a string. Atomid can now be a list or an integer. 1880 The list is a list of atom's id values and the single integer is atom's atomtype""" 1881 #moves the molecule's atoms bonded to the particle to the particle's position 1882 #need to do this step first before the atom id's and row indexes get out of sync 1883 #which will occur after atoms get deleted. 1884 print 'modifying the bonded molecules atom information' 1885 for i in range(len(self.bondinginformation)): 1886 #patom=particle.getsurfatom(self.bondinginformation[i][1]) 1887 #if self.atomtype=='molecular':#3-5->x,y,z[molecular] 1888 # for j in range(3,6): 1889 # self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here 1890 #elif self.atomtype=='full':#4-6->x,y,z[full] 1891 # for j in range(4,7): 1892 # self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here 1893 #alters the atomtype to newid of the molecule's atoms bonded to the particle. 1894 self.modifyatom(self.bondinginformation[i][0]+1,2,newid) #uses the atomid here 1895 if self.atomtype=='full': 1896 # Alters the charge of the molecule's atoms bonded to the particle to newcharge 1897 self.modifyatom(self.bondinginformation[i][0]+1,3,newcharge) 1898 1899 #This is a separate loop so atom id's and row indexes for previous steps wont get out of sync 1900 #create atomnumbers to begin deletion process. 1901 atomnumbers=[] 1902 for i in range(len(self.bondinginformation)): 1903 atomnumbers.append(self.bondinginformation[i][0]+1)#using actual atomnumbers rather than the row index 1904 1905 #Call algorithm to find all atoms bonded to atomnumbers in the direction of atomid. 1906 #Than delete those atoms and all molecule structures which contain those atoms. 1907 self.deleteatoms(atomnumbers,atomid) 1908 1909 print 'beginning deletion process of the surface atoms for which the molecule atoms have replaced' 1910 #Goes through the bondinginformation and superficially removes the surfaceatom 1911 for i in range(len(self.bondinginformation)): 1912 particle.removesurfatom(self.bondinginformation[i][1]) #uses the row number 1913 #Allows the particle extract method used outside this class to remove this atom. 1914 1915 def createxyz(self,file,data,routine='mass', values=None): 1916 """Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user. 1917 The mass version is assessed by setting the routine to 'mass' which is the default method. 1918 The other version is assessed by setting the routine to 'atomtype'. 1919 The other version takes values which is a list containing the value the user wants to assign those atomtypes to. 1920 The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1. 1921 This makes it easier to find the atomtype and assign the value from the list 1922 All atom data is assumed to have image flags in the data.""" 1923 f=open(file,'w') 1924 f.write('{0}\n'.format(len(self.atoms))) 1925 f.write('atoms\n') 1926 if routine=='mass': 1927 for line in self.atoms: 1928 type=int(line[2])-1 1929 mass=data.masses[type][1] 1930 if mass=='12.0107': elementnum=6 1931 elif mass=='15.9994': elementnum=8 1932 elif mass=='26.9815':elementnum=13 1933 else: 1934 print 'no matching mass value. The method will exit' 1935 return 1936 l=len(line)-1 1937 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 1938 f.close() 1939 elif routine=='atomtype': 1940 for line in self.atoms: 1941 type=int(line[2])-1 1942 elementnum=values[type] 1943 l=len(line)-1 1944 f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3])) 1945 f.close() 1946