MolKit.amberPrmTop

30 """class to hold parameters for Amber Force Field calcuations 31 """

32 - def __init__(self, allDictList = [all_amino94_dat], ntDictList = [all_aminont94_dat], 33 ctDictList = [all_aminoct94_dat]):

34 #amber parameter reference dictionaries: 35 if len(allDictList)==0: 36 allDict = all_amino94_dat 37 else: 38 allDict = allDictList[0] 39 if type(allDict)==types.StringType: 40 allDict = self.getDictObj(allDict) 41 if len(allDictList)>1: 42 for d in allDictList: 43 if type(d)== types.StringType: 44 d = self.getDictObj(d) 45 allDict.update(d) 46 #allDict.extend(d) 47 self.allDict = allDict 48 if len(ntDictList)==0: 49 ntDict = all_aminont94_dat 50 else: 51 ntDict = ntDictList[0] 52 if type(ntDict)==types.StringType: 53 ntDict = self.getDictObj(ntDict) 54 if len(ntDictList)>1: 55 for d in ntDictList: 56 if type(d)== types.StringType: 57 d = self.getDictObj(d) 58 ntDict.update(d) 59 #ntDict.extend(d) 60 self.ntDict = ntDict 61 if len(ctDictList)==0: 62 ctDict = all_aminoct94_dat 63 else: 64 ctDict = ctDictList[0] 65 if type(ctDict)==types.StringType: 66 ctDict = self.getDictObj(ctDict) 67 if len(ctDictList)>1: 68 for d in ctDictList: 69 if type(d)== types.StringType: 70 d = self.getDictObj(d) 71 ctDict.update(d) 72 #ctDict.extend(d) 73 self.ctDict = ctDict 74 #formatD is used for write method 75 formatD = {} 76 for k in ['Iac', 'Iblo', 'Cno', 'Ipres', 'ExclAt']: 77 formatD[k] = ('%6d', 12, 0) 78 for k in ['Charges', 'Masses', 'Rk', 'Req', 'Tk', 'Teq',\ 79 'Pk', 'Pn', 'Phase', 'Solty', 'Cn1', 'Cn2']: 80 formatD[k] = ('%16.8E', 5, 0) 81 for k in ['AtomNames', 'ResNames', 'AtomSym', 'AtomTree']: 82 formatD[k] = ('%-4.4s', 20, 0) 83 for k in ['allHBnds', 'allBnds']: 84 formatD[k] = ('%6d', 12, 3) 85 for k in ['allHAngs', 'allAngs']: 86 formatD[k] = ('%6d', 12, 4) 87 for k in ['allHDihe', 'allDihe']: 88 formatD[k] = ('%6d', 12, 5) 89 self.formatD = formatD 90 #processAtoms results are built in this dictionary 91 self.prmDict = {}

92 93

94 - def getDictObj(self, nmstr):

95 #mod = __import__('MolKit/data/' + nmstr) 96 #dict = eval('mod.'+ nmstr) 97 mod = __import__('MolKit') 98 b = getattr(mod.data, nmstr) 99 dict = getattr(b, nmstr) 100 return dict

101 102 103 104

105 - def loadFromFile(self, filename):

106 """reads a parmtop file""" 107 self.prmDict = self.py_read(filename) 108 self.createSffCdataStruct(self.prmDict)

109 110

111 - def processAtoms(self, atoms, parmDict=None, reorder=1):

112 """finds all Amber parameters for the given set of atoms 113 parmDict is parm94_dat """ 114 115 116 if atoms: 117 self.build(atoms, parmDict, reorder) 118 self.createSffCdataStruct(self.prmDict) 119 print 'after call to createSffCdataStruct'

120 121

122 - def checkSanity(self):

123 d = self.prmDict 124 #length checks: 125 Natom = d['Natom'] 126 assert len(d['Charges']) == Natom 127 assert len(d['Masses']) == Natom 128 assert len(d['Iac']) == Natom 129 assert len(d['Iblo']) == Natom 130 assert len(d['AtomRes']) == Natom 131 assert len(d['N14pairs']) == Natom 132 assert len(d['TreeJoin']) == Natom 133 134 Nres = d['Nres'] 135 assert len(d['Ipres']) == Nres + 1 136 137 assert len(d['AtomNames']) == Natom * 4 + 81 138 assert len(d['AtomSym']) == Natom * 4 + 81 139 assert len(d['AtomTree']) == Natom * 4 + 81 140 assert len(d['ResNames']) == Nres * 4 + 81 141 142 #Ntypes is number of unique amber_types w/equiv replacement 143 Ntypes = d['Ntypes'] 144 assert len(d['Cno']) == Ntypes**2 145 assert len(d['ExclAt']) == d['Nnb'] 146 assert len(d['Cn1']) == Ntypes*(Ntypes+1)/2. 147 assert len(d['Cn2']) == Ntypes*(Ntypes+1)/2. 148 149 #Numbnd is number of bnd types 150 Numbnd = d['Numbnd'] 151 assert len(d['Rk']) == Numbnd 152 assert len(d['Req']) == Numbnd 153 #Numang is number of angle types 154 Numang = d['Numang'] 155 assert len(d['Tk']) == Numang 156 assert len(d['Teq']) == Numang 157 #Numptra is number of dihe types 158 Nptra = d['Nptra'] 159 assert len(d['Pk']) == Nptra 160 assert len(d['Pn']) == Nptra 161 assert len(d['Phase']) == Nptra 162 163 assert len(d['Solty']) == d['Natyp'] 164 165 #Nbona is number of bonds w/out H 166 Nbona = d['Nbona'] 167 assert len(d['BondAt1']) == Nbona 168 assert len(d['BondAt2']) == Nbona 169 assert len(d['BondNum']) == Nbona 170 #Nbonh is number of bonds w/ H 171 Nbonh = d['Nbonh'] 172 assert len(d['BondHAt1']) == Nbonh 173 assert len(d['BondHAt2']) == Nbonh 174 assert len(d['BondHNum']) == Nbonh 175 176 #Ntheta is number of angles w/out H 177 Ntheta = d['Ntheta'] 178 assert len(d['AngleAt1']) == Ntheta 179 assert len(d['AngleAt2']) == Ntheta 180 assert len(d['AngleAt3']) == Ntheta 181 assert len(d['AngleNum']) == Ntheta 182 183 #Ntheth is number of angles w/ H 184 Ntheth = d['Ntheth'] 185 assert len(d['AngleHAt1']) == Ntheth 186 assert len(d['AngleHAt2']) == Ntheth 187 assert len(d['AngleHAt3']) == Ntheth 188 assert len(d['AngleHNum']) == Ntheth 189 190 #Nphia is number of dihedrals w/out H 191 Nphia = d['Nphia'] 192 assert len(d['DihAt1']) == Nphia 193 assert len(d['DihAt2']) == Nphia 194 assert len(d['DihAt3']) == Nphia 195 assert len(d['DihAt4']) == Nphia 196 assert len(d['DihNum']) == Nphia 197 198 #Nphih is number of dihedrals w/ H 199 Nphih = d['Nphih'] 200 assert len(d['DihHAt1']) == Nphih 201 assert len(d['DihHAt2']) == Nphih 202 assert len(d['DihHAt3']) == Nphih 203 assert len(d['DihHAt4']) == Nphih 204 assert len(d['DihHNum']) == Nphih 205 206 ##WHAT ABOUT HB10, HB12, N14pairs, N14pairlist 207 208 #value based on length checks: 209 #all values of BondNum and BondHNum in range (1, Numbnd) 210 for v in d['BondNum']: 211 assert v >0 and v < Numbnd + 1 212 for v in d['BondHNum']: 213 assert v >0 and v < Numbnd + 1 214 #all values of AngleNum and AngleHNum in range (1, Numang) 215 for v in d['AngleNum']: 216 assert v >0 and v < Numang + 1 217 for v in d['AngleHNum']: 218 assert v >0 and v < Numang + 1 219 #all values of DihNum and DihHNum in range (1, Nptra) 220 for v in d['DihNum']: 221 assert v >0 and v < Nptra + 1 222 for v in d['DihHNum']: 223 assert v >0 and v < Nptra + 1

224 225

226 - def createSffCdataStruct(self, dict):

227 """Create a C prm data structure""" 228 print 'in createSffCdataStruct' 229 self.ambPrm = AmberParm('test1', parmdict=dict) 230 print 'after call to init'

231 232

233 - def build(self, allAtoms, parmDict, reorder):

234 235 # find out amber special residue name and 236 # order the atoms inside a residue to follow the Amber convention 237 self.residues = allAtoms.parent.uniq() 238 239 self.residues.sort() 240 241 self.fixResNamesAndOrderAtoms(reorder) 242 243 # save ordered chains 244 self.chains = self.residues.parent.uniq() 245 self.chains.sort() 246 247 # save ordered atoms 248 self.atoms = self.residues.atoms 249 250 # renumber them 251 self.atoms.number = range(1, len(allAtoms)+1) 252 print 'after call to checkRes' 253 254 self.getTopology(self.atoms, parmDict) 255 print 'after call to getTopology' 256 257 if reorder: 258 self.checkSanity() 259 print 'passed sanity check' 260 else: 261 print 'skipping sanity check' 262 return

263 264

265 - def reorderAtoms(self, res, atList):

266 ats = [] 267 rlen = len(res.atoms) 268 if rlen!=len(atList): 269 print "atoms missing in residue", res 270 print "expected:", atList 271 print "found :", res.atoms.name 272 for i in range(rlen): 273 a = atList[i] 274 for j in range(rlen): 275 b = res.atoms[j] 276 # DON'T rename HN atom H, HN1->H1, etc... 277 # use editCommands instead 278 #if b.name=='HN': b.name='H' 279 #elif len(b.name)==3 and b.name[:2]=='HN': 280 #b.name ='H'+b.name[2] 281 if b.name==a: 282 ats.append(b) 283 break 284 if len(ats)==len(res.atoms): 285 res.children.data = ats 286 res.atoms.data = ats

287 288

289 - def fixResNamesAndOrderAtoms(self, reorder):

290 # level list of atom names used to rename residues 291 # check is HIS is HIS, HID, HIP, HIE, etc... 292 293 residues = self.residues 294 last = len(residues)-1 295 for i in range(len(residues)): 296 residue = residues[i] 297 chNames = residue.atoms.name 298 299 amberResType = residue.type 300 301 if amberResType=='CYS': 302 returnVal = 'CYS' 303 #3/21: 304 if 'HSG' in chNames or 'HG' in chNames: 305 amberResType ='CYS' 306 elif 'HN' in chNames: 307 amberResType = 'CYM' 308 else: 309 amberResType = 'CYX' 310 elif amberResType=='LYS': 311 # THIS DOESN'T SUPPORT LYH assigned in all.in 312 returnVal = 'LYS' 313 if 'HZ1' in chNames or 'HZN1' in chNames: 314 amberResType ='LYS' 315 else: 316 amberResType ='LYN' 317 elif amberResType=='ASP': 318 returnVal = 'ASP' 319 #3/21 320 if 'HD' in chNames or 'HD2' in chNames: 321 amberResType ='ASH' 322 else: 323 amberResType ='ASP' 324 elif amberResType=='GLU': 325 returnVal = 'GLU' 326 #3/21 327 if 'HE' in chNames or 'HE2' in chNames: 328 amberResType ='GLH' 329 else: 330 amberResType ='GLU' 331 elif amberResType=='HIS': 332 returnVal = 'HIS' 333 hasHD1 = 'HD1' in chNames 334 hasHD2 = 'HD2' in chNames 335 hasHE1 = 'HE1' in chNames 336 hasHE2 = 'HE2' in chNames 337 if hasHD1 and hasHE1: 338 if hasHD2 and not hasHE2: 339 amberResType = 'HID' 340 elif hasHD2 and hasHE2: 341 amberResType = 'HIP' 342 elif (not hasHD1) and (hasHE1 and hasHD2 and hasHE2): 343 amberResType = 'HIE' 344 else: 345 print 'unknown HISTIDINE config' 346 raise ValueError 347 348 residue.amber_type = amberResType 349 if residue == residue.parent.residues[0]: 350 residue.amber_dict = self.ntDict[amberResType] 351 elif residue == residue.parent.residues[-1]: 352 residue.amber_dict = self.ctDict[amberResType] 353 else: 354 residue.amber_dict = self.allDict[amberResType] 355 356 if reorder: 357 self.reorderAtoms(residue, residue.amber_dict['atNameList'])

358 359

360 - def processChain(self, residues, parmDict):

361 #this should be called with a list of residues representing a chain 362 # NOTE: self.parmDict is parm94 which was parsed by Ruth while parmDict is 363 # MolKit.parm94_dat.py 364 dict = self.prmDict 365 #residues = self.residues 366 367 # initialize 368 atNames = '' 369 atSym = '' 370 atTree = '' 371 resname = '' 372 373 masses = dict['Masses'] 374 charges = dict['Charges'] 375 uniqList = [] 376 uniqTypes = {} # used to build list with equivalent names removed 377 atypTypes = {} # used to build list without equivalent names removed 378 allTypeList = [] # list of all types 379 380 last = len(residues)-1 381 dict['Nres'] = dict['Nres'] + last + 1 382 atres = dict['AtomRes'] 383 ipres = dict['Ipres'] 384 maxResLen = 0 385 386 for i in range(last+1): 387 res = residues[i] 388 atoms = res.atoms 389 390 nbat = len(atoms) 391 if nbat > maxResLen: maxResLen = nbat 392 ipres.append(ipres[-1]+nbat) 393 resname = resname + res.amber_type + ' ' 394 395 ad = res.amber_dict 396 pdm = parmDict.atomTypes 397 398 for a in atoms: 399 # get the amber atom type 400 name = a.name 401 atres.append(i+1) 402 atNames = atNames+'%-4s'%name 403 atD = ad[name] 404 a.amber_type = newtype = '%-2s'%atD['type'] 405 chg = a._charges['amber'] = atD['charge']*18.2223 406 charges.append(chg) 407 mas = a.mass = pdm[newtype][0] 408 masses.append(mas) 409 atTree = atTree+'%-4.4s'%atD['tree'] 410 allTypeList.append(newtype) 411 atSym = atSym+'%-4s'%newtype 412 symb = newtype[0] 413 if symb in parmDict.AtomEquiv.keys(): 414 if newtype in parmDict.AtomEquiv[symb]: 415 newsym = symb + ' ' 416 uniqTypes[symb+' '] = 0 417 a.amber_symbol = symb+' ' 418 if newsym not in uniqList: 419 uniqList.append(newsym) 420 else: 421 uniqTypes[newtype] = 0 422 a.amber_symbol = newtype 423 if newtype not in uniqList: 424 uniqList.append(newtype) 425 else: 426 uniqTypes[newtype] = 0 427 a.amber_symbol = newtype 428 if newtype not in uniqList: uniqList.append(newtype) 429 # to get uniq list of all types w/out equiv replacement 430 atypTypes[newtype] = 0 431 432 # post processing of some variable 433 dict['AtomNames'] = dict['AtomNames'] + atNames 434 dict['AtomSym'] = dict['AtomSym'] + atSym 435 dict['AtomTree'] = dict['AtomTree'] + atTree 436 dict['ResNames'] = dict['ResNames'] + resname 437 438 # save list of unique types for later use 439 ###1/10: 440 #self.uniqTypeList = uniqList 441 uL = self.uniqTypeList 442 for t in uniqList: 443 if t not in uL: 444 uL.append(t) 445 #self.uniqTypeList = uniqTypes.keys() 446 self.uniqTypeList = uL 447 448 ntypes = len(uL) 449 dict['Ntypes'] = ntypes 450 451 aL = self.atypList 452 for t in atypTypes.keys(): 453 if t not in aL: 454 aL.append(t) 455 self.atypList = aL 456 dict['Natyp'] = len(aL) 457 458 dict['Ntype2d'] = ntypes*ntypes 459 dict['Nttyp'] = ntypes * (ntypes+1)/2 460 if maxResLen > dict['Nmxrs']: 461 dict['Nmxrs'] = maxResLen 462 463 464 newtypelist = [] 465 for t in residues.atoms.amber_symbol: 466 # Iac is 1-based 467 newtypelist.append( self.uniqTypeList.index(t) + 1 ) 468 ###1/10: 469 #dict['Iac'] = newtypelist 470 dict['Iac'].extend( newtypelist)

471 472

473 - def processBonds(self, bonds, parmDict):

474 # NOTE: self,parmDict is parm94 parsed by Ruth while parmDict is 475 # MolKit.parm94_dat.py): 476 477 dict = self.prmDict 478 bat1 = dict['BondAt1'] 479 bat2 = dict['BondAt2'] 480 bnum = dict['BondNum'] 481 batH1 = dict['BondHAt1'] 482 batH2 = dict['BondHAt2'] 483 bHnum = dict['BondHNum'] 484 rk = dict['Rk'] 485 req = dict['Req'] 486 487 bndTypes = {} # used to build a unique list of bond types 488 btDict = parmDict.bondTypes #needed to check for wildcard * in type 489 490 for b in bonds: 491 a1 = b.atom1 492 #t1 = a1.amber_symbol 493 t1 = a1.amber_type 494 a2 = b.atom2 495 #t2 = a2.amber_symbol 496 t2 = a2.amber_type 497 if t1<t2: 498 newtype = '%-2.2s-%-2.2s'%(t1,t2) 499 else: 500 newtype = '%-2.2s-%-2.2s'%(t2,t1) 501 bndTypes[newtype] = 0 502 503 n1 = (a1.number-1)*3 504 n2 = (a2.number-1)*3 505 if n2<n1: tmp=n1; n1=n2; n2=tmp 506 507 if a1.element=='H' or a2.element=='H': 508 bHnum.append(newtype) 509 batH1.append(n1) 510 batH2.append(n2) 511 else: 512 bnum.append(newtype) 513 bat1.append(n1) 514 bat2.append(n2) 515 516 dict['Numbnd'] = len(bndTypes) 517 btlist = bndTypes.keys() 518 519 for bt in btlist: 520 rk.append( btDict[bt][0] ) 521 req.append( btDict[bt][1] ) 522 523 newbnum = [] 524 for b in bnum: 525 newbnum.append( btlist.index(b) + 1 ) 526 dict['BondNum'] = newbnum 527 528 newbnum = [] 529 for b in bHnum: 530 newbnum.append( btlist.index(b) + 1 ) 531 dict['BondHNum'] = newbnum 532 533 return

534 535

536 - def processAngles(self, allAtoms, parmDict):

537 dict = self.prmDict 538 aa1 = dict['AngleAt1'] 539 aa2 = dict['AngleAt2'] 540 aa3 = dict['AngleAt3'] 541 anum = dict['AngleNum'] 542 aHa1 = dict['AngleHAt1'] 543 aHa2 = dict['AngleHAt2'] 544 aHa3 = dict['AngleHAt3'] 545 aHnum = dict['AngleHNum'] 546 tk = dict['Tk'] 547 teq = dict['Teq'] 548 549 angTypes = {} 550 atdict = parmDict.bondAngles 551 552 for a1 in allAtoms: 553 t1 = a1.amber_type 554 for b in a1.bonds: 555 a2 = b.atom1 556 if id(a2)==id(a1): a2=b.atom2 557 t2 = a2.amber_type 558 for b2 in a2.bonds: 559 a3 = b2.atom1 560 if id(a3)==id(a2): a3=b2.atom2 561 if id(a3)==id(a1): continue 562 if a1.number > a3.number: continue 563 564 t3 = a3.amber_type 565 566 nn1 = n1 = (a1.number-1)*3 567 nn2 = n2 = (a2.number-1)*3 568 nn3 = n3 = (a3.number-1)*3 569 if n3<n1: 570 nn3 = n1 571 nn1 = n3 572 573 rev = 0 574 if (t1==t3 and a1.name > a3.name) or t3 < t1: 575 rev = 1 576 577 if rev: 578 newtype = '%-2.2s-%-2.2s-%-2.2s'%(t3,t2,t1) 579 else: 580 newtype = '%-2.2s-%-2.2s-%-2.2s'%(t1, t2, t3) 581 #have to check for wildcard * 582 angTypes[newtype] = 0 583 584 if a1.element=='H' or a2.element=='H' or a3.element=='H': 585 aHa1.append( nn1 ) 586 aHa2.append( nn2 ) 587 aHa3.append( nn3 ) 588 aHnum.append(newtype) 589 else: 590 aa1.append( nn1 ) 591 aa2.append( nn2 ) 592 aa3.append( nn3 ) 593 anum.append(newtype) 594 595 atlist = angTypes.keys() 596 597 torad = pi / 180.0 598 atKeys = atdict.keys() 599 for t in atlist: 600 tk.append( atdict[t][0] ) 601 teq.append( atdict[t][1]*torad ) 602 603 anewlist = [] 604 for a in anum: 605 anewlist.append( atlist.index( a ) + 1 ) 606 dict['AngleNum'] = anewlist 607 608 anewlist = [] 609 for a in aHnum: 610 anewlist.append( atlist.index( a ) + 1 ) 611 dict['AngleHNum'] = anewlist 612 613 dict['Numang'] = len(atlist) 614 dict['Ntheth'] = len(aHa1) 615 dict['Mtheta'] = len(aa1) 616 dict['Ntheta'] = len(aa1) 617 618 return

619 620

621 - def checkDiheType(self, t, t2, t3, t4, dict):

622 #zero X 623 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%(t,t2,t3,t4) 624 if dict.has_key(newtype): return newtype 625 626 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%(t4,t3,t2,t) 627 if dict.has_key(newtype): return newtype 628 629 #X 630 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t2,t3,t4) 631 if dict.has_key(newtype): return newtype 632 633 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t3,t2,t) 634 if dict.has_key(newtype): return newtype 635 636 #2X 637 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t2,t3,'X') 638 if dict.has_key(newtype): return newtype 639 640 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t3,t2,'X') 641 if dict.has_key(newtype): return newtype 642 643 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X','X',t3,t4) 644 if dict.has_key(newtype): return newtype 645 646 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X','X',t2,t) 647 if dict.has_key(newtype): return newtype 648 649 raise RuntimeError('dihedral type not in dictionary')

650 651 ## it is slower to check a list if the key is in there than to ask a 652 ## dictionanry if it has this key 653 ## 654 ## keys = dict.keys() 655 ## #zero X 656 ## newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%(t,t2,t3,t4) 657 ## if newtype in keys: 658 ## return newtype 659 ## newtype2 = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%(t4,t3,t2,t) 660 ## if newtype2 in keys: 661 ## return newtype2 662 ## #X 663 ## newtypeX = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t2,t3,t4) 664 ## if newtypeX in keys: 665 ## return newtypeX 666 ## newtype2X = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t3,t2,t) 667 ## if newtype2X in keys: 668 ## return newtype2X 669 ## #2X 670 ## newtypeX_X = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t2,t3,'X') 671 ## if newtypeX_X in keys: 672 ## return newtypeX_X 673 ## newtype2X_X = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X',t3,t2,'X') 674 ## if newtype2X_X in keys: 675 ## return newtype2X_X 676 ## newtypeXX = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X','X',t3,t4) 677 ## if newtypeXX in keys: 678 ## return newtypeXX 679 ## newtype2XX = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%('X','X',t2,t) 680 ## if newtype2XX in keys: 681 ## return newtype2XX 682 ## raise RuntimError('dihedral type not in dictionary') 683 684

685 - def processTorsions(self, allAtoms, parmDict):

686 # find torsions and also excuded atoms 687 dict = self.prmDict 688 foundDihedTypes = {} 689 ta1 = dict['DihAt1'] 690 ta2 = dict['DihAt2'] 691 ta3 = dict['DihAt3'] 692 ta4 = dict['DihAt4'] 693 tnum = dict['DihNum'] 694 taH1 = dict['DihHAt1'] 695 taH2 = dict['DihHAt2'] 696 taH3 = dict['DihHAt3'] 697 taH4 = dict['DihHAt4'] 698 tHnum = dict['DihHNum'] 699 700 nb14 = dict['N14pairs'] 701 n14list = dict['N14pairlist'] 702 703 iblo = dict['Iblo'] 704 exclAt = dict['ExclAt'] 705 706 dihedTypes = parmDict.dihedTypes 707 708 for a1 in allAtoms: 709 n14 = [] 710 excl = [] 711 t1 = a1.amber_type 712 restyp = a1.parent.type 713 if restyp in ['PRO', 'TRP', 'HID', 'HIE', 'HIP']: 714 ringlist = self.AA5rings[restyp] 715 else: 716 ringlist = None 717 718 for b in a1.bonds: 719 a2 = b.atom1 720 if id(a2)==id(a1): a2=b.atom2 721 t2 = a2.amber_type 722 if a2.number > a1.number: excl.append(a2.number) 723 for b2 in a2.bonds: 724 a3 = b2.atom1 725 if id(a3)==id(a2): a3=b2.atom2 726 if id(a3)==id(a1): continue 727 if a3.number > a1.number: excl.append(a3.number) 728 t3 = a3.amber_type 729 for b3 in a3.bonds: 730 a4 = b3.atom1 731 if id(a4)==id(a3): a4=b3.atom2 732 if id(a4)==id(a2): continue 733 if id(a4)==id(a1): continue 734 if a1.number > a4.number: continue 735 excl.append(a4.number) 736 t4 = a4.amber_type 737 738 newtype = '%-2.2s-%-2.2s-%-2.2s-%-2.2s'%(t1,t2,t3,t4) 739 dtype = self.checkDiheType(t1,t2,t3,t4,dihedTypes) 740 741 for i in range(len(dihedTypes[dtype])): 742 tname = dtype+'_'+str(i) 743 foundDihedTypes[tname] = 0 744 sign3 = 1 745 period = dihedTypes[dtype][i][3] 746 if period < 0.0: sign3= -1 747 if a4.parent==a1.parent: 748 if ringlist and a4.name in ringlist \ 749 and a1.name in ringlist: 750 sign3= -1 751 if a1.element=='H' or a2.element=='H' or \ 752 a3.element=='H' or a4.element=='H': 753 taH1.append( (a1.number-1)*3 ) 754 taH2.append( (a2.number-1)*3 ) 755 taH3.append( sign3*(a3.number-1)*3 ) 756 taH4.append( (a4.number-1)*3 ) 757 tHnum.append( tname ) 758 else: 759 ta1.append( (a1.number-1)*3 ) 760 ta2.append( (a2.number-1)*3 ) 761 ta3.append( sign3*(a3.number-1)*3 ) 762 ta4.append( (a4.number-1)*3 ) 763 tnum.append( tname ) 764 if sign3>0.0: 765 # this trick work only for 6 rings and 766 # prevents from adding 14 interactions 767 # twice between atoms in the ring 768 # PRO, TRP and HIS and cp. have to be handle 769 # separately 770 num = a4.number-1 771 if num not in n14: 772 n14.append( num ) 773 else: # make 3rd atom in torsion negative 774 ta3[-1] = -ta3[-1] 775 if len(excl): 776 # excl can contain duplicated values (pro tyr phe cycles) 777 # we also sort the values (probably only comsetics) 778 excl.sort() 779 last = excl[0] 780 uexcl = [last] 781 for i in range(1,len(excl)): 782 if excl[i]!=last: 783 last = excl[i] 784 uexcl.append(last) 785 iblo.append(len(uexcl)) 786 exclAt.extend(uexcl) 787 else: 788 iblo.append( 1 ) 789 exclAt.append( 0 ) 790 791 nb14.append(len(n14)) 792 ##!##1/28: n14.sort() 793 n14list.extend(n14) 794 795 # remember how many proper diehedrals 796 lastProper = len(tnum) 797 lastHProper = len(tHnum) 798 799 # loop over residues to add improper torsions 800 sumAts = 0 801 foundImproperDihedTypes = {} 802 for res in self.residues: 803 foundImproperDihedTypes = self.getImpropTors( 804 res, sumAts, foundImproperDihedTypes, parmDict) 805 sumAts = sumAts + len(res.atoms) 806 807 #typeDict = foundDihedTypes.copy() 808 #typeDict.update(foundImproperDihedTypes) 809 #print typeDict.keys() 810 811 dict['Nptra'] = len(foundDihedTypes) + len(foundImproperDihedTypes) 812 dict['Mphia'] = dict['Nphia'] = len(ta1) 813 dict['Nphih'] = len(taH1) 814 815 pn = dict['Pn'] 816 pk = dict['Pk'] 817 phase = dict['Phase'] 818 dtlist = foundDihedTypes.keys() 819 torad = pi/180. 820 for t in dtlist: 821 index = int(t[-1]) 822 val = dihedTypes[t[:-2]][index] # remove the '_x' 823 pk.append(val[1]/val[0]) 824 phase.append(val[2]*torad) 825 pn.append(fabs(val[3])) 826 827 dihedTypes = parmDict.improperDihed 828 dtlist1 = foundImproperDihedTypes.keys() 829 830 for t in dtlist1: 831 val = dihedTypes[t] 832 pk.append(val[0]) 833 phase.append(val[1]*torad) 834 pn.append(val[2]) 835 836 typenum = [] 837 dtlist = dtlist + dtlist1 838 for t in tnum: 839 typenum.append( dtlist.index(t) + 1 ) # types are 1-based 840 dict['DihNum'] = typenum 841 842 typenum = [] 843 for t in tHnum: 844 typenum.append( dtlist.index(t) + 1 ) # types are 1-based 845 dict['DihHNum'] = typenum 846 847 dict['Nnb'] = len(dict['ExclAt']) 848 #print len(tnum), len(dict['DihNum']) 849 850 return

851 852 853

854 - def getImpropTors(self, res, sumAts, foundDihedTypes, parmDict):

855 #eg tList:[['CA','+M','C','0'],['-M','CA','N','H']] 856 dict = self.prmDict 857 offset = sumAts * 3 858 nameList = res.atoms.name 859 typeList = res.atoms.amber_type 860 861 ta1 = dict['DihAt1'] 862 ta2 = dict['DihAt2'] 863 ta3 = dict['DihAt3'] 864 ta4 = dict['DihAt4'] 865 tnum = dict['DihNum'] 866 taH1 = dict['DihHAt1'] 867 taH2 = dict['DihHAt2'] 868 taH3 = dict['DihHAt3'] 869 taH4 = dict['DihHAt4'] 870 tHnum = dict['DihHNum'] 871 872 dihedTypes = parmDict.improperDihed 873 atNameList = res.amber_dict['atNameList'] 874 resat = res.atoms 875 for item in res.amber_dict['impropTors']: 876 atomNum = [] 877 atomType = [] 878 newTors = [] 879 offset = res.atoms[0].number 880 #use hasH to detect 'HZ2' etc 881 hasH = 0 882 for t in item: 883 if t[0]=='H': hasH = 1 884 if len(t)==2 and t[1]=='M': 885 if t[0]=='-': 886 atomType.append('C ') 887 atomNum.append(offset - 2) 888 else: 889 atomType.append('N ') 890 atomNum.append(offset + len(res.atoms) ) 891 else: 892 atIndex = atNameList.index(t) 893 atom = resat[atIndex] 894 atomType.append(atom.amber_type) 895 atomNum.append( atom.number ) 896 897 newType = self.checkDiheType(atomType[0], atomType[1], 898 atomType[2], atomType[3], 899 dihedTypes) 900 foundDihedTypes[newType] = 0 901 902 if hasH: 903 taH1.append( (atomNum[0]-1)*3 ) 904 taH2.append( (atomNum[1]-1)*3 ) 905 taH3.append(-(atomNum[2]-1)*3 ) 906 taH4.append(-(atomNum[3]-1)*3 ) 907 tHnum.append(newType) 908 else: 909 ta1.append( (atomNum[0]-1)*3 ) 910 ta2.append( (atomNum[1]-1)*3 ) 911 ta3.append(-(atomNum[2]-1)*3 ) 912 ta4.append(-(atomNum[3]-1)*3 ) 913 tnum.append(newType) 914 915 return foundDihedTypes

916 917

918 - def getTopology(self, allAtoms, parmDict):

919 920 dict = self.prmDict 921 dict['ititl'] = allAtoms.top.uniq()[0].name + '.prmtop\n' 922 923 natom = dict['Natom'] = len(allAtoms) 924 dict['Nat3'] = natom * 3 925 926 dict['AtomNames'] = '' 927 dict['AtomSym'] = '' 928 dict['AtomTree'] = '' 929 dict['Ntypes'] = 0 930 dict['Natyp'] = 0 931 dict['Ntype2d'] = 0 932 dict['Nttyp'] = 0 933 dict['Masses'] = [] 934 dict['Charges'] = [] 935 dict['Nres'] = 0 936 dict['AtomRes'] = [] 937 dict['ResNames'] = '' 938 dict['Ipres'] = [1] 939 dict['Nmxrs'] = 0 940 ###1/10: 941 dict['Iac'] = [] 942 self.uniqTypeList = [] 943 #used for construction of Natyp 944 self.atypList = [] 945 946 # fill get all arrays that are of len natom 947 # we have to call for each chain 948 for ch in self.chains: 949 self.processChain( ch.residues, parmDict) 950 951 #PAD AtomNames with 81 spaces 952 dict['AtomNames'] = dict['AtomNames'] + 81*' ' 953 dict['AtomSym'] = dict['AtomSym'] + 81*' ' 954 dict['AtomTree'] = dict['AtomTree'] + 81*' ' 955 dict['ResNames'] = dict['ResNames'] + 81*' ' 956 957 # create Iac list 958 #iac = [] 959 #tl = self.uniqTypeList 960 #for a in allAtoms: 961 # iac.append( tl.index(a.amber_symbol) + 1 ) 962 # delattr(a, 'amber_symbol') 963 #dict['Iac'] = iac 964 965 # to find out the number of bonds with hydrogen we simply count the 966 # number of hydrogen atoms 967 hlist = allAtoms.get(lambda x: x.element=='H') 968 if hlist is not None and len(hlist): 969 dict['Nbonh'] = numHs = len(hlist) 970 else: 971 numHs = 0 972 973 # number of bonds not involving an H atom 974 bonds = allAtoms.bonds[0] 975 dict['Mbona'] = len(bonds) - numHs 976 977 # since no bonds are constrined, Nbona==Mbona 978 dict['Nbona'] = dict['Mbona'] 979 980 print 'after call to processChain' 981 982 # new process bond info 983 dict['BondAt1'] = [] 984 dict['BondAt2'] = [] 985 dict['BondNum'] = [] 986 dict['BondHAt1'] = [] 987 dict['BondHAt2'] = [] 988 dict['BondHNum'] = [] 989 dict['Rk'] = [] 990 dict['Req'] = [] 991 self.processBonds(bonds, parmDict) 992 print 'after call to processBonds' 993 994 # now process the angles 995 dict['AngleAt1'] = [] 996 dict['AngleAt2'] = [] 997 dict['AngleAt3'] = [] 998 dict['AngleNum'] = [] 999 dict['AngleHAt1'] = [] 1000 dict['AngleHAt2'] = [] 1001 dict['AngleHAt3'] = [] 1002 dict['AngleHNum'] = [] 1003 dict['Tk'] = [] 1004 dict['Teq'] = [] 1005 self.processAngles(allAtoms, parmDict) 1006 print 'after call to processAngles' 1007 1008 # now handle the torsions 1009 dict['Nhparm'] = 0 1010 dict['Nparm'] = 0 1011 1012 dict['DihAt1'] = [] 1013 dict['DihAt2'] = [] 1014 dict['DihAt3'] = [] 1015 dict['DihAt4'] = [] 1016 dict['DihNum'] = [] 1017 dict['DihHAt1'] = [] 1018 dict['DihHAt2'] = [] 1019 dict['DihHAt3'] = [] 1020 dict['DihHAt4'] = [] 1021 dict['DihHNum'] = [] 1022 dict['Pn'] = [] 1023 dict['Pk'] = [] 1024 dict['Phase'] = [] 1025 1026 dict['Nphih'] = dict['Mphia'] = dict['Nphia'] = dict['Nptra'] = 0 1027 dict['N14pairs'] = [] 1028 dict['N14pairlist'] = [] 1029 1030 dict['Nnb'] =0 1031 dict['Iblo'] = [] 1032 dict['ExclAt'] = [] 1033 # FIXME 1034 self.AA5rings ={ 1035 'PRO':['N', 'CA', 'CB', 'CG', 'CD'], 1036 'TRP':['CG', 'CD1', 'CD2', 'NE1', 'CE2'], 1037 'HID':['CG', 'ND1', 'CE1', 'NE2', 'CD2'], 1038 'HIE':['CG', 'ND1', 'CE1', 'NE2', 'CD2'], 1039 'HIP':['CG', 'ND1', 'CE1', 'NE2', 'CD2'] 1040 } 1041 self.processTorsions(allAtoms, parmDict) 1042 print 'after call to processTorsions' 1043 1044 # some unused values 1045 dict['Nspm'] = 1 1046 dict['Box'] = [0., 0., 0.] 1047 dict['Boundary'] = [natom] 1048 dict['TreeJoin'] = range(natom) 1049 dict['Nphb'] = 0 1050 dict['HB12'] = [] 1051 dict['HB10'] = [] 1052 llist = ['Ifpert', 'Nbper','Ngper','Ndper','Mbper', 'Mgper', 1053 'Mdper','IfBox', 'IfCap', 'Cutcap', 'Xcap', 'Ycap', 1054 'Zcap', 'Natcap','Ipatm', 'Nspsol','Iptres'] 1055 for item in llist: 1056 dict[item] = 0 1057 1058 dict['Cno'] = self.getICO( dict['Ntypes'] ) 1059 dict['Solty'] = self.getSOLTY( dict['Natyp'] ) 1060 1061 dict['Cn1'], dict['Cn2'] = self.getCNList(parmDict) 1062 1063 return

1064 1065

1066 - def getICO(self, ntypes):

1067 ct = 1 1068 icoArray = Numeric.zeros((ntypes, ntypes), 'i') 1069 for i in range(1, ntypes+1): 1070 for j in range(1, i+1): 1071 icoArray[i-1,j-1]=ct 1072 icoArray[j-1,i-1]=ct 1073 ct = ct+1 1074 return icoArray.flat.tolist()

1075 1076

1077 - def getSOLTY(self, ntypes):

1078 soltyList = [] 1079 for i in range(ntypes): 1080 soltyList.append(0.) 1081 return soltyList

1082 1083

1084 - def getCN(self, type1, type2, pow, parmDict, factor=1):

1085 #pow is 12 or 6 1086 #factor is 1 except when pow is 6 1087 d = parmDict.potParam 1088 if type1=='N3': type1='N ' 1089 if type2=='N3': type2='N ' 1090 r1, eps1 = d[type1][:2] 1091 r2, eps2 = d[type2][:2] 1092 eps = sqrt(eps1*eps2) 1093 rij = r1 + r2 1094 newval = factor*eps*rij**pow 1095 return newval

1096 1097

1098 - def getCNList(self, parmDict):

1099 ntypes = len(self.uniqTypeList) 1100 ct = 1 1101 ## size = self.prmDict['Nttyp'] 1102 ## cn1List = [0]*size 1103 ## cn2List = [0]*size 1104 ## iac = self.prmDict['Iac'] 1105 ## cno = self.prmDict['Cno'] 1106 ## for i in range(ntypes): 1107 ## indi = i*ntypes 1108 ## ival = self.uniqTypeList[i] 1109 ## for j in range(i, ntypes): 1110 ## jval = self.uniqTypeList[j] 1111 ## ind = cno[indi+j]-1 1112 ## cn1List[ind] = self.getCN(jval, ival, 12, parmDict) 1113 ## cn2List[ind] = self.getCN(jval, ival, 6, parmDict, 2) 1114 nttyp = self.prmDict['Nttyp'] 1115 cn1List = [] 1116 cn2List = [] 1117 for j in range(ntypes): 1118 jval = self.uniqTypeList[j] 1119 for i in range(j+1): 1120 ival = self.uniqTypeList[i] 1121 cn1List.append(self.getCN(ival, jval, 12, parmDict)) 1122 cn2List.append(self.getCN(ival, jval, 6, parmDict, 2)) 1123 1124 return cn1List, cn2List

1125 1126

1127 - def readSummary(self, allLines, dict):

1128 #set summary numbers 1129 ll = split(allLines[1]) 1130 assert len(ll)==12 1131 #FIX THESE NAMES! 1132 natom = dict['Natom'] = int(ll[0]) 1133 ntypes = dict['Ntypes'] = int(ll[1]) 1134 nbonh = dict['Nbonh'] = int(ll[2]) 1135 dict['Mbona'] = int(ll[3]) 1136 ntheth = dict['Ntheth'] = int(ll[4]) 1137 dict['Mtheta'] = int(ll[5]) 1138 nphih = dict['Nphih'] = int(ll[6]) 1139 dict['Mphia'] = int(ll[7]) 1140 dict['Nhparm'] = int(ll[8]) 1141 dict['Nparm'] = int(ll[9]) 1142 #called 'next' in some documentation 1143 #NEXT-> Nnb 1144 next = dict['Nnb'] = int(ll[10]) 1145 dict['Nres'] = int(ll[11]) 1146 ll = split(allLines[2]) 1147 assert len(ll)==12 1148 nbona = dict['Nbona'] = int(ll[0]) 1149 ntheta = dict['Ntheta'] = int(ll[1]) 1150 nphia = dict['Nphia'] = int(ll[2]) 1151 numbnd = dict['Numbnd'] = int(ll[3]) 1152 numang = dict['Numang'] = int(ll[4]) 1153 numptra = dict['Nptra'] = int(ll[5]) 1154 natyp = dict['Natyp'] = int(ll[6]) 1155 dict['Nphb'] = int(ll[7]) 1156 dict['Ifpert'] = int(ll[8]) 1157 dict['Nbper'] = int(ll[9]) 1158 dict['Ngper'] = int(ll[10]) 1159 dict['Ndper'] = int(ll[11]) 1160 ll = split(allLines[3]) 1161 assert len(ll)==6 1162 dict['Mbper'] = int(ll[0]) 1163 dict['Mgper'] = int(ll[1]) 1164 dict['Mdper'] = int(ll[2]) 1165 dict['IfBox'] = int(ll[3]) 1166 dict['Nmxrs'] = int(ll[4]) 1167 dict['IfCap'] = int(ll[5]) 1168 return dict

1169 1170

1171 - def readIGRAPH(self, allLines, numIGRAPH, ind=3):

1172 #the names are not necessarily whitespace delimited 1173 igraph = [] 1174 for i in range(numIGRAPH): 1175 ind = ind + 1 1176 l = allLines[ind] 1177 for k in range(20): 1178 it = l[k*4:k*4+4] 1179 igraph.append(strip(it)) 1180 #igraph.extend(split(l)) 1181 return igraph, ind

1182 1183

1184 - def readCHRG(self, allLines, ind, numCHRG, natom):

1185 chrg = [] 1186 ct = 0 1187 for i in range(numCHRG): 1188 ind = ind + 1 1189 l = allLines[ind] 1190 newl = [] 1191 # build 5 charges per line if enough are left 1192 #otherwise, build the last line's worth 1193 if natom - ct >=5: 1194 rct = 5 1195 else: 1196 rct = natom - ct 1197 for q in range(rct): 1198 lindex = q*16 1199 item = l[lindex:lindex+16] 1200 newl.append(float(item)) 1201 ct = ct + 1 1202 chrg.extend(newl) 1203 return chrg, ind

1204 1205

1206 - def readNUMEX(self, allLines, ind, numIAC):

1207 numex = [] 1208 NumexSUM = 0 1209 for i in range(numIAC): 1210 ind = ind + 1 1211 ll = split(allLines[ind]) 1212 newl = [] 1213 for item in ll: 1214 newent = int(item) 1215 newl.append(newent) 1216 NumexSUM = NumexSUM + newent 1217 numex.extend(newl) 1218 return numex, ind, NumexSUM

1219 1220

1221 - def readLABRES(self, allLines, ind):

1222 done = 0 1223 labres = [] 1224 while not done: 1225 ind = ind + 1 1226 ll = split(allLines[ind]) 1227 try: 1228 int(ll[0]) 1229 done = 1 1230 break 1231 except ValueError: 1232 labres.extend(ll) 1233 #correct for 1 extra line read here 1234 ind = ind - 1 1235 return labres, ind

1236 1237

1238 - def readFList(self, allLines, ind, numITEMS):

1239 v = [] 1240 for i in range(numITEMS): 1241 ind = ind + 1 1242 ll = split(allLines[ind]) 1243 newl = [] 1244 for item in ll: 1245 newl.append(float(item)) 1246 v.extend(newl) 1247 return v, ind

1248 1249

1250 - def readIList(self, allLines, ind, numITEMS):

1251 v = [] 1252 for i in range(numITEMS): 1253 ind = ind + 1 1254 ll = split(allLines[ind]) 1255 newl = [] 1256 for item in ll: 1257 newl.append(int(item)) 1258 v.extend(newl) 1259 return v, ind

1260 1261

1262 - def readILList(self, allLines, ind, numITEMS, n):

1263 bhlist = [] 1264 for i in range(n): 1265 bhlist.append([]) 1266 ct = 0 1267 for i in range(numITEMS): 1268 ind = ind + 1 1269 ll = split(allLines[ind]) 1270 for j in range(len(ll)): 1271 item = ll[j] 1272 newl = bhlist[ct%n] 1273 newl.append(int(item)) 1274 ct = ct + 1 1275 return bhlist, ind

1276 1277

1278 - def py_read(self, filename, **kw):

1279 #??dict['Iptres'] #dict['Nspsol'] #dict['Ipatm'] #dict['Natcap'] 1280 f = open(filename, 'r') 1281 allLines = f.readlines() 1282 f.close() 1283 dict = {} 1284 #set title 1285 dict['ititl'] = allLines[0] 1286 1287 #get summary numbers 1288 dict = self.readSummary(allLines, dict) 1289 1290 #set up convenience fields: 1291 natom = dict['Natom'] 1292 ntypes = dict['Ntypes'] 1293 dict['Nat3'] = natom * 3 1294 dict['Ntype2d'] = ntypes ** 2 1295 nttyp = dict['Nttyp'] = ntypes * (ntypes+1)/2 1296 1297 # read IGRAPH->AtomNames 1298 numIGRAPH = int(ceil((natom*1.)/20.)) 1299 anames, ind = self.readIGRAPH(allLines, numIGRAPH) 1300 dict['AtomNames'] = join(anames) 1301 1302 # read CHRG->Charges 1303 numCHRG = int(ceil((natom*1.)/5.)) 1304 dict['Charges'], ind = self.readCHRG(allLines, ind, numCHRG, natom) 1305 1306 # read AMASS **same number of lines as charges->Masses 1307 dict['Masses'], ind = self.readFList(allLines, ind, numCHRG) 1308 1309 # read IAC **NOT same number of lines as IGRAPH 12!! 1310 numIAC = int(ceil((natom*1.)/12.)) 1311 dict['Iac'], ind = self.readIList(allLines, ind, numIAC) 1312 1313 # read NUMEX **same number of lines as IAC 1314 dict['Iblo'], ind, NumexSUM = self.readNUMEX(allLines, ind, numIAC) 1315 1316 # read ICO *Ntype2d/12 1317 numICO = int(ceil((ntypes**2*1.0)/12.)) 1318 dict['Cno'], ind = self.readIList(allLines, ind, numICO) 1319 1320 ##NB this should be half of a matrix 1321 # read LABRES....no way to know how many 1322 dict['ResNames'], ind = self.readLABRES(allLines, ind) 1323 labres = dict['ResNames'] 1324 1325 # read IPRES....depends on len of LABRES 1326 numIPRES = int(ceil((len(labres)*1.)/20.)) 1327 dict['Ipres'], ind = self.readIList(allLines, ind, numIPRES) 1328 1329 # read RK + REQ-> depend on numbnd 1330 numbnd = dict['Numbnd'] 1331 numRK = int(ceil((numbnd*1.)/5.)) 1332 dict['Rk'], ind = self.readFList(allLines, ind, numRK) 1333 dict['Req'], ind = self.readFList(allLines, ind, numRK) 1334 1335 # read TK + TEQ-> depend on numang 1336 numang = dict['Numang'] 1337 numTK = int(ceil((numang*1.)/5.)) 1338 dict['Tk'], ind = self.readFList(allLines, ind, numTK) 1339 dict['Teq'], ind = self.readFList(allLines, ind, numTK) 1340 1341 # read PK, PN + PHASE-> depend on numptra 1342 nptra = dict['Nptra'] 1343 numPK = int(ceil((nptra*1.)/5.)) 1344 dict['Pk'], ind = self.readFList(allLines, ind, numPK) 1345 dict['Pn'], ind = self.readFList(allLines, ind, numPK) 1346 dict['Phase'], ind = self.readFList(allLines, ind, numPK) 1347 1348 # read SOLTY 1349 natyp = dict['Natyp'] 1350 numSOLTY = int(ceil((natyp*1.)/5.)) 1351 dict['Solty'], ind = self.readFList(allLines, ind, numSOLTY) 1352 1353 # read CN1 and CN2 1354 numCN = int(ceil((nttyp*1.)/5.)) 1355 dict['Cn1'], ind = self.readFList(allLines, ind, numCN) 1356 dict['Cn2'], ind = self.readFList(allLines, ind, numCN) 1357 1358 # read IBH, JBH, ICBH 12 1359 nbonh = dict['Nbonh'] 1360 numIBH = int(ceil((nbonh*3.0)/12.)) 1361 [dict['BondHAt1'], dict['BondHAt2'], dict['BondHNum']], ind = \ 1362 self.readILList(allLines, ind, numIBH, 3) 1363 # read IB, JB, ICB 12 1364 nbona = dict['Nbona'] 1365 numIB = int(ceil((nbona*3.0)/12.)) 1366 [dict['BondAt1'], dict['BondAt2'], dict['BondNum']], ind = \ 1367 self.readILList(allLines, ind, numIB, 3) 1368 1369 # read ITH, JTH, KTH, ICTH 12 1370 ntheth = dict['Ntheth'] 1371 numITH = int(ceil((ntheth*4.0)/12.)) 1372 [dict['AngleHAt1'], dict['AngleHAt2'], dict['AngleHAt3'],\ 1373 dict['AngleHNum']], ind = self.readILList(allLines, ind, numITH, 4) 1374 # read IT, JT, KT, ICT 12 1375 ntheta = dict['Ntheta'] 1376 numIT = int(ceil((ntheta*4.0)/12.)) 1377 [dict['AngleAt1'], dict['AngleAt2'], dict['AngleAt3'],\ 1378 dict['AngleNum']], ind = self.readILList(allLines, ind, numIT, 4) 1379 1380 # read IPH, JPH, KPH, LPH, ICPH 12 1381 nphih = dict['Nphih'] 1382 numIPH = int(ceil((nphih*5.0)/12.)) 1383 [dict['DihHAt1'], dict['DihHAt2'], dict['DihHAt3'], dict['DihHAt4'],\ 1384 dict['DihHNum']], ind = self.readILList(allLines, ind, numIPH, 5) 1385 # read IP, JP, KP, LP, ICP 12 1386 nphia = dict['Nphia'] 1387 numIP = int(ceil((nphia*5.0)/12.)) 1388 [dict['DihAt1'], dict['DihAt2'], dict['DihAt3'], dict['DihAt4'],\ 1389 dict['DihNum']], ind = self.readILList(allLines, ind, numIP, 5) 1390 1391 # read NATEX 12 1392 #FIX THIS: has to be the sum of previous entries 1393 numATEX = int(ceil((NumexSUM*1.0)/12.)) 1394 dict['ExclAt'], ind = self.readIList(allLines, ind, numATEX) 1395 1396 # read CN1 and CN2 1397 # skip ASOL 1398 # skip BSOL 1399 # skip HBCUT 1400 ind = ind + 3 1401 # read ISYMBL 20 1402 asym, ind = self.readIGRAPH(allLines, numIGRAPH, ind) 1403 dict['AtomSym'] = join(asym) 1404 1405 # read ITREE 20 1406 atree, ind = self.readIGRAPH(allLines, numIGRAPH, ind) 1407 dict['AtomTree'] = join(atree) 1408 return dict

1409 1410

1411 - def makeList(self, llist, num):

1412 newL = [] 1413 for i in range(len(llist[0])): 1414 ni = [] 1415 for j in range(num): 1416 ni.append(llist[j][i]) 1417 newL.append(ni) 1418 return newL

1419 1420 1421 #functions to write self

1422 - def write(self, filename, **kw):

1423 fptr = open(filename, 'w') 1424 dict = self.prmDict 1425 self.writeItitl(fptr, dict['ititl']) 1426 self.writeSummary(fptr) 1427 #WHAT ABOUT SOLTY??? 1428 self.writeString(fptr,dict['AtomNames'][:-81]) 1429 for k in ['Charges', 'Masses', 'Iac','Iblo','Cno']: 1430 item = dict[k] 1431 f = self.formatD[k] 1432 if f[2]: 1433 self.writeTupleList(fptr, item, f[0], f[1], f[2]) 1434 else: 1435 self.writeList(fptr, item, f[0], f[1]) 1436 self.writeString(fptr,dict['ResNames'][:-81]) 1437 self.writeList(fptr, dict['Ipres'][:-1], '%6d', 12 ) 1438 for k in ['Rk', 'Req', 'Tk', 'Teq', 1439 'Pk', 'Pn', 'Phase', 'Solty', 'Cn1','Cn2']: 1440 item = dict[k] 1441 f = self.formatD[k] 1442 if f[2]: 1443 self.writeTupleList(fptr, item, f[0], f[1], f[2]) 1444 else: 1445 self.writeList(fptr, item, f[0], f[1]) 1446 #next write bnds angs and dihe 1447 allHBnds = zip(dict['BondHAt1'], dict['BondHAt2'], 1448 dict['BondHNum']) 1449 self.writeTupleList(fptr, allHBnds, "%6d", 12, 3) 1450 allBnds = zip(dict['BondAt1'], dict['BondAt2'], 1451 dict['BondNum']) 1452 self.writeTupleList(fptr, allBnds, "%6d", 12, 3) 1453 1454 allHAngs = zip(dict['AngleHAt1'], dict['AngleHAt2'], 1455 dict['AngleHAt3'], dict['AngleHNum']) 1456 self.writeTupleList(fptr, allHAngs, "%6d", 12,4) 1457 allAngs = zip(dict['AngleAt1'], dict['AngleAt2'], 1458 dict['AngleAt3'], dict['AngleNum']) 1459 self.writeTupleList(fptr, allAngs, "%6d", 12, 4) 1460 1461 allHDiHe = zip(dict['DihHAt1'], dict['DihHAt2'], 1462 dict['DihHAt3'], dict['DihHAt4'], dict['DihHNum']) 1463 self.writeTupleList(fptr, allHDiHe, "%6d", 12,5) 1464 allDiHe = zip(dict['DihAt1'], dict['DihAt2'], 1465 dict['DihAt3'], dict['DihAt4'], dict['DihNum']) 1466 self.writeTupleList(fptr, allDiHe, "%6d", 12, 5) 1467 self.writeList(fptr, dict['ExclAt'], '%6d', 12) 1468 fptr.write('\n') 1469 fptr.write('\n') 1470 fptr.write('\n') 1471 for k in ['AtomSym', 'AtomTree']: 1472 item = dict[k][:-81] 1473 self.writeString(fptr, item) 1474 zList = [] 1475 for i in range(dict['Natom']): 1476 zList.append(0) 1477 self.writeList(fptr, zList, "%6d", 12) 1478 self.writeList(fptr, zList, "%6d", 12) 1479 fptr.close()

1480 1481

1482 - def writeString(self, fptr, item):

1483 n = int(ceil(len(item)/80.)) 1484 for p in range(n): 1485 if p!=n-1: 1486 fptr.write(item[p*80:(p+1)*80]+'\n') 1487 else: 1488 #write to the end, whereever it is 1489 fptr.write(item[p*80:]+'\n')

1490 1491

1492 - def writeList(self, fptr, outList, formatStr="%4.4s", lineCt=12):

1493 ct = 0 1494 s = "" 1495 nlformatStr = formatStr+'\n' 1496 lenList = len(outList) 1497 for i in range(lenList): 1498 #do something with outList[i] 1499 s = s + formatStr%outList[i] 1500 #ct is how many item are in s 1501 ct = ct + 1 1502 #if line is full, write it and reset s and ct 1503 if ct%lineCt==0: 1504 s = s + '\n' 1505 fptr.write(s) 1506 s = "" 1507 ct = 0 1508 #if last entry write it and exit 1509 elif i == lenList-1: 1510 s = s + '\n' 1511 fptr.write(s) 1512 break

1513 1514

1515 - def writeTupleList(self, fptr, outList, formatStr="%4.4s", lineCt=12, ll=2):

1516 ct = 0 1517 s = "" 1518 nlformatStr = formatStr+'\n' 1519 for i in range(len(outList)): 1520 if i==len(outList)-1: 1521 for k in range(ll): 1522 s = s + formatStr%outList[i][k] 1523 ct = ct + 1 1524 if ct%lineCt==0: 1525 s = s + '\n' 1526 fptr.write(s) 1527 s = "" 1528 ct = 0 1529 #after adding last entry, if anything left, print it 1530 if ct!=0: 1531 s = s + '\n' 1532 fptr.write(s) 1533 else: 1534 for k in range(ll): 1535 s = s + formatStr%outList[i][k] 1536 ct = ct + 1 1537 if ct%lineCt==0: 1538 s = s + '\n' 1539 fptr.write(s) 1540 s = "" 1541 ct = 0

1542

1543 - def writeItitl(self, fptr, ititl):

1544 fptr.write(ititl)

1545 1546

1547 - def writeSummary(self, fptr):

1548 #SUMMARY 1549 #fptr.write('SUMMARY\n') 1550 ##FIX THESE NAMES!!! 1551 kL1 = ['Natom','Ntypes','Nbonh','Mbona',\ 1552 'Ntheth','Mtheta','Nphih','Mphia','Nhparm',\ 1553 'Nparm','Nnb','Nres'] 1554 kL2 = ['Nbona','Ntheta','Nphia','Numbnd',\ 1555 'Numang','Nptra','Natyp','Nphb','Ifpert',\ 1556 'Nbper','Ngper','Ndper'] 1557 kL3 = ['Mbper','Mgper','Mdper','IfBox','Nmxrs',\ 1558 'IfCap'] 1559 1560 for l in [kL1, kL2, kL3]: 1561 newL = [] 1562 for item in l: 1563 newL.append(self.prmDict[item]) 1564 #print 'newL=', newL 1565 self.writeList(fptr, newL, "%6d", 12)

Source Code for Module MolKit.amberPrmTop