######################################################################### # # # This is a module that contains functions generally useful for the # # MMPBSA.py script. A full list of functions/subroutines is shown below.# # It must be included to insure proper functioning of MMPBSA.py # # # # Last updated: 05/21/2010 # # # ######################################################################### ########################## GPL LICENSE INFO ############################ # Copyright (C) 2009 Dwight McGee, Billy Miller III, and Jason Swails # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. ######################################################################## # List of functions and a brief description of their purpose # # which: finds an executable in the PATH much like the unix command # fileexists: determines whether a passed file is found in the current # directory # fileexists_noprint: determines whether a passed file is found in the # current directory, but does not print an error # message # linenum: returns the number of lines in a given file # resnum: returns the number of residues defined in a topology file # numatm: returns the number of atoms defined in a topology file # remove: deletes files created by MMPBSA.py. Allows 3 levels of # removal (total, some, none) # getResStrings: for decomp - generate strings needed for mdin files # gbmdin: writes the mdin file for the gb MMPBSA simulations # pbmdin: writes the mdin file for the pb MMPBSA simulations # pbmdin_old: writes the mdin file for the pb MMPBSA simulations (old pbsa) # decomp_gbmdin: writes the mdin file for the decomp gb MMPBSA simulations # decomp_pbmdin: writes the mdin file for the decomp pb MMPBSA simulations # getenergies: calculates each component of an energy in a gb mdout # file as well as the total # getpbenergies: same as getenergies above, but for a pb mdout file # round: rounds a given number to a given decimal point # average: returns the average of a list # stdev: returns the standard deviation of a list # stdev_cume: runs a cumulative standard deviation # printusage: prints the command line usage of MMPBSA.py to the screen # digit: returns the given digit of a number. # CheckIncomp: Checks for incompatibilities in the input variables. # CheckRadiis: Checks that the radii set is compatible with igb choice # and either warns or quits # getresinfo: Gets the residue name of a specific residue number in a # given prmtop file # getbondiset: Gets the radii set used to set up a given prmtop # getnumatms: Returns the number of atoms corresponding to a specific # residue name # getressymbol: Returns the 1-letter symbol for a specified amino acid # Format: Formats a floating point decimal to be used in output # (specify decimals and leading spaces) # ScaleDistance: Gets 6 coords (2 atoms) and adjusts bond length to the # specified distance, adjusts only 2nd set of 3 coords. # getallresinfo: Takes a given flag and returns every value in given # prmtop corresponding to that flag # CheckMutParm: Checks to make sure the complex mutant top file and # regular complex top file differ in one specific residue # CheckRecLigParms: Checks to make sure receptor and ligand mutant top # files are also consistent # getligmask: Finds the ligand mask based on given prm files # getrecmask: Finds receptor mask from ligand mask and complex prmtop # maskToDec: Converts an amber mask to a decomposition string # GetMasks: returns the ligand and receptor masks # PtrajNumFrames: Finds out how many frames have been written based on # a ptraj output. # nmodecalc: sets up and executes an nmode calculation on a series of 3 # simulations # gbcalc: sets up and executes gb calculations # pbcalc: sets up and executes pb calculations # decomp_out: wrapper for decomp output. calls internal functions # gboutput: prints final results from gb calculations # pboutput: prints final results from pb calculations # decoutput12: parses and prints decomp output idecomp=1or2 # decoutput34: parses and prints decomp output idecomp=3or4 # separate: prints field separator in final output file # GetPtrajEntropy: parses the ptraj entropy calculation file for the # parameters # GetNmodeEntropy: parses the nmode entropy calculation files for the # parameters # PrintPtrajEnt: prints the ptraj entropy results to the final output # PrintNmodeEnt: prints the nmode entropy results to the final output # PrintFinalResults: prints the final results of calculations # ######################################################################## #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def which(program): import os, sys def is_exe(fpath): # determines if a file is executable or not return os.path.exists(fpath) and os.access(fpath, os.X_OK) fpath, fname = os.path.split(program) if fpath: if is_exe(program): return program # if found, return program else: # if above doesn't work, try looking for program in every folder in PATH for path in os.environ["PATH"].split(os.pathsep): exe_file = os.path.join(path, program) if is_exe(exe_file): return exe_file # if it's executable, return the file try: # now look through AMBERHOME/exe/_program_, but catch exception if AMBERHOME not defined if is_exe(os.environ["AMBERHOME"] + "/exe/" + program): return os.environ["AMBERHOME"] + "/exe/" + program except KeyError: print >> sys.stderr, "Error: AMBERHOME is not defined!" return 'none' # if program can still not be found... return "none" #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def fileexists(file): import sys try: f = open(file,'r') # try to open the file for reading except IOError: # catch exception if it's not there print >> sys.stderr, 'Error: Specified file ({0}) does not exist!'.format(file) return -1 # return error code f.close() # if it exists, close it up and return 0 return 0 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def fileexists_noprint(file): try: f = open(file,'r') # try to open the file for reading except IOError: # catch exception if it's not there and return error code return -1 f.close() # if it's there, close up the file and return 0 return 0 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def linenum(file): try: f = open(file,'r') # open the file for reading, catching if it's not there except IOError: return -1 # count number of lines by literally incrementing a counter for each line in the file i = 0 for line in f: i += 1 f.close() return i #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def resnum(topfile): import sys try: top = open(topfile, 'r') except IOError: print >> sys.stderr, 'Error: topology file {0} does not exist!'.format(topfile) return -1 lines = top.readlines() top.close() # read the number of residues for x in range(len(lines)): if lines[x].startswith('%FLAG POINTERS'): # once we hit this, we know where to find it x = x + 1 # go to the next line while not lines[x].startswith('%FORMAT'): # skip over comments x = x + 1 words = lines[x+2].split() # skip over FORMAT line, and it's the 2nd number on the next line return int(words[1]) # this is only reached if FLAG POINTERS was never found print >> sys.stderr, 'Error: Could not find %FLAG POINTERS in topology file {0}!'.format(topfile) return -1 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def numatm(topfile): # returns number of atoms in a prmtop import sys try: top = open(topfile, 'r') except IOError: return -1 lines = top.readlines() # read lines into memory top.close() for x in range(len(lines)): if lines[x].startswith('%FLAG POINTERS'): # it is in this block x = x + 1 # go to next line while not lines[x].startswith('%FORMAT'): # skip over comments x = x + 1 words = lines[x+1].split() # line after %FORMAT return int(words[0]) # it's the first number # this is only reached if FLAG POINTERS was not found print >> sys.stderr, 'Error: Could not find %FLAG POINTERS in topology file {0}!'.format(topfile) return -1 # return error code #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def remove(flag,mpi_size): import os, sys if flag == -3: # internal -- keep all mdin files. cmd = 'tar -cf tmp.tar _MMPBSA_*.mdin; rm -f _MMPBSA_* io.mc restrt; tar -xf tmp.tar; rm tmp.tar' os.system(cmd) if flag == -2: # internal -- if program exits in error, keep files and say so print >> sys.stderr, 'NOTE: All files have been retained for ' + \ 'debugging purposes. Type MMPBSA.py --clean to erase these files.' sys.exit() elif flag == 0 or flag == -1: os.system('rm -f _MMPBSA_* io.mc restrt') # remove all temp files. -1 is depreciated and now equivalent to 0 elif flag == 1: # keep mdcrds, mdouts, and other relevant output files. os.system('rm -f _MMPBSA_*.inpcrd _MMPBSA_*.in _MMPBSA_*mdin io.mc restrt') for x in range(mpi_size): # remove all temp files created by individual threads in MMPBSA.py.MPI os.system('rm -f _MMPBSA_*.{0}'.format(x)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getResStrings(comprm, recprm, ligprm, selection, ligstart): def improperuse(code): import sys if code == 1: print >> sys.stderr, "Error: Improper definition of print_res!" elif code == 2: print >> sys.stderr, "Error: You specified a residue multiple times in print_res!" return -1 lignum = resnum(ligprm) # number of residues in ligand prmtop recnum = resnum(recprm) # number of residues in receptor prmtop comnum = resnum(comprm) # number of residues in complex prmtop comsel = '' # string containing the complex selection recsel = '' # string containing the receptor selection ligsel = '' # string containing the ligand selection # parse the selection if ligstart == -1: return ["EDIT","EDIT","EDIT"] # if we didn't guess the mask, we can't go through this if selection == "none": return ["1 " + str(comnum), "1 " + str(recnum), "1 " + str(lignum)] # default, all residues else: items = selection.split(';') # otherwise, split along delimiter try: for x in range(len(items)): # loop through all delimited fields items[x] = items[x].strip() # get rid of whitespace if len(items[x]) == 0: # skip over any blank fields continue individual = items[x].split('-') # break each field up by - comsel += individual[0].strip() + ' ' # the first number is always in complex if len(individual) > 1: comsel += individual[1].strip() + ' ' # put 2nd number if it's there else: comsel += individual[0].strip() + ' ' # otherwise repeat the 1st number # Now it's time to add the residues to the receptor/ligand selections if len(individual) == 1: # if there is only 1 res in this field individual[0] = int(individual[0].strip()) # turn it into an integer if individual[0] < 0 or individual[0] > comnum: # catch blatantly illegal choices return improperuse(1) # return error code if individual[0] >= ligstart and individual[0] < ligstart + lignum: # if it's in ligand individual[0] = individual[0] - ligstart + 1 # adjust number to its value in ligand prmtop ligsel += str(individual[0]) + ' ' + str(individual[0]) + ' ' # add it to ligand selection elif individual[0] < ligstart: # if it comes before ligand in receptor ... recsel += str(individual[0]) + ' ' + str(individual[0]) + ' ' # add it to receptor_selection as-is else: # otherwise... recsel += str(individual[0] - lignum) + ' ' + str(individual[0] - lignum) + ' ' # adjust for absence of ligand else: # if it's a field like " num-num " individual[0] = int(individual[0].strip()) # make both numbers integers individual[1] = int(individual[1].strip()) if individual[0] > individual[1] or individual[0] <= 0 or individual[1] <= 0 or \ individual[0] > comnum or individual[1] > comnum: # catch blatantly illegal choices return improperuse(1) # return error code if individual[0] < ligstart and individual[1] < ligstart: # if it starts and ends before start of ligand recsel += str(individual[0]) + ' ' + str(individual[1]) + ' ' # add it to receptor selection elif individual[0] < ligstart and \ (individual[1] >= ligstart and individual[1] < ligstart + lignum): # starts before ligand ends in ligand recsel += str(individual[0]) + ' ' + str(ligstart - 1) + ' ' ligsel += '1 ' + str(individual[1] - ligstart + 1) elif (individual[0] >= ligstart and individual[0] < ligstart + lignum) and \ (individual[1] >= ligstart and individual[1] < ligstart + lignum): # starts and ends in ligand ligsel += str(individual[0] - ligstart + 1) + ' ' + \ str(individual[1] - ligstart + 1) + ' ' elif (individual[0] >= ligstart and individual[0] < ligstart + lignum) and \ (individual[1] >= ligstart + lignum): # starts in ligand ends after ligand (in receptor) recsel += str(ligstart) + ' ' + str(individual[1] - lignum) + ' ' ligsel += str(individual[0] - ligstart + 1) + ' ' + str(lignum) + ' ' else: # starts and ends after ligand recsel += str(individual[0] - lignum) + ' ' + str(individual[1] - lignum) + ' ' except ValueError: # catch any use of non-integers return improperuse(1) except IndexError: # catch any other weird errors coming from poorly-defined fields return improperuse(1) return [comsel, recsel, ligsel] # return the selections #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def gbmdin(INPUT): mdin = open('_MMPBSA_gb.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, idecomp=0, cut=999.0, nsnb=99999, imin=5, maxcyc=0, ncyc=0, igb={0}, saltcon={1}, intdiel={2}, extdiel={3}, gbsa={4}, surften={5}, / """.format(INPUT[9], INPUT[11], INPUT[13], INPUT[12], INPUT[10], INPUT[14])) mdin.close() #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def pbmdin(INPUT): # INPUT[37] is sander_apbs - 0: use PBSA, otherwise: use sander.APBS if INPUT[37] == '0': # we need to convert ionic strength from M to milliM for internal PBSA istrng = str(float(INPUT[21]) * 1000) mdin = open('_MMPBSA_pb.mdin','w') mdin.write("""File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, ntpr=1, imin=5, maxcyc=0, igb=10, ipb=1, inp={0}, / &pb dbfopt=1, epsin={1}, epsout={2}, istrng={3}, radiopt={4}, sprob={5}, space={6}, maxitn={7}, cavity_surften={8}, cavity_offset={9}, fillratio={10}, / """.format(INPUT[22],INPUT[16],INPUT[17],istrng,INPUT[28],INPUT[20],INPUT[18],INPUT[19],INPUT[23],\ INPUT[24],INPUT[25])) mdin.close() else: # we need to convert surften into kJ/(mol A^2) for sander.APBS. # Default amber energy is kcal. surften = float(INPUT[23]) * 4.184 mdin = open('_MMPBSA_pb.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl maxcyc=0, imin=5, cut=999.0, igb=6, ntb=0, ntpr=1, ntmin=2, / &apbs apbs_print=1, calc_type=0, cmeth=1, bcfl=2, srfm=1, chgm=1, pdie={0}, sdie={1}, srad={2}, nion=2, ionq=1.0, -1.0, ionc={3}, {3}, ionrr=2.0, 2.0, radiopt={4}, calcforce=0, calcnpenergy=1, grid={5},{5},{5}, gamma={6} / """.format(INPUT[16],INPUT[17],INPUT[20],INPUT[21],INPUT[28],INPUT[18],surften)) mdin.close() #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def pbmdin_old(INPUT): # we need to convert ionic strength from M to milliM for internal PBSA istrng = str(float(INPUT[21]) * 1000) mdin = open('_MMPBSA_pb.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=0, cutcap=-1, xcap=0, ycap=0, zcap=0, idecomp=0, ntpr=1, ivcap=0, ntmin=2, igb=10, / &pb dbfopt=1, epsin={0}, epsout={1}, istrng={2}, radiopt={3}, sprob={4}, space={5}, maxitn={6}, npopt={7}, cavity_surften={8}, cavity_offset={9}, fillratio={10} / """.format(INPUT[16],INPUT[17],istrng,INPUT[28],INPUT[20],INPUT[18],INPUT[19],INPUT[22],INPUT[23], INPUT[24],INPUT[25])) mdin.close() #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def decomp_gbmdin(INPUT, ligstart): # get the RES strings used for printing in the input files printstrings = getResStrings(INPUT[2], INPUT[3], INPUT[4], INPUT[55], ligstart) if printstrings == -1: return -1 mdin = open('_MMPBSA_gb_decomp_com.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=1, ncyc=0, igb={0}, saltcon={1}, intdiel={2}, extdiel={3}, gbsa={4}, surften={5}, idecomp={6}, offset=0.09, / Residues considered as REC RRES {7} END Residues considered as LIG LRES {8} END Residues to print RES {9} END END """.format(INPUT[9],INPUT[11],INPUT[13],INPUT[12],INPUT[10],INPUT[14],INPUT[54],maskToDec(INPUT[35]), maskToDec(INPUT[36]),printstrings[0])) mdin.close() mdin = open('_MMPBSA_gb_decomp_rec.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=1, ncyc=0, igb={0}, saltcon={1}, intdiel={2}, extdiel={3}, gbsa={4}, surften={5}, idecomp={6}, offset=0.09, / Residues considered as REC RRES 1 {7} END Residues to print RES {8} END END """.format(INPUT[9],INPUT[11],INPUT[13],INPUT[12],INPUT[10],INPUT[14],INPUT[54], resnum(INPUT[3]),printstrings[1])) mdin.close() mdin = open('_MMPBSA_gb_decomp_lig.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=1, ncyc=0, igb={0}, saltcon={1}, intdiel={2}, extdiel={3}, gbsa={4}, surften={5}, idecomp={6}, offset=0.09, / Residues considered as LIG LRES 1 {7} END Residues to print RES {8} END END """.format(INPUT[9],INPUT[11],INPUT[13],INPUT[12],INPUT[10],INPUT[14],INPUT[54], resnum(INPUT[4]),printstrings[2])) mdin.close() #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def decomp_pbmdin(INPUT, ligstart): # get the RES strings used for printing in the input files printstrings = getResStrings(INPUT[2], INPUT[3], INPUT[4], INPUT[55], ligstart) if printstrings == -1: return -1 # we need to convert ionic strength from M to milliM for internal PBSA istrng = str(float(INPUT[21]) * 1000) mdin = open('_MMPBSA_pb_decomp_com.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=0, cutcap=-1, xcap=0, ycap=0, zcap=0, idecomp=0, ntpr=1, ivcap=0, ntmin=2, igb=10, ipb=2, inp={0}, idecomp={1}, / &pb dbfopt=1, epsin={2}, epsout={3}, istrng={4}, radiopt={5}, sprob={6}, space={7}, maxitn={8}, npopt={9}, cavity_surften={10}, cavity_offset={11}, fillratio={12}, / Residues considered as REC RRES {13} END Residues considered as LIG LRES {14} END Residues to print RES {15} END END """.format(INPUT[22],INPUT[54],INPUT[16],INPUT[17],istrng,INPUT[28],INPUT[20],INPUT[18],INPUT[19],INPUT[22], INPUT[23],INPUT[24],INPUT[25],maskToDec(INPUT[35]),maskToDec(INPUT[36]),printstrings[0])) mdin.close() # npopt = inp = 0 --> npolar solvation is non-decomposable mdin = open('_MMPBSA_pb_decomp_rec.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=0, cutcap=-1, xcap=0, ycap=0, zcap=0, idecomp=0, ntpr=1, ivcap=0, ntmin=2, igb=10, ipb=2, inp={0}, idecomp={1}, / &pb dbfopt=1, epsin={2}, epsout={3}, istrng={4}, radiopt={5}, sprob={6}, space={7}, maxitn={8}, npopt={9}, cavity_surften={10}, cavity_offset={11}, fillratio={12}, / Residues considered as REC RRES 1 {13} END Residues to print RES {14} END END """.format(INPUT[22],INPUT[54],INPUT[16],INPUT[17],istrng,INPUT[28],INPUT[20],INPUT[18],INPUT[19],INPUT[22], INPUT[23],INPUT[24],INPUT[25],resnum(INPUT[3]),printstrings[1])) mdin.close() # npopt = inp = 0 --> npolar solvation is non-decomposable mdin = open('_MMPBSA_pb_decomp_lig.mdin','w') mdin.write( """File generated by MMPBSA.py &cntrl ntb=0, cut=999.0, nsnb=99999, imin=5, maxcyc=0, cutcap=-1, xcap=0, ycap=0, zcap=0, idecomp=0, ntpr=1, ivcap=0, ntmin=2, igb=10, ipb=2, inp={0}, idecomp={1}, / &pb dbfopt=1, epsin={2}, epsout={3}, istrng={4}, radiopt={5}, sprob={6}, space={7}, maxitn={8}, npopt={9}, cavity_surften={10}, cavity_offset={11}, fillratio={12}, / Residues considered as LIG LRES 1 {13} END Residues to print RES {14} END END """.format(INPUT[22],INPUT[54],INPUT[16],INPUT[17],istrng,INPUT[28],INPUT[20],INPUT[18],INPUT[19],INPUT[22], INPUT[23],INPUT[24],INPUT[25],resnum(INPUT[4]),printstrings[2])) mdin.close() # npopt = inp = 0 --> npolar solvation is non-decomposable #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getenergies(bond, angle, dihed, vdwaals, eel, egb, vdwof, eelof, esurf, total, outputfile, onetraj): import sys marker = 0 # a marker used to keep track of where we are in the mdout file # Parsing the mdout file. We are parsing potential terms that look like: # # NSTEP ENERGY RMS GMAX NAME NUMBER # 1 -5.0085E+03 1.7869E+01 1.0861E+02 C 3395 # # BOND = 747.3975 ANGLE = 2159.7999 DIHED = 2674.1632 # VDWAALS = -2025.6491 EEL = -17003.5331 EGB = -3086.6385 # 1-4 VDW = 905.5446 1-4 EEL = 10529.1160 RESTRAINT = 0.0000 # ESURF = 91.3192 for outputline in outputfile: # outputfile was passed as an open file object. Loop through the lines if 'error' in outputline.lower(): # look for errors in the sander mdout print >> sys.stderr, 'Error: Sander error ({0})'.format(outputline.strip()) return -1 if 'BOND' in outputline: # if the BOND potential term label is in this line: words = outputline.split() # split along whitespace marker = 1 # we are in the first line we are parsing values from try: bond.append(float(words[2])) # bond value is 3rd in array after BOND and = angle.append(float(words[5])) # angle is 6th dihed.append(float(words[8])) # dihedral is 9th except ValueError: # this will occur if the values are ********** (incompatible prmtop/mdcrds) print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 elif marker == 1: # this occurs the cycle after BOND was found. It is the next line words = outputline.split() # split along whitespace marker += 1 # move the marker along (should be 2 now) try: vdwaals.append(float(words[2])) # vdw is 3rd word after VDWAALS and = eel.append(float(words[5])) # eel is 6th egb.append(float(words[8])) # egb is 9th except ValueError: # find ********* in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # and the offending line and return error code return -1 elif marker == 2: # this occurs 2 lines after BOND was found. words = outputline.split() # split along whitespace marker = marker + 1 # move marker along try: vdwof.append(float(words[3])) # 1-4vdw is the 4th word (after 1-4, VDW, and =) eelof.append(float(words[7])) # 1-4 eel is the 8th word except ValueError: # find ********** in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 elif marker == 3: # this occurs 3 lines after BOND was found words = outputline.split() # split along whitespace marker = marker + 1 # move marker along try: esurf.append(float(words[2])) # esurf is the 3rd word after ESURF and = except ValueError: # find ******* in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 if len(bond) == 0: # if we went through whole mdout and didn't find BOND once... print >> sys.stderr, 'Error: No potential terms in sander output! Check output files.' return -1 # return error code after printing error message for x in range(len(bond)): # all arrays will be same length if script made it this far try: # add all potential terms sum = bond[x] + angle[x] + dihed[x] + vdwaals[x] + eel[x] + egb[x] + vdwof[x] + eelof[x] + esurf[x] except IndexError: # catch this very strange error if it occurs... it never should print >> sys.stderr, 'Error: Inconsistent number of potential terms in mdouts. This is a strange error' return -1 total.append(sum) # append sum to the total array return 1 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getpbenergies(bond, angle, dihed, vdwaals, eel, egb, vdwof, eelof, ecavity, edisper, total, outputfile, apbs, onetraj): import sys marker = 0 # the marker tracks where in the output block we are (0 for BOND, 1 for VDWAALS, etc.) # We are parsing potential terms from blocks that look like # NSTEP ENERGY RMS GMAX NAME NUMBER # 1 -5.2363E+03 1.7793E+01 1.1102E+02 C 3395 # # BOND = 747.3975 ANGLE = 2159.7999 DIHED = 2674.1632 # VDWAALS = -2025.6491 EEL = -17003.5331 EPB = -3290.2645 # 1-4 VDW = 905.5446 1-4 EEL = 10529.1160 RESTRAINT = 0.0000 # ECAVITY = 67.1578 EDISPER = 0.0000 # NOTE: for sander.APBS, ECAVITY is named ENPOLAR and EDISPER does not exist for outputline in outputfile: words = outputline.split() # split along whitespace if len(words) > 0: # only analyze non-blank lines if words[0] == 'BOND': # we have reached the bond marker = 1 # update marker to reflect this try: bond.append(float(words[2])) # bond term is 3rd word after BOND and = angle.append(float(words[5])) # angle is 6th word dihed.append(float(words[8])) # dihedral is 9th word except ValueError: # catch ******** in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 elif marker == 1: # this is hit on VDWAALS line marker = marker + 1 # move marker along try: vdwaals.append(float(words[2])) # vdwaals is 3rd word eel.append(float(words[5])) # eel is 6th word egb.append(float(words[8])) # egb is 9th word except ValueError: # catch ****** in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 elif marker == 2: # this is hit on 1-4 VDW line marker = marker + 1 # move marker along try: vdwof.append(float(words[3])) # VDW One-Four is 4th word after 1-4, VDW, and = eelof.append(float(words[7])) # EEL One-Four is 8th word (ignore restraint) except ValueError: # catch ****** in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 elif marker == 3: # this is hit on the ECAVITY/ENPOLAR line marker = marker + 1 # move marker along try: ecavity.append(float(words[2])) # ecavity/enpolar is the 3rd word after ECAVITY and = if apbs != 1: edisper.append(float(words[5])) # if we're not using APBS, get EDISPER except ValueError: # catc ******* in mdout print >> sys.stderr, 'Error: Sander output is missing values!' # print error message print >> sys.stderr, outputline.strip() # print offending line and return error code return -1 if len(bond) == 0: # make sure there's at least 1 potential term. All should be same length print >> sys.stderr, 'Error: No potential terms in sander output! Check output files.' return -1 # return error code after printing error message for x in range(len(bond)): # add up all terms to get the total try: if apbs != 1: sum = vdwaals[x] + eel[x] + egb[x] + vdwof[x] + eelof[x] + ecavity[x] + edisper[x] + \ bond[x] + angle[x] + dihed[x] else: sum = vdwaals[x] + eel[x] + egb[x] + vdwof[x] + eelof[x] + ecavity[x] + bond[x] + \ angle[x] + dihed[x] except IndexError: # this would be a very strange error to get... print >> sys.stderr, 'Error: Inconsistent number of potential terms in mdouts. This is a strange error' return -1 # return error code if this error is hit... total.append(sum) # tack the sum onto the end of the total array. return 1 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def round(number, decimals): import math try: num = float(number) # make sure the object passed is a float except ValueError: return -1 # else return error code try: dec = int(decimals) # make sure the second object is an integer except ValueError: return -1 # else return error code # move the desired decimal to the ones place and subtract off everything so only decimals are left test = (num * 10 ** decimals) - math.floor(num * 10 ** decimals) test = math.floor(test * 10) # first decimal we're NOT interested in so we can test to round up or down if test < 5: # if less than 5, round down (floor) toreturn = math.floor(num * 10 ** decimals) / 10 ** decimals else: # otherwise, round up (ceil) toreturn = math.ceil(num * 10 ** decimals) / 10 ** decimals return toreturn # return it #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def average(list): sum = 0 for x in range(len(list)): try: sum = sum + float(list[x]) except ValueError: return -1 if len(list) == 0: return -1 return sum / len(list) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def stdev(list, avg): import math if len(list) < 2: return 0 try: ave = float(avg) except ValueError: ave = average(list) sum = 0 for x in range(len(list)): sum = sum + (list[x] - ave) ** 2 return math.sqrt(sum / len(list)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def printusage(): import sys print "Usage: MMPBSA.py {-O} -i inputfile \\" print " -o output_file \\" print " -sp solvated_prmtop \\" print " -cp complex_prmtop \\" print " -rp receptor_prmtop \\" print " -lp ligand_prmtop \\" print " -y mdcrd1 mdcrd2 mdcrd3 ... mdcrdN \\" print " {-do decomp_output_file} \\" print " {-yr receptor_mdcrd1 ... receptor_mdcrdN} \\" print " {-yl ligand_mdcrd1 ... ligand_mdcrdN} \\ " print " {-mc mutant_complex_prmtop} \\" print " {-mr mutant_receptor_prmtop} \\" print " {-ml mutant_ligand_prmtop} \\" print " {-srp solvated_receptor_prmtop} \\" print " {-slp solvated_ligand_prmtop} \\" print " {-make-mdins} || {-use-mdins} || {-rewrite-output}" #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def digit(num, powten): import math tmp = math.floor(float(num) / 10 ** powten) # move desired digit to the ones place return int(tmp % 10) # return modulus of 10 so only ones place is returned, which is desired digit #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def CheckIncomp(INPUT, num): import sys, os # this does not have to be tested/exceptions caught because that was done in MMPBSA.py before this function ligand_prmtop = INPUT[4] startframe = int(INPUT[5]) endframe = int(INPUT[6]) interval = int(INPUT[7]) igb = int(INPUT[9]) gbsa = int(INPUT[10]) saltcon = float(INPUT[11]) extdiel = float(INPUT[12]) intdiel = float(INPUT[13]) surften = float(INPUT[14]) surfoff = float(INPUT[15]) indi = float(INPUT[16]) exdi = float(INPUT[17]) scale = float(INPUT[18]) linit = float(INPUT[19]) prbrad = float(INPUT[20]) istrng = float(INPUT[21]) npopt = int(INPUT[22]) cavity_surften = float(INPUT[23]) cavity_offset = float(INPUT[24]) fill_ratio = float(INPUT[25]) nproc = int(INPUT[26]) mpi_cmd = INPUT[26] radiopt = float(INPUT[28]) verbose = int(INPUT[29]) initial_traj = int(INPUT[30]) mutant_residue = int(INPUT[31]) ligand_mask = INPUT[36] sander_apbs = int(INPUT[37]) entropy = int(INPUT[39]) dielc = float(INPUT[40]) drms = float(INPUT[41]) maxcyc = int(INPUT[42]) nmstartframe = int(INPUT[43]) nmendframe = int(INPUT[44]) nminterval = int(INPUT[45]) idecomp = int(INPUT[54]) nmode_igb = int(INPUT[58]) nmode_istrng = float(INPUT[59]) inerr = 0 # return value. takes error code if hard error is encountered if startframe > endframe: print >> sys.stderr, 'Input Error: Your STARTFRAME, {0} must be lower than your ENDFRAME {1}'.\ format(startframe,endframe) inerr = -1 if startframe < 1 or endframe < 1: print >> sys.stderr, 'Input Error: STARTFRAME and ENDFRAME must both be positive!' inerr = -1 if igb != 1 and igb != 2 and igb != 5 and igb != 7: print >> sys.stderr, 'Input Error: Invalid value for IGB ({0}). It must be 1, 2, 5, or 7'.format(igb) inerr = -1 if saltcon < 0: print >> sys.stderr, 'Input Error: SALTCON must be non-negative!' inerr = -1 if extdiel < 0: print >> sys.stderr, 'Input Error: EXTDIEL must be non-negative!' inerr = -1 if intdiel < 0: print >> sys.stderr, 'Input Error: INTDIEL must be non-negative!' inerr = -1 if surften < 0: print >> sys.stderr, 'Input Error: SURFTEN must be non-negative!' inerr = -1 if indi < 0: print >> sys.stderr, 'Input Error: INDI must be non-negative!' inerr = -1 if exdi < 0: print >> sys.stderr, 'Input Error: EXDI must be non-negative!' inerr = -1 if scale <= 0: print >> sys.stderr, 'Input Error: SCALE must be postive!' inerr = -1 if linit < 1: print >> sys.stderr, 'Input Error: LINIT must be positive!' inerr = -1 if prbrad != 1.4 and prbrad != 1.6: print >> sys.stderr, 'Input Error: PRBRAD ({0}) only has support for values 1.4 and 1.6'.format(prbrad) inerr = -1 if istrng < 0: print >> sys.stderr, 'Input Error: ISTRNG must be non-negative!' inerr = -1 if npopt != 0 and npopt != 1 and npopt != 2: print >> sys.stderr, 'Input Error: INP ({0}) must be either 0, 1, or 2!'.format(npopt) inerr = -1 if cavity_surften < 0: print >> sys.stderr, 'Input Error: CAVITY_SURFTEN must be non-negative!' inerr = -1 if fill_ratio <= 0: print >> sys.stderr, 'Input Error: FILL_RATIO must be positive!' inerr = -1 if radiopt != 0 and radiopt != 1: print >> sys.stderr, 'Input Error: RADIOPT ({0}) must be 0 or 1!'.format(radiopt) inerr = -1 if dielc < 0: print >> sys.stderr, 'Input Error: dielc must be non-negative!' inerr = -1 if maxcyc < 1: print >> sys.stderr, 'Input Error: maxcyc ({0}) must be a positive integer!'.format(maxcyc) inerr = -1 if nmstartframe < 1 or nmendframe < 1 or nmendframe < nmstartframe: print >> sys.stderr, 'Input Error: nmstartframe and nmendframe must both be positive integers, ' print >> sys.stderr, ' and the first cannot be larger!' inerr = -1 if nminterval < 1: print >> sys.stderr, 'Input Error: nminterval ({0}) must be a positive integer!'.format(nminterval) inerr = -1 if idecomp < 0 or idecomp > 4: print >> sys.stderr, 'Input Error: IDECOMP must be 1, 2, 3, or 4!' inerr = -1 if idecomp != 0 and gbsa != 2: print >> sys.stderr, 'Input Error: When IDECOMP is set, GBSA must be set to 2!' inerr = -1 if idecomp != 0 and sander_apbs != 0: print >> sys.stderr, 'Input Error: IDECOMP cannot be used with sander.APBS!' inerr = -1 if entropy != 0 and entropy != 1: print >> sys.stderr, 'Input Error: ENTROPY ({0}) must be 0 or 1!'.format(entropy) inerr = -1 return inerr #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def CheckRadiis(igb, prmtops, warnings): import sys firstradii = getbondiset(prmtops[0]) # all radii to be checked must have the same radii set. Check this now if igb == 7 and firstradii != 'Bondi radii (bondi)': # this causes error in sander, so check it here print >> sys.stderr, 'Error: igb=7 requires the bondi radius set! Remake your topology file(s) with the command' print >> sys.stderr, ' set default pbradii bondi' return -1 # fatal error for x in range(1,len(prmtops)): if firstradii != getbondiset(prmtops[x]): # if there is a radii set that is not the same as the first... print >> sys.stderr, 'Error: Your topology files have different pbradii sets!' # print error message and... return -1 # return fatal error if firstradii == 'NO RADII SET SPECIFIED': # if prmtop doesn't have radius set, as in some versions of sleap print >> sys.stderr, 'Warning: No radii flag was specified in any of your topology files!' warnings.append('Warning: No radii set was specified in any of your topology files!') if igb == 1 and firstradii != 'amber6 modified Bondi radii (amber6)' and firstradii != 'modified Bondi radii (mbondi)': print >> sys.stderr, 'Warning: igb=1 should be used with either amber6 or mbondi pbradii set. Yours are ' + firstradii warnings.append('igb=1 should be used with either amber6 or mbondi pbradii set. Yours are ' + firstradii) elif igb == 2 and firstradii != 'H(N)-modified Bondi radii (mbondi2)': print >> sys.stderr, 'Warning: igb=2 should be used with mbondi2 pbradii set. Yours are ' + firstradii warnings.append('igb=2 should be used with mbondi2 pbradii set. Yours are ' + firstradii) elif igb==5 and firstradii != 'H(N)-modified Bondi radii (mbondi2)' \ and firstradii != 'Bondi radii (bondi)': print >> sys.stderr, 'Warning: igb=5 should be used with either mbondi2 or bondi pbradii set. Yours are ' + firstradii warnings.append('igb=5 should be used with either mbondi2 or bondi pbradii set. Yours are ' + firstradii) return 1 # all warnings are non-fatal #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getresinfo(res, topname, flag): infos = getallresinfo(topname, flag) # get all residue info return infos[res-1] # and simply return the residue of interest #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getbondiset(topname): import inputparse, sys try: topfile = open(topname, 'r') except IOError: print >> sys.stderr, 'Error: Cannot open prmtop {0}!'.format(topname) return -1 toplines = topfile.readlines() # load topology lines into memory topfile.close() # close topology file for x in range(len(toplines)): if '%FLAG RADIUS_SET' in toplines[x]: # look for this flag x = x + 1 # skip over the line while not toplines[x].startswith('%FORMAT'): x = x + 1 # skip over any comments return inputparse.Check(toplines[x+1]) # return the line describing radii set. return 'NO RADII SET SPECIFIED' # if it couldn't find radii set... #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getnumatms(resname): import sys # return number of atoms in specified residue if resname in 'GLY': return 7 if resname in 'ALA CYM CYX': return 10 if resname in 'CYS SER': return 11 if resname in 'ASP': return 12 if resname in 'ASN PRO THR': return 14 if resname in 'GLU': return 15 if resname in 'GLH VAL': return 16 if resname in 'GLN HID HIE MET': return 17 if resname in 'HIP': return 18 if resname in 'ILE LEU': return 19 if resname in 'PHE': return 20 if resname in 'LYN TYR': return 21 if resname in 'LYS': return 22 if resname in 'ARG TRP': return 24 else: print >> sys.stderr, 'Error: Unrecognized residue. Add {0} to getnumatms(resname).'.format(resname) return -1 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getressymbol(resname): # return 1-letter symbol of given amino acid if resname in 'ALA': return 'A' if resname in 'ARG': return 'R' if resname in 'ASN': return 'N' if resname in 'ASP ASH': return 'D' if resname in 'CYS CYX CYM': return 'C' if resname in 'GLU GLH': return 'E' if resname in 'GLN': return 'Q' if resname in 'GLY': return 'G' if resname in 'HIP HID HIE': return 'H' if resname in 'ILE': return 'I' if resname in 'LEU': return 'L' if resname in 'LYN LYS': return 'K' if resname in 'MET': return 'M' if resname in 'PHE': return 'F' if resname in 'PRO': return 'P' if resname in 'SER': return 'S' if resname in 'THR': return 'T' if resname in 'TRP': return 'W' if resname in 'TYR': return 'Y' if resname in 'VAL': return 'V' return 'ONE LETTER SYMBOL FOR {0} NOT FOUND'.format(resname) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def Format(num, uno, dos): number = str(round(num,dos)) # round number to given number of decimal places (dos) decimals = number.split('.') # split along decimal while len(decimals[0]) < uno: # pad number with spaces until it reaches desired length decimals[0] = ' ' + decimals[0] while len(decimals[1]) < dos: decimals[1] = decimals[1] + '0' # add zeroes until it reaches desired number of decimals return decimals[0] + '.' + decimals[1] # return result (length uno + dos + 1) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def ScaleDistance(coords, dist): import math if len(coords) != 6: # needs x1,y1,z1, x2,y2,z2. If not 6 numbers, return error code return -1 coords[3] = coords[3] - coords[0] # transform coordinates so that x1,y1,z1 are at the origin coords[4] = coords[4] - coords[1] coords[5] = coords[5] - coords[2] # the actual distance then becomes sqrt of x2**2 + y2**2 + z2**2 actualdist = math.sqrt(coords[3] * coords[3] + coords[4] * coords[4] + coords[5] * coords[5]) scalefactor = dist / actualdist # scale it so it reaches the desired distance coords[3] = scalefactor * coords[3] # scale all factors by desired amount coords[4] = scalefactor * coords[4] coords[5] = scalefactor * coords[5] coords[3] = coords[3] + coords[0] # move coordinates back so x1,y1,z1 are back to their original place coords[4] = coords[4] + coords[1] coords[5] = coords[5] + coords[2] return coords # return the coordinates #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getallresinfo(prmtop, flag): import sys # open topology file, read into memory, and close it top = open(prmtop,'r') toplines = top.readlines() top.close() items = [] # this is the array that holds all pieces of information gathered from a specific part of top file for x in range(len(toplines)): if flag in toplines[x]: # if flag is found in the topology line... x = x + 1 # skip to next line while not 'FORMAT' in toplines[x]: # skip over comments x = x + 1 x = x + 1 # skip over FORMAT line while not 'FLAG' in toplines[x]: # read until %FLAG is reached (next prmtop section) words = toplines[x].split() # split along whitespace (will not work for every block) for y in range(len(words)): items.append(words[y]) # add on each item in order x = x + 1 break # jump out of for loop return items # return array with all the information #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def CheckMutParm(prmtop1, prmtop2): import sys residues1 = [] # holder for residue names of prmtop1 residues2 = [] # holder for residue names of prmtop2 diffs = 0 # count the number of differences between residues1 and residues2 mutres = -1 # holds the location of the mutant residue; starts as error code residues1 = getallresinfo(prmtop1,'RESIDUE_LABEL') # load the residue names residues2 = getallresinfo(prmtop2,'RESIDUE_LABEL') if len(residues1) != len(residues2): # the residue numbers must be identical print >> sys.stderr, 'Error: {0} and {1} differ in their number of residues!'.format(prmtop1, prmtop2) return -1 # return error code for x in range(len(residues1)): # loop through all residues if residues1[x] != residues2[x]: diffs += 1 # increase number of differences if residues2[x] != 'ALA': # the 2nd residue must be ALA for alanine scanning print >> sys.stderr, 'Error: Mutant residue {0} is {1} but must be ALA!'.format(x+1,residues2[x]) return -1 # if not, print error message and return error code mutres = x + 1 # we've located the mutant residue number if diffs == 0: # if there are no differences, print error message print >> sys.stderr, 'Error: Your mutant complex topolgy file is the same as your original!' elif diffs > 1: # if there is more than one difference, print error message print >> sys.stderr, 'Error: Your mutant topology file {0} differs from {1} in more than one place!'.\ format(prmtop2, prmtop1) mutres = -1 # change mutres back to the error code return mutres # return the location of the mutant residue #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def CheckRecLigParms(recprm1, recprm2, ligprm1, ligprm2): import sys recresidues1 = [] # holds names of receptor residues for prmtop1 recresidues2 = [] # holds names of receptor residues for prmtop2 ligresidues1 = [] # holds names of ligand residues for prmtop1 ligresidues2 = [] # holds names of ligand residues for prmtop2 diffs = 0 # count the number of differences if recprm1 != recprm2: # only do this if the mutation is in the receptor recresidues1 = getallresinfo(recprm1,'RESIDUE_LABEL') # load receptor residues for prmtop 1 recresidues2 = getallresinfo(recprm2,'RESIDUE_LABEL') # load receptor residues for prmtop 2 if len(recresidues1) != len(recresidues2): # test that they're the same length print >> sys.stderr, 'Error: Mutant receptor prmtop does not have the same number of residues' print >> sys.stderr, ' as normal receptor prmtop!' return -1 # print error message and return error code if not for x in range(len(recresidues1)): # loop through all residues if recresidues1[x] != recresidues2[x]: # if they're not equal... diffs += 1 # this is a difference (count it) if recresidues2[x] != 'ALA': # make sure the second is an alanine print >> sys.stderr, 'Error: Receptor residue {0} differs, but mutant isn\'t ALA!'.format(x+1) return -1 # quit if it's not if ligprm1 != ligprm2: # if the mutation is in the ligand, do the same thing as for the receptor above ligresidues1 = getallresinfo(ligprm1,'RESIDUE_LABEL') ligresidues2 = getallresinfo(ligprm2,'RESIDUE_LABEL') if len(ligresidues1) != len(ligresidues2): print >> sys.stderr, 'Error: Mutant ligand prmtop does not have the same number of residues' print >> sys.stderr, ' as normal ligand prmtop!' return -1 for x in range(len(ligresidues1)): if ligresidues1[x] != ligresidues2[x]: diffs = diffs + 1 if ligresidues2[x] != 'ALA': print >> sys.stderr, 'Error: Ligand residue {0} differs, but mutant isn\'t ALA!'.format(x+1) return -1 if diffs == 0: # check for no differences print >> sys.stderr, 'Error: No difference between mutant/regular topology files!' return -1 elif diffs > 1: # check for multiple differences print >> sys.stderr, 'Error: More than one mutated residue!' return -1 return 1 # return success if it makes it this far #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getligmask(comprm, recprm, ligprm): import sys # returns [ligand position, ligand mask] comres = getallresinfo(comprm, 'RESIDUE_LABEL') # fill comres with complex residue names recres = getallresinfo(recprm, 'RESIDUE_LABEL') # fill recres with receptor residue names ligres = getallresinfo(ligprm, 'RESIDUE_LABEL') # fill ligres with ligand residue names cresnum = len(comres) # this is the number of residues in the complex rresnum = len(recres) # this is the number of residues in the receptor lresnum = len(ligres) # this is the number of residues in the ligand if cresnum != rresnum + lresnum: # complex residue number must equal sum of ligand and recptor resnums return -1 # if not, return error code cres = '' # residue from the Complex to compare ores = '' # residue from the Other prmtop to compare (rec or lig based on where in the complex it is) # The way this is supposed to work: start from the end of the receptor, and insert all ligand residues # at each location in the receptor residue sequence. Compare that resulting sequence to the complex # sequence. If it matches, we know the ligand_mask and receptor_mask. for x in range(rresnum): # loop through all positions in the receptor prmtop z = rresnum - x - 1 # z is the location of the ligand -- start from the end rather than the beginning rightpos = True # we will assume it's the RIGHT POSition until proven otherwise for y in range(cresnum): # now loop through all complex resnames to compare with rec/lig resnames cres = comres[y] # the residue name from the complex... if y <= z: # if our complex resname counter is less than our ligand location... ores = recres[y] # ores is from receptor at location y elif y > z and y - z <= lresnum: # otherwise, if it's greater than z and less than z + num residues in ligand... ores = ligres[y-z-1] # residue is in ligand, but subtract out number of receptor residues that aren't in ligand else: # otherwise, it's still in receptor, but after the ligand appears ores = recres[y-lresnum] # it's in receptor, but subtract number of ligand residues as they're not in receptor if cres != ores: # if they're not equal... rightpos = False # it's not the right position if rightpos: # we've now finished comparison for that particular ligand location, so if it's right position... if lresnum == 1: # if the ligand is only one residue: return [z+2,':{0}'.format(z+2)] # ligand residue is z+2 (+1 for index starting at 0, +1 for lig starting AFTER rec ends) else: # if ligand is more than 1 residue, we need to specify a range return [z+2,':{0}-{1}'.format(z+2,z+lresnum+1)] # same ligand start, but goes through lresnum (+1 for index at 0) # the previous loop will not check for ligand in position 1, do that now if we reach this far rightpos = True for x in range(cresnum): cres = comres[x] if x < lresnum: ores = ligres[x] else: ores = recres[x-lresnum] if cres != ores: rightpos = False if rightpos: if lresnum == 1: return [1,':1'] else: return [1,':1-{0}'.format(lresnum)] else: # if none of these work, we couldn't find the correct masks: return -1 # return error code #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def getrecmask(rresnum, ligpos, lresnum): # get receptor mask after we've gotten ligand mask and know where ligand starts recmask = ':' # receptor mask always starts with : for specifying residue if ligpos != 1: # take the case when ligand does not start in the first location if ligpos == 2: # if ligand is in position 2, receptor does not have a range for first location... recmask += '1' # it's just 1 else: # otherwise it does have a range recmask += '1-{0}'.format(ligpos-1) # it spans 1 through 1 position before the ligand starts else: # if ligpos is 1, receptor is the range after ligand ends recmask += '{0}-{1}'.format(1+lresnum,rresnum+lresnum) return recmask # now we know this is the end of the receptor mask -- return it if ligpos == rresnum + 1: # if the ligand position is 1 after receptor res number, it's at the end... return recmask # so return the receptor mask now else: # otherwise, the ligand is in the middle, and we have to add the second chunk of the receptor mask recmask += ':{0}'.format(ligpos+lresnum) if ligpos+lresnum != lresnum+rresnum: # if there's more than just 1 residue in receptor after ligand... recmask += '-{0}'.format(lresnum+rresnum) # specify the range to the end return recmask # return the final receptor mask #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def GetMasks(comprm, recprm, ligprm, debug): import sys # getligmask returns array with [lig starting position, lig mask] liginfo = getligmask(comprm, recprm, ligprm) # if it returned an error code if liginfo == -1: print >> sys.stderr, 'Error: Could not predict mask from topology files! Make sure ligand residues are' print >> sys.stderr, ' sequential or specify receptor_mask and ligand_mask in the input file.' remove(debug,0) # getrecmask returns string with receptor mask recinfo = getrecmask(resnum(recprm), liginfo[0], resnum(ligprm)) # return array with [ligand mask, rec mask, lig starting position] return [liginfo[1], recinfo, liginfo[0]] #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def maskToDec(mask): # split a mask at each : maskcomponents = mask.split(':') # array of :-delimited fields totstring = '' # total string for x in range(len(maskcomponents)): # loop through all fields maskcomponents[x] = maskcomponents[x].strip() # get rid of whitespace if len(maskcomponents[x]) == 0: # skip empty fields continue if '-' in maskcomponents[x]: # if multiple residues specified maskcomponents[x] = maskcomponents[x].replace('-',' ') # replace - with space else: # if single residue specified maskcomponents[x] += ' ' + maskcomponents[x] # duplicate residue as needed totstring += maskcomponents[x] + ' ' # add this to the total string return totstring # return the total string #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def PtrajNumFrames(filename): try: file = open(filename, 'r') # try to open the file except IOError: return -1 # return error code if we can't for line in file: # parse each line if 'Coordinate processing' in line: # this is the line that has number of frames run words = line.split() return words[5] # it is the 6th word in the line return -1 # return error code if we couldn't find the number of threads processed #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def nmodecalc(prefix, istrng, maxcyc, drms, gasorvac, dielc, complex_prmtop, receptor_prmtop, ligand_prmtop, suffix, exe, sfx, temp): import os,sys, math kappa = math.sqrt(0.10806*istrng) # calculate kappa from the ionic strength if gasorvac == 1: # do nmode calc in GB implicit solvent string = '"ntpr=10000, diel=C, kappa={0}, cut=1000, gb=1, dielc={1}, temp0={2}"'.format(kappa,dielc,temp) elif gasorvac == 0: # do nmode calc in gas phase string = '"ntpr=10000, diel=R, kappa={0}, cut=1000, gb=1, dielc={1}, temp0={2}"'.format(kappa,dielc, temp) print >> sys.stdout, ' master thread is calculating {0} frames\n'.format(len(suffix)) for x in range(len(suffix)): # loop through all frames to be done by this processor and calculate com, rec, and lig entropies print >> sys.stdout, ' calculating complex for frame number {0}'.format(x) os.system('{0} {1}complex_nm.pdb.{2} {3} {4} {5} {6} >> {1}complex_nm.out{7}'.\ format(exe,prefix,suffix[x],complex_prmtop,maxcyc,drms,string,sfx)) print >> sys.stdout, ' calculating receptor for frame number {0}'.format(x) os.system('{0} {1}receptor_nm.pdb.{2} {3} {4} {5} {6} >> {1}receptor_nm.out{7}'.\ format(exe,prefix,suffix[x],receptor_prmtop,maxcyc,drms,string,sfx)) print >> sys.stdout, ' calculating ligand for frame number {0}'.format(x) os.system('{0} {1}ligand_nm.pdb.{2} {3} {4} {5} {6} >> {1}ligand_nm.out{7}'.\ format(exe,prefix,suffix[x],ligand_prmtop,maxcyc,drms,string,sfx)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def gbcalc(sander, suffix, prefix, complex_prmtop, receptor_prmtop, ligand_prmtop, idecomp): import os, sys if idecomp == 0: # if we aren't doing decomp, all mdins are the same mdin1 = "_MMPBSA_gb.mdin" mdin2 = "_MMPBSA_gb.mdin" mdin3 = "_MMPBSA_gb.mdin" else: # otherwise each simulation has its own mdin mdin1 = "_MMPBSA_gb_decomp_com.mdin" mdin2 = "_MMPBSA_gb_decomp_rec.mdin" mdin3 = "_MMPBSA_gb_decomp_lig.mdin" # run the calculations print >> sys.stdout, ' calculating ligand contribution...' os.system(('{0} -O -i {1} -o {2}ligand_gb.mdout{3} -p {4} -c {2}' + 'dummyligand.inpcrd -y {2}ligand.mdcrd{3}').format(sander,mdin3, prefix, suffix, ligand_prmtop)) print >> sys.stdout, ' calculating receptor contribution...' os.system(('{0} -O -i {1} -o {2}receptor_gb.mdout{3} -p {4} -c {2}' + 'dummyreceptor.inpcrd -y {2}receptor.mdcrd{3}').format(sander,mdin2, prefix, suffix, receptor_prmtop)) print >> sys.stdout, ' calculating complex contribution...' os.system(('{0} -O -i {1} -o {2}complex_gb.mdout{3} -p {4} -c {2}' + 'dummycomplex.inpcrd -y {2}complex.mdcrd{3}').format(sander,mdin1, prefix, suffix, complex_prmtop)) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def pbcalc(sander, suffix, prefix, complex_prmtop, receptor_prmtop, ligand_prmtop, idecomp): import os, sys if idecomp == 0: # if idecomp=0, all mdins are the same mdin1 = "_MMPBSA_pb.mdin" mdin2 = "_MMPBSA_pb.mdin" mdin3 = "_MMPBSA_pb.mdin" else: # otherwise, each simulation has its own mdin mdin1 = "_MMPBSA_pb_decomp_com.mdin" mdin2 = "_MMPBSA_pb_decomp_rec.mdin" mdin3 = "_MMPBSA_pb_decomp_lig.mdin" # run the calculations cmd = ('{0} -O -i {1} -o {2}ligand_pb.mdout{3} -p {4} -c {2}' + 'dummyligand.inpcrd -y {2}ligand.mdcrd{3}').format(sander,mdin3, prefix, suffix, ligand_prmtop) if sander.endswith('APBS'): cmd += " >> _MMPBSA_pbsanderoutput.junk" # if we're doing APBS, redirect apbs stdoutput to another file print >> sys.stdout, ' calculating ligand contribution...' os.system(cmd) cmd = ('{0} -O -i {1} -o {2}receptor_pb.mdout{3} -p {4} -c {2}' + 'dummyreceptor.inpcrd -y {2}receptor.mdcrd{3}').format(sander,mdin2, prefix, suffix, receptor_prmtop) if sander.endswith('APBS'): cmd += " >> _MMPBSA_pbsanderoutput.junk" # if we're doing APBS, redirect apbs stdoutput to another file print >> sys.stdout, ' calculating receptor contribution...' os.system(cmd) cmd = ('{0} -O -i {1} -o {2}complex_pb.mdout{3} -p {4} -c {2}' + 'dummycomplex.inpcrd -y {2}complex.mdcrd{3}').format(sander,mdin1, prefix, suffix, complex_prmtop) if sander.endswith('APBS'): cmd += " >> _MMPBSA_pbsanderoutput.junk" # if we're doing APBS, redirect apbs stdoutput to another file print >> sys.stdout, ' calculating complex contribution...' os.system(cmd) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def decomp_out(finaloutput, idecomp, debug, verbose, type, frames, comprm, recprm, ligprm, ligstart, surften): import os # this module controls the writing of each type of decomp output outputfile = open(finaloutput,'a') # open the output file and append to it if idecomp == 1 or idecomp == 2: # non-pairwise, call non-pairwise parser return decoutput12(outputfile,debug,verbose,frames,comprm,recprm,ligprm,ligstart,type,idecomp,surften) elif idecomp == 3 or idecomp == 4: # pairwise: call pairwise parser return decoutput34(outputfile,debug,verbose,frames,comprm,recprm,ligprm,ligstart,type,idecomp,surften,finaloutput) outputfile.close() # close the output file return -1 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def gboutput(mdout1, mdout2, mdout3, header, finaloutput, debug, numframes, onetraj, verbose): import sys, math bonddif = [] # array containing the differences in the BOND term for each frame angledif = [] # array containing the differences in the ANGLE term for each frame diheddif = [] # array containing the differences in the DIHED term for each frame vdwaalsdif = [] # array containing the differences in the VDWAALS term for each frame eeldif = [] # array containing the differences in the EEL term for each frame egbdif = [] # array containing the differences in the EGB term for each frame vdwofdif = [] # array containing the differences in the 1-4 VDW term for each frame eelofdif = [] # array containing the differences in the 1-4 EEL term for each frame esurfdif = [] # array containing the differences in the ESURF term for each frame deltag = [] # array containing the delta G for each frame bond = [] # array containing the BOND term for each frame: temp for each mdout file angle = [] # array containing the ANGLE term for each frame: temp for each mdout file dihed = [] # array containing the DIHED term for each frame: temp for each mdout file vdwaals = [] # array containing the VDWAALS term for each frame: temp for each mdout file eel = [] # array containing the EEL term for each frame: temp for each mdout file egb = [] # array containing the EGB term for each frame: temp for each mdout file vdwof = [] # array containing the 1-4 VDWAALS term for each frame: temp for each mdout file eelof = [] # array containing the 1-4 EEL term for each frame: temp for each mdout file esurf = [] # array containing the ESURF term for each frame: temp for each mdout file total = [] # array containing the TOTAL of all terms for each frame: temp for each mdout file # arrays to hold averages and std devs for complex, receptor, and ligand: they are arrays of 3-item # arrays that hold the value for the complex, receptor, and ligand, in that order. bondavg = [] # 3-item array that holds average BOND angleavg = [] # 3-item array that holds average ANGLE dihedavg = [] # 3-item array that holds average DIHED vdwaalsavg = [] # 3-item array that holds average VDWAALS eelavg = [] # 3-item array that holds average EEL egbavg = [] # 3-item array that holds average EGB vdwofavg = [] # 3-item array that holds average 1-4 VDWAALS eelofavg = [] # 3-item array that holds average 1-4 EEL esurfavg = [] # 3-item array that holds average ESURF totalavg = [] # 3-item array that holds average TOTAL sum of all terms bondstdev = [] # 3-item array that holds standard deviation BOND anglestdev = [] # 3-item array that holds standard deviation ANGLE dihedstdev = [] # 3-item array that holds standard deviation DIHED vdwaalsstdev = [] # 3-item array that holds standard deviation VDWAALS eelstdev = [] # 3-item array that holds standard deviation EEL egbstdev = [] # 3-item array that holds standard deviation EGB vdwofstdev = [] # 3-item array that holds standard deviation 1-4 VDWAALS eelofstdev = [] # 3-item array that holds standard deviation 1-4 EEL esurfstdev = [] # 3-item array that holds standard deviation ESURF totalstdev = [] # 3-item array that holds standard deviation TOTAL sum of all terms # get the energies and such from the complex try: outputfile = open(mdout1,'r') # try to open up the first mdout file (complex) except IOError: print >> sys.stderr, 'Error: Sander error! GB mdout file was not created.' remove(debug,0) # catch error message and quit if no mdout file is created # send the output file along with arrays to getenergies. The arrays come back filled isinerr = getenergies(bond, angle, dihed, vdwaals, eel, egb, vdwof, eelof, esurf, total, outputfile, onetraj) outputfile.close() # we have what we need from the file, close it if isinerr == -1: # if getenergies returned an error code, it already printed, so just quit remove(debug,0) # for the complex, append all totals to the deltag. Receptor and ligand will each be subtracted from this for x in range(len(total)): deltag.append(total[x]) for x in range(len(bond)): # the complex value is appended to each difference. rec/lig will be subtracted later bonddif.append(bond[x]) angledif.append(angle[x]) diheddif.append(dihed[x]) vdwaalsdif.append(vdwaals[x]) eeldif.append(eel[x]) egbdif.append(egb[x]) vdwofdif.append(vdwof[x]) eelofdif.append(eelof[x]) esurfdif.append(esurf[x]) # add the averages for the complex to each array (1st element) bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) esurfavg.append(average(esurf)) totalavg.append(average(total)) # add the standard deviations for the complex to each array (1st element) bondstdev.append(stdev(bond, bondavg[0])) anglestdev.append(stdev(angle, angleavg[0])) dihedstdev.append(stdev(dihed, dihedavg[0])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[0])) eelstdev.append(stdev(eel, eelavg[0])) egbstdev.append(stdev(egb, egbavg[0])) vdwofstdev.append(stdev(vdwof, vdwofavg[0])) eelofstdev.append(stdev(eelof, eelofavg[0])) esurfstdev.append(stdev(esurf, esurfavg[0])) totalstdev.append(stdev(total, totalavg[0])) # clear the energy arrays and move to the receptor bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] esurf = [] total = [] # now try and do the same for receptor: I will only highlight differences here. For full comments, see section for # for complex above try: outputfile = open(mdout2,'r') except IOError: print >> sys.stderr, 'Error: Sander error! GB mdout file was not created.' remove(debug,0) isinerr = getenergies(bond, angle, dihed, vdwaals, eel, egb, vdwof, eelof, esurf, total, outputfile, onetraj) outputfile.close() if isinerr == -1: remove(debug,0) for x in range(len(deltag)): deltag[x] = deltag[x] - total[x] for x in range(len(bond)): # subtract values of receptor from those already present from complex bonddif[x] = bonddif[x] - bond[x] angledif[x] = angledif[x] - angle[x] diheddif[x] = diheddif[x] - dihed[x] vdwaalsdif[x] = vdwaalsdif[x] - vdwaals[x] eeldif[x] = eeldif[x] - eel[x] egbdif[x] = egbdif[x] - egb[x] vdwofdif[x] = vdwofdif[x] - vdwof[x] eelofdif[x] = eelofdif[x] - eelof[x] esurfdif[x] = esurfdif[x] - esurf[x] # set the averages of the receptor bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) esurfavg.append(average(esurf)) totalavg.append(average(total)) # set the standard deviations of the receptor bondstdev.append(stdev(bond, bondavg[1])) anglestdev.append(stdev(angle, angleavg[1])) dihedstdev.append(stdev(dihed, dihedavg[1])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[1])) eelstdev.append(stdev(eel, eelavg[1])) egbstdev.append(stdev(egb, egbavg[1])) vdwofstdev.append(stdev(vdwof, vdwofavg[1])) eelofstdev.append(stdev(eelof, eelofavg[1])) esurfstdev.append(stdev(esurf, esurfavg[1])) totalstdev.append(stdev(total, totalavg[1])) # clear the energy arrays and move to the ligand bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] esurf = [] total = [] try: outputfile = open(mdout3,'r') except IOError: print >> sys.stderr, 'Error: Sander error! GB mdout file was not created.' remove(debug,0) isinerr = getenergies(bond, angle, dihed, vdwaals, eel, egb, vdwof, eelof, esurf, total, outputfile, onetraj) outputfile.close() if isinerr == -1: remove(debug,0) for x in range(len(deltag)): deltag[x] = deltag[x] - total[x] for x in range(len(bond)): # subtract ligand values from the differences (already complex-receptor) bonddif[x] = bonddif[x] - bond[x] angledif[x] = angledif[x] - angle[x] diheddif[x] = diheddif[x] - dihed[x] vdwaalsdif[x] = vdwaalsdif[x] - vdwaals[x] eeldif[x] = eeldif[x] - eel[x] egbdif[x] = egbdif[x] - egb[x] vdwofdif[x] = vdwofdif[x] - vdwof[x] eelofdif[x] = eelofdif[x] - eelof[x] esurfdif[x] = esurfdif[x] - esurf[x] # differences are now complete # set the averages of the ligand bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) esurfavg.append(average(esurf)) totalavg.append(average(total)) deltagavg = average(deltag) # set the standard deviations of the ligand bondstdev.append(stdev(bond, bondavg[2])) anglestdev.append(stdev(angle, angleavg[2])) dihedstdev.append(stdev(dihed, dihedavg[2])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[2])) eelstdev.append(stdev(eel, eelavg[2])) egbstdev.append(stdev(egb, egbavg[2])) vdwofstdev.append(stdev(vdwof, vdwofavg[2])) eelofstdev.append(stdev(eelof, eelofavg[2])) esurfstdev.append(stdev(esurf, esurfavg[2])) totalstdev.append(stdev(total, totalavg[2])) deltagstdev = stdev(deltag, deltagavg) # clear unneeded arrays to free up memory bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] esurf = [] total = [] # get the averages of the differences bonddifavg = average(bonddif) # average of BOND differences angledifavg = average(angledif) # average of ANGLE differences diheddifavg = average(diheddif) # average of DIHED differences vdwaalsdifavg = average(vdwaalsdif) # average of VDWAALS differences eeldifavg = average(eeldif) # average of EEL differences egbdifavg = average(egbdif) # average of EGB differences vdwofdifavg = average(vdwofdif) # average of 1-4 VDW differences eelofdifavg = average(eelofdif) # average of 1-4 EEL differences esurfdifavg = average(esurfdif) # average of ESURF differences bonddifstdev = stdev(bonddif, bonddifavg) # standard deviation of BOND differences angledifstdev = stdev(angledif, angledifavg) # standard deviation of ANGLE differences diheddifstdev = stdev(diheddif, diheddifavg) # standard deviation of DIHED differences vdwaalsdifstdev = stdev(vdwaalsdif, vdwaalsdifavg) # standard deviation of VDWAALS differences eeldifstdev = stdev(eeldif, eeldifavg) # standard deviation of EEL differences egbdifstdev = stdev(egbdif, egbdifavg) # standard deviation of EGB differences vdwofdifstdev = stdev(vdwofdif, vdwofdifavg) # standard deviation of 1-4 VDW differences eelofdifstdev = stdev(eelofdif, eelofdifavg) # standard deviation of 1-4 EEL differences esurfdifstdev = stdev(esurfdif, esurfdifavg) # standard deviation of ESURF differences gascomavg = vdwaalsavg[0] + eelavg[0] # vdwaals and eel are always part of DELTA G gas (for complex) gasrecavg = vdwaalsavg[1] + eelavg[1] # vdwaals and eel are always part of DELTA G gas (for receptor) gasligavg = vdwaalsavg[2] + eelavg[2] # vdwaals and eel are always part of DELTA G gas (for ligand) gasdifavg = vdwaalsdifavg + eeldifavg # vdwaals and eel are always part of DELTA G gas (for differences) gascomstd = vdwaalsstdev[0] ** 2 + eelstdev[0] ** 2 # propogate standard deviation (don't take sqrt yet...) gasrecstd = vdwaalsstdev[1] ** 2 + eelstdev[1] ** 2 # propogate standard deviation (don't take sqrt yet...) gasligstd = vdwaalsstdev[2] ** 2 + eelstdev[2] ** 2 # propogate standard deviation (don't take sqrt yet...) gasdifstd = vdwaalsdifstdev ** 2 + eeldifstdev ** 2 # propogate standard deviation (don't take sqrt yet...) solvcomavg = egbavg[0] + esurfavg[0] # DELTA G solv is egb + esurf (pol + npol) for complex solvrecavg = egbavg[1] + esurfavg[1] # for receptor solvligavg = egbavg[2] + esurfavg[2] # for ligand solvdifavg = egbdifavg + esurfdifavg # for differences solvcomstd = math.sqrt(egbstdev[0] ** 2 + esurfstdev[0] ** 2) # propagate stdevs. no other terms exist: take sqrt solvrecstd = math.sqrt(egbstdev[1] ** 2 + esurfstdev[1] ** 2) solvligstd = math.sqrt(egbstdev[2] ** 2 + esurfstdev[2] ** 2) solvdifstd = math.sqrt(egbdifstdev ** 2 + esurfdifstdev ** 2) if onetraj: # if we have only one trajectory, check the internal potential terms for consistency failed = False # it hasn't failed until we prove it has (bond, angle, dihed) failedof = False # it hasn't failed until we prove it has (1-4 interactions) for x in range(len(bonddif)): # for all differences, check that the diffs are not greater than 0.001 cutoff if (abs(bonddif[x]) > 0.001 or abs(angledif[x]) > 0.001 or abs(diheddif[x]) > 0.001) and not failed: finaloutput.write(""" WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS (BOND, ANGLE, AND/OR DIHED). CHECK YOUR INPUT FILES AND SYSTEM SETUP. THESE RESULTS MAY NOT BE RELIABLE (check differences)! """) failed = True # we've now shown it's failed elif (abs(vdwofdif[x]) > 0.001 or abs(eelofdif[x]) > 0.001) and not failedof: finaloutput.write(""" WARNING: INCONSISTENCIES EXIST WITHIN 1-4 NON-BONDED TERMS. CHECK YOUR INPUT FILES AND SYSTEM SETUP. THESE RESULTS MAY NOT BE RELIABLE (check differences)! The most common cause of this is inconsistent charge definitions across topology files. """) failedof = True # we've now shown that 1-4 interactions failed test if not onetraj or failed or verbose == 2: # add bond/angle/dihed to gas terms if it failed, it's asked for, or not onetraj gascomavg += bondavg[0] + angleavg[0] + dihedavg[0] gasrecavg += bondavg[1] + angleavg[1] + dihedavg[1] gasligavg += bondavg[2] + angleavg[2] + dihedavg[2] gasdifavg += bonddifavg + angledifavg + diheddifavg gascomstd += bondstdev[0]**2 + anglestdev[0]**2 + dihedstdev[0]**2 gasrecstd += bondstdev[1]**2 + anglestdev[1]**2 + dihedstdev[1]**2 gasligstd += bondstdev[2]**2 + anglestdev[2]**2 + dihedstdev[2]**2 gasdifstd += bonddifstdev**2 + angledifstdev**2 + diheddifstdev**2 else: # if we don't add them, then subtract them from the total totalavg[0] = totalavg[0]-bondavg[0]-angleavg[0]-dihedavg[0] totalavg[1] = totalavg[1]-bondavg[1]-angleavg[1]-dihedavg[1] totalavg[2] = totalavg[2]-bondavg[2]-angleavg[2]-dihedavg[2] if not onetraj or failedof or verbose == 2: # add 1-4 interactions to gas terms if it failed, it's asked for, or not onetraj gascomavg += eelofavg[0] + vdwofavg[0] gasrecavg += eelofavg[1] + vdwofavg[1] gasligavg += eelofavg[2] + vdwofavg[2] gasdifavg += eelofdifavg + vdwofdifavg gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 else: # if we don't add them, then subtract them from the total totalavg[0] = totalavg[0]-eelofavg[0]-vdwofavg[0] totalavg[1] = totalavg[1]-eelofavg[1]-vdwofavg[1] totalavg[2] = totalavg[2]-eelofavg[2]-vdwofavg[2] gascomstd = math.sqrt(gascomstd) # NOW take sqrt of stdevs since we're done adding stuff gasrecstd = math.sqrt(gasrecstd) gasligstd = math.sqrt(gasligstd) gasdifstd = math.sqrt(gasdifstd) finaloutput.write(header + '\nGENERALIZED BORN:\n\n') if verbose >= 1 or failed or failedof: # only print complex if asked for (verbose >= 1) or if we failed finaloutput.write('Complex:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: # only print out internals if we have multi trajectories or it failed or asked for finaloutput.write('BOND ' + Format(bondavg[0],15,4) + Format(bondstdev[0],17,4) + Format(bondstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[0],15,4) + Format(anglestdev[0],17,4) + Format(anglestdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[0],15,4) + Format(dihedstdev[0],17,4) + Format(dihedstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[0],15,4) + Format(vdwaalsstdev[0],17,4) + Format(vdwaalsstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[0],15,4) + Format(eelstdev[0],17,4) + Format(eelstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: # only print out 1-4s if we have multi trajectories or it failed or asked for finaloutput.write('1-4 VDW ' + Format(vdwofavg[0],15,4) + Format(vdwofstdev[0],17,4) + Format(vdwofstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[0],15,4) + Format(eelofstdev[0],17,4) + Format(eelofstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EGB ' + Format(egbavg[0],15,4) + Format(egbstdev[0],17,4) + Format(egbstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ESURF ' + Format(esurfavg[0],15,4) + Format(esurfstdev[0],17,4) + Format(esurfstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gascomavg,15,4) + Format(gascomstd,17,4) + Format(gascomstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvcomavg,15,4) + Format(solvcomstd,17,4) + Format(solvcomstd/ math.sqrt(numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[0],15,4) + Format(totalstdev[0],17,4) + Format(totalstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') finaloutput.write('Receptor:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: # only print out internals if failed, multi traj, or asked for finaloutput.write('BOND ' + Format(bondavg[1],15,4) + Format(bondstdev[1],17,4) + Format(bondstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[1],15,4) + Format(anglestdev[1],17,4) + Format(anglestdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[1],15,4) + Format(dihedstdev[1],17,4) + Format(dihedstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[1],15,4) + Format(vdwaalsstdev[1],17,4) + Format(vdwaalsstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[1],15,4) + Format(eelstdev[1],17,4) + Format(eelstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: # only print out 1-4s if failed, multi traj, or asked for finaloutput.write('1-4 VDW ' + Format(vdwofavg[1],15,4) + Format(vdwofstdev[1],17,4) + Format(vdwofstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[1],15,4) + Format(eelofstdev[1],17,4) + Format(eelofstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EGB ' + Format(egbavg[1],15,4) + Format(egbstdev[1],17,4) + Format(egbstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ESURF ' + Format(esurfavg[1],15,4) + Format(esurfstdev[1],17,4) + Format(esurfstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gasrecavg,15,4) + Format(gasrecstd,17,4) + Format(gasrecstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvrecavg,15,4) + Format(solvrecstd,17,4) + Format(solvrecstd/ math.sqrt(numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[1],15,4) + Format(totalstdev[1],17,4) + Format(totalstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') finaloutput.write('Ligand:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: # only print out internals if failed, multi traj, or asked for finaloutput.write('BOND ' + Format(bondavg[2],15,4) + Format(bondstdev[2],17,4) + Format(bondstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[2],15,4) + Format(anglestdev[2],17,4) + Format(anglestdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[2],15,4) + Format(dihedstdev[2],17,4) + Format(dihedstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[2],15,4) + Format(vdwaalsstdev[2],17,4) + Format(vdwaalsstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[2],15,4) + Format(eelstdev[2],17,4) + Format(eelstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: # only print out 1-4s if failed, multi traj or asked for finaloutput.write('1-4 VDW ' + Format(vdwofavg[2],15,4) + Format(vdwofstdev[2],17,4) + Format(vdwofstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[2],15,4) + Format(eelofstdev[2],17,4) + Format(eelofstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EGB ' + Format(egbavg[2],15,4) + Format(egbstdev[2],17,4) + Format(egbstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ESURF ' + Format(esurfavg[2],15,4) + Format(esurfstdev[2],17,4) + Format(esurfstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gasligavg,15,4) + Format(gasligstd,17,4) + Format(gasligstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvligavg,15,4) + Format(solvligstd,17,4) + Format(solvligstd/math.sqrt( numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[2],15,4) + Format(totalstdev[2],17,4) + Format(totalstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') # now print out differences finaloutput.write('Differences (Complex - Receptor - Ligand):') finaloutput.write('\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: # only print internals if multi traj, failed, or asked for finaloutput.write('BOND ' + Format(bonddifavg,15,4) + Format(bonddifstdev,17,4) + Format(bonddifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angledifavg,15,4) + Format(angledifstdev,17,4) + Format(angledifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(diheddifavg,15,4) + Format(diheddifstdev,17,4) + Format(diheddifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsdifavg,15,4) + Format(vdwaalsdifstdev,17,4) + Format(vdwaalsdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eeldifavg,15,4) + Format(eeldifstdev,17,4) + Format(eeldifstdev/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: # only print out 1-4s if multi traj, failed, or asked for finaloutput.write('1-4 VDW ' + Format(vdwofdifavg,15,4) + Format(vdwofdifstdev,17,4) + Format(vdwofdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofdifavg,15,4) + Format(eelofdifstdev,17,4) + Format(eelofdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EGB ' + Format(egbdifavg,15,4) + Format(egbdifstdev,17,4) + Format(egbdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ESURF ' + Format(esurfdifavg,15,4) + Format(esurfdifstdev,17,4) + Format(esurfdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nDELTA G gas ' + Format(gasdifavg,15,4) + Format(gasdifstd,17,4) + Format(gasdifstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('DELTA G solv ' + Format(solvdifavg,15,4) + Format(solvdifstd,17,4) + Format(solvdifstd/ math.sqrt(numframes),15,4) + '\n\n') finaloutput.write('\n DELTA G binding = ' + Format(deltagavg,10,4) + ' +/-' + Format(deltagstdev,7,4) + Format(deltagstdev/math.sqrt(float(numframes)),18,4) + '\n') return [deltagavg,deltagstdev] # return the average and stdev of DELTA G #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def pboutput(mdout1, mdout2, mdout3, header, finaloutput, debug, numframes, apbs, onetraj, verbose): import math, sys # the PB parser is virtually the same as the GB parser, only it has edisper and ecavity instead of esurf bonddif = [] angledif = [] diheddif = [] vdwaalsdif = [] eeldif = [] egbdif = [] vdwofdif = [] eelofdif = [] ecavitydif = [] edisperdif = [] deltag = [] bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] ecavity = [] edisper = [] total = [] # arrays to hold the avgs and std devs for complex, receptor, and ligand bondavg = [] angleavg = [] dihedavg = [] vdwaalsavg = [] eelavg = [] egbavg = [] vdwofavg = [] eelofavg = [] ecavityavg = [] edisperavg = [] totalavg = [] bondstdev = [] anglestdev = [] dihedstdev = [] vdwaalsstdev = [] eelstdev = [] egbstdev = [] vdwofstdev = [] eelofstdev = [] ecavitystdev = [] edisperstdev = [] totalstdev = [] # get the energies and such from the complex try: outputfile = open(mdout1,'r') except IOError: print >> sys.stderr, 'Error: Sander error! PB mdout file was not created.' remove(debug,0) isinerr = getpbenergies(bond,angle,dihed,vdwaals,eel,egb,vdwof,eelof,ecavity,edisper,total,outputfile,apbs,onetraj) outputfile.close() if isinerr == -1: remove(debug,0) for x in range(len(total)): deltag.append(total[x]) for x in range(len(bond)): bonddif.append(bond[x]) angledif.append(angle[x]) diheddif.append(dihed[x]) vdwaalsdif.append(vdwaals[x]) eeldif.append(eel[x]) egbdif.append(egb[x]) vdwofdif.append(vdwof[x]) eelofdif.append(eelof[x]) ecavitydif.append(ecavity[x]) if apbs != 1: edisperdif.append(edisper[x]) # set the averages of the complex bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) ecavityavg.append(average(ecavity)) if apbs != 1: edisperavg.append(average(edisper)) totalavg.append(average(total)) # set the standard deviations of the complex bondstdev.append(stdev(bond, bondavg[0])) anglestdev.append(stdev(angle, angleavg[0])) dihedstdev.append(stdev(dihed, dihedavg[0])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[0])) eelstdev.append(stdev(eel, eelavg[0])) egbstdev.append(stdev(egb, egbavg[0])) vdwofstdev.append(stdev(vdwof, vdwofavg[0])) eelofstdev.append(stdev(eelof, eelofavg[0])) ecavitystdev.append(stdev(ecavity, ecavityavg[0])) if apbs != 1: edisperstdev.append(stdev(edisper, edisperavg[0])) totalstdev.append(stdev(total, totalavg[0])) # clear the energy arrays and move to the receptor bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] ecavity = [] edisper = [] total = [] try: outputfile = open(mdout2,'r') except IOError: print >> sys.stderr, 'Error: Sander error! PB mdout file was ' + \ 'not created.' remove(debug,0) isinerr = getpbenergies(bond,angle,dihed,vdwaals,eel,egb,vdwof,eelof,ecavity,edisper,total,outputfile,apbs,onetraj) outputfile.close() if isinerr == -1: remove(debug,0) for x in range(len(deltag)): deltag[x] = deltag[x] - total[x] for x in range(len(bond)): bonddif[x] = bonddif[x] - bond[x] angledif[x] = angledif[x] - angle[x] diheddif[x] = diheddif[x] - dihed[x] vdwaalsdif[x] = vdwaalsdif[x] - vdwaals[x] eeldif[x] = eeldif[x] - eel[x] egbdif[x] = egbdif[x] - egb[x] vdwofdif[x] = vdwofdif[x] - vdwof[x] eelofdif[x] = eelofdif[x] - eelof[x] ecavitydif[x] = ecavitydif[x] - ecavity[x] if apbs != 1: edisperdif[x] = edisperdif[x] - edisper[x] # set the averages of the receptor bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) ecavityavg.append(average(ecavity)) if apbs != 1: edisperavg.append(average(edisper)) totalavg.append(average(total)) # set the standard deviations of the receptor bondstdev.append(stdev(bond, bondavg[1])) anglestdev.append(stdev(angle, angleavg[1])) dihedstdev.append(stdev(dihed, dihedavg[1])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[1])) eelstdev.append(stdev(eel, eelavg[1])) egbstdev.append(stdev(egb, egbavg[1])) vdwofstdev.append(stdev(vdwof, vdwofavg[1])) eelofstdev.append(stdev(eelof, eelofavg[1])) ecavitystdev.append(stdev(ecavity, ecavityavg[1])) if apbs != 1: edisperstdev.append(stdev(edisper, edisperavg[1])) totalstdev.append(stdev(total, totalavg[1])) # clear the energy arrays and move to the ligand bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] ecavity = [] edisper = [] total = [] try: outputfile = open(mdout3,'r') except IOError: print >> sys.stderr, 'Error: Sander error! PB mdout file was not created.' remove(debug,0) isinerr = getpbenergies(bond,angle,dihed,vdwaals,eel,egb,vdwof,eelof,ecavity,edisper,total,outputfile,apbs,onetraj) outputfile.close() if isinerr == -1: remove(debug,0) for x in range(len(deltag)): deltag[x] = deltag[x] - total[x] for x in range(len(bond)): bonddif[x] = bonddif[x] - bond[x] angledif[x] = angledif[x] - angle[x] diheddif[x] = diheddif[x] - dihed[x] vdwaalsdif[x] = vdwaalsdif[x] - vdwaals[x] eeldif[x] = eeldif[x] - eel[x] egbdif[x] = egbdif[x] - egb[x] vdwofdif[x] = vdwofdif[x] - vdwof[x] eelofdif[x] = eelofdif[x] - eelof[x] ecavitydif[x] = ecavitydif[x] - ecavity[x] if apbs != 1: edisperdif[x] = edisperdif[x] - edisper[x] # set the averages of the ligand bondavg.append(average(bond)) angleavg.append(average(angle)) dihedavg.append(average(dihed)) vdwaalsavg.append(average(vdwaals)) eelavg.append(average(eel)) egbavg.append(average(egb)) vdwofavg.append(average(vdwof)) eelofavg.append(average(eelof)) ecavityavg.append(average(ecavity)) if apbs != 1: edisperavg.append(average(edisper)) totalavg.append(average(total)) # set the standard deviations of the ligand bondstdev.append(stdev(bond, bondavg[2])) anglestdev.append(stdev(angle, angleavg[2])) dihedstdev.append(stdev(dihed, dihedavg[2])) vdwaalsstdev.append(stdev(vdwaals, vdwaalsavg[2])) eelstdev.append(stdev(eel, eelavg[2])) egbstdev.append(stdev(egb, egbavg[2])) vdwofstdev.append(stdev(vdwof, vdwofavg[2])) eelofstdev.append(stdev(eelof, eelofavg[2])) ecavitystdev.append(stdev(ecavity, ecavityavg[2])) if apbs != 1: edisperstdev.append(stdev(edisper, edisperavg[2])) totalstdev.append(stdev(total, totalavg[2])) deltagavg = average(deltag) deltagstdev = stdev(deltag,deltagavg) # clear unneeded arrays bond = [] angle = [] dihed = [] vdwaals = [] eel = [] egb = [] vdwof = [] eelof = [] ecavity = [] edisper = [] total = [] # get the averages and stdevs of the differences bonddifavg = average(bonddif) angledifavg = average(angledif) diheddifavg = average(diheddif) vdwaalsdifavg = average(vdwaalsdif) eeldifavg = average(eeldif) egbdifavg = average(egbdif) vdwofdifavg = average(vdwofdif) eelofdifavg = average(eelofdif) ecavitydifavg = average(ecavitydif) bonddifstdev = stdev(bonddif, bonddifavg) angledifstdev = stdev(angledif, angledifavg) diheddifstdev = stdev(diheddif, diheddifavg) vdwaalsdifstdev = stdev(vdwaalsdif, vdwaalsdifavg) eeldifstdev = stdev(eeldif, eeldifavg) egbdifstdev = stdev(egbdif, egbdifavg) vdwofdifstdev = stdev(vdwofdif, vdwofdifavg) eelofdifstdev = stdev(eelofdif, eelofdifavg) ecavitydifstdev = stdev(ecavitydif, ecavitydifavg) if apbs != 1: edisperdifavg = average(edisperdif) edisperdifstdev = stdev(edisperdif, edisperdifavg) gascomavg = vdwaalsavg[0] + eelavg[0] gasrecavg = vdwaalsavg[1] + eelavg[1] gasligavg = vdwaalsavg[2] + eelavg[2] gasdifavg = vdwaalsdifavg + eeldifavg gascomstd = vdwaalsstdev[0] ** 2 + eelstdev[0] ** 2 gasrecstd = vdwaalsstdev[1] ** 2 + eelstdev[1] ** 2 gasligstd = vdwaalsstdev[2] ** 2 + eelstdev[2] ** 2 gasdifstd = vdwaalsdifstdev ** 2 + eeldifstdev ** 2 solvcomavg = egbavg[0] + ecavityavg[0] solvrecavg = egbavg[1] + ecavityavg[1] solvligavg = egbavg[2] + ecavityavg[2] solvdifavg = egbdifavg + ecavitydifavg solvcomstd = math.sqrt(egbstdev[0] ** 2 + ecavitystdev[0] ** 2) solvrecstd = math.sqrt(egbstdev[1] ** 2 + ecavitystdev[1] ** 2) solvligstd = math.sqrt(egbstdev[2] ** 2 + ecavitystdev[2] ** 2) solvdifstd = math.sqrt(egbdifstdev ** 2 + ecavitydifstdev ** 2) if onetraj: failed = False failedof = False for x in range(len(bonddif)): if (abs(bonddif[x]) > 0.001 or abs(angledif[x]) > 0.001 or abs(diheddif[x]) > 0.001) and not failed: finaloutput.write(""" WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS (BOND, ANGLE, AND/OR DIHED). CHECK YOUR INPUT FILES AND SYSTEM SETUP. THESE RESULTS MAY NOT BE RELIABLE (check differences)! """) failed = True elif (abs(vdwofdif[x]) > 0.001 or abs(eelofdif[x]) > 0.001) and not failedof: finaloutput.write(""" WARNING: INCONSISTENCIES EXIST WITHIN 1-4 NON-BONDED TERMS CHECK YOUR INPUT FILES AND SYSTEM SETUP. THESE RESULTS MAY NOT BE RELIABLE (check differences)! The most common cause of this is inconsistent charge definitions across topology files. """) failedof = True if not onetraj or failed or verbose == 2: gascomavg += bondavg[0] + angleavg[0] + dihedavg[0] gasrecavg += bondavg[1] + angleavg[1] + dihedavg[1] gasligavg += bondavg[2] + angleavg[2] + dihedavg[2] gasdifavg += bonddifavg + angledifavg + diheddifavg gascomstd += bondstdev[0]**2 + anglestdev[0]**2 + dihedstdev[0]**2 gasrecstd += bondstdev[1]**2 + anglestdev[1]**2 + dihedstdev[1]**2 gasligstd += bondstdev[2]**2 + anglestdev[2]**2 + dihedstdev[2]**2 gasdifstd += bonddifstdev**2 + angledifstdev**2 + diheddifstdev**2 else: totalavg[0] = totalavg[0]-bondavg[0]-angleavg[0]-dihedavg[0] totalavg[1] = totalavg[1]-bondavg[1]-angleavg[1]-dihedavg[1] totalavg[2] = totalavg[2]-bondavg[2]-angleavg[2]-dihedavg[2] if not onetraj or failedof or verbose == 2: gascomavg += eelofavg[0] + vdwofavg[0] gasrecavg += eelofavg[1] + vdwofavg[1] gasligavg += eelofavg[2] + vdwofavg[2] gasdifavg += eelofdifavg + vdwofdifavg gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 gascomstd += eelofstdev[0]**2 + vdwofstdev[0]**2 else: totalavg[0] = totalavg[0]-eelofavg[0]-vdwofavg[0] totalavg[1] = totalavg[1]-eelofavg[1]-vdwofavg[1] totalavg[2] = totalavg[2]-eelofavg[2]-vdwofavg[2] gascomstd = math.sqrt(gascomstd) gasrecstd = math.sqrt(gasrecstd) gasligstd = math.sqrt(gasligstd) gasdifstd = math.sqrt(gasdifstd) finaloutput.write(header + '\nPOISSON BOLTZMANN:\n\n') if verbose >= 1 or failed or failedof: finaloutput.write('Complex:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: finaloutput.write('BOND ' + Format(bondavg[0],15,4) + Format(bondstdev[0],17,4) + Format(bondstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[0],15,4) + Format(anglestdev[0],17,4) + Format(anglestdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[0],15,4) + Format(dihedstdev[0],17,4) + Format(dihedstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[0],15,4) + Format(vdwaalsstdev[0],17,4) + Format(vdwaalsstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[0],15,4) + Format(eelstdev[0],17,4) + Format(eelstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: finaloutput.write('1-4 VDW ' + Format(vdwofavg[0],15,4) + Format(vdwofstdev[0],17,4) + Format(vdwofstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[0],15,4) + Format(eelofstdev[0],17,4) + Format(eelofstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EPB ' + Format(egbavg[0],15,4) + Format(egbstdev[0],17,4) + Format(egbstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') if apbs != 1: finaloutput.write('ECAVITY ' + Format(ecavityavg[0],15,4) + Format(ecavitystdev[0],17,4) + Format(ecavitystdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') if edisperavg[0] > 0.0001 and edisperavg[1] > 0.0001 and edisperavg[2] > 0.0001: finaloutput.write('EDISPER ' + Format(edisperavg[0],15,4) + Format(edisperstdev[0],17,4) + Format(edisperstdev[0]/math.sqrt(float(numframes)),15,4) + '\n\n') else: finaloutput.write('ENPOLAR ' + Format(ecavityavg[0],15,4) + Format(ecavitystdev[0],17,4) + Format(ecavitystdev[0]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gascomavg,15,4) + Format(gascomstd,17,4) + Format(gascomstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvcomavg,15,4) + Format(solvcomstd,17,4) + Format(solvcomstd/math.sqrt( numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[0],15,4) + Format(totalstdev[0],17,4) + Format(totalstdev[0]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') finaloutput.write('Receptor:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: finaloutput.write('BOND ' + Format(bondavg[1],15,4) + Format(bondstdev[1],17,4) + Format(bondstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[1],15,4) + Format(anglestdev[1],17,4) + Format(anglestdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[1],15,4) + Format(dihedstdev[1],17,4) + Format(dihedstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[1],15,4) + Format(vdwaalsstdev[1],17,4) + Format(vdwaalsstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[1],15,4) + Format(eelstdev[1],17,4) + Format(eelstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: finaloutput.write('1-4 VDW ' + Format(vdwofavg[1],15,4) + Format(vdwofstdev[1],17,4) + Format(vdwofstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[1],15,4) + Format(eelofstdev[1],17,4) + Format(eelofstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EPB ' + Format(egbavg[1],15,4) + Format(egbstdev[1],17,4) + Format(egbstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') if apbs != 1: finaloutput.write('ECAVITY ' + Format(ecavityavg[1],15,4) + Format(ecavitystdev[1],17,4) + Format(ecavitystdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') if edisperavg[0] > 0.0001 and edisperavg[1] > 0.0001 and edisperavg[2] > 0.0001: finaloutput.write('EDISPER ' + Format(edisperavg[1],15,4) + Format(edisperstdev[1],17,4) + Format(edisperstdev[1]/math.sqrt(float(numframes)),15,4) + '\n\n') else: finaloutput.write('ENPOLAR ' + Format(ecavityavg[1],15,4) + Format(ecavitystdev[1],17,4) + Format(ecavitystdev[1]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gasrecavg,15,4) + Format(gasrecstd,17,4) + Format(gasrecstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvrecavg,15,4) + Format(solvrecstd,17,4) + Format(solvrecstd/math.sqrt( numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[1],15,4) + Format(totalstdev[1],17,4) + Format(totalstdev[1]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') finaloutput.write('Ligand:\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: finaloutput.write('BOND ' + Format(bondavg[2],15,4) + Format(bondstdev[2],17,4) + Format(bondstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angleavg[2],15,4) + Format(anglestdev[2],17,4) + Format(anglestdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(dihedavg[2],15,4) + Format(dihedstdev[2],17,4) + Format(dihedstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsavg[2],15,4) + Format(vdwaalsstdev[2],17,4) + Format(vdwaalsstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eelavg[2],15,4) + Format(eelstdev[2],17,4) + Format(eelstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: finaloutput.write('1-4 VDW ' + Format(vdwofavg[2],15,4) + Format(vdwofstdev[2],17,4) + Format(vdwofstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofavg[2],15,4) + Format(eelofstdev[2],17,4) + Format(eelofstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EPB ' + Format(egbavg[2],15,4) + Format(egbstdev[2],17,4) + Format(egbstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') if apbs != 1: finaloutput.write('ECAVITY ' + Format(ecavityavg[2],15,4) + Format(ecavitystdev[2],17,4) + Format(ecavitystdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') if edisperavg[0] > 0.0001 and edisperavg[1] > 0.0001 and edisperavg[2] > 0.0001: finaloutput.write('EDISPER ' + Format(edisperavg[2],15,4) + Format(edisperstdev[2],17,4) + Format(edisperstdev[2]/math.sqrt(float(numframes)),15,4) + '\n\n') else: finaloutput.write('ENPOLAR ' + Format(ecavityavg[2],15,4) + Format(ecavitystdev[2],17,4) + Format(ecavitystdev[2]/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nG gas ' + Format(gasligavg,15,4) + Format(gasligstd,17,4) + Format(gasligstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('G solv ' + Format(solvligavg,15,4) + Format(solvligstd,17,4) + Format(solvligstd/math.sqrt( numframes),15,4) + '\n\n') finaloutput.write('TOTAL ' + Format(totalavg[2],15,4) + Format(totalstdev[2],17,4) + Format(totalstdev[2]/ math.sqrt(float(numframes)),15,4) + '\n\n\n') finaloutput.write('Differences (Complex - Receptor - Ligand):') finaloutput.write('\nEnergy Component Average Std. Dev. Std. Err. of Mean\n') finaloutput.write('-------------------------------------------------------------------------------\n') if not onetraj or failed or verbose == 2: finaloutput.write('BOND ' + Format(bonddifavg,15,4) + Format(bonddifstdev,17,4) + Format(bonddifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('ANGLE ' + Format(angledifavg,15,4) + Format(angledifstdev,17,4) + Format(angledifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('DIHED ' + Format(diheddifavg,15,4) + Format(diheddifstdev,17,4) + Format(diheddifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('VDWAALS ' + Format(vdwaalsdifavg,15,4) + Format(vdwaalsdifstdev,17,4) + Format(vdwaalsdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EEL ' + Format(eeldifavg,15,4) + Format(eeldifstdev,17,4) + Format(eeldifstdev/ math.sqrt(float(numframes)),15,4) + '\n') if not onetraj or failedof or verbose == 2: finaloutput.write('1-4 VDW ' + Format(vdwofdifavg,15,4) + Format(vdwofdifstdev,17,4) + Format(vdwofdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('1-4 EEL ' + Format(eelofdifavg,15,4) + Format(eelofdifstdev,17,4) + Format(eelofdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('EPB ' + Format(egbdifavg,15,4) + Format(egbdifstdev,17,4) + Format(egbdifstdev/ math.sqrt(float(numframes)),15,4) + '\n') if apbs != 1: finaloutput.write('ECAVITY ' + Format(ecavitydifavg,15,4) + Format(ecavitydifstdev,17,4) + Format(ecavitydifstdev/ math.sqrt(float(numframes)),15,4) + '\n') if edisperdifavg > 0.0001: finaloutput.write('EDISPER ' + Format(edisperdifavg,15,4) + Format(edisperdifstdev,17,4) + Format(edisperdifstdev/math.sqrt(float(numframes)),15,4) + '\n\n') else: finaloutput.write('ENPOLAR ' + Format(ecavitydifavg,15,4) + Format(ecavitydifstdev,17,4) + Format(ecavitydifstdev/ math.sqrt(float(numframes)),15,4) + '\n') finaloutput.write('\nDELTA G gas ' + Format(gasdifavg,15,4) + Format(gasdifstd,17,4) + Format(gasdifstd/ math.sqrt(numframes),15,4) + '\n') finaloutput.write('DELTA G solv ' + Format(solvdifavg,15,4) + Format(solvdifstd,17,4) + Format(solvdifstd/ math.sqrt(numframes),15,4) + '\n\n') finaloutput.write('\n\n DELTA G binding = ' + Format(deltagavg,10,4) + ' +/-' + Format(deltagstdev,7,4) + Format(deltagstdev/math.sqrt(float(numframes)),18,4) + '\n') return [deltagavg,deltagstdev] #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def decoutput12(output, debug, verbose, frames, comprm, recprm, ligprm, ligstart, type, decomp, surften): import sys, math frames = float(frames) mdout = open("_MMPBSA_complex_" + type + ".mdout",'r') first = True # Logical variable -- is this the first frame we've analyzed? intot = False # if we are currently analyzing the total energy decomp inbac = False # if we are currently analyzing the backbone energy decomp insc = False # if we are currently analyzing the sidechain energy decomp skip = False # if we've just seen "PRINT" and must now skip lines until we reach data onnum = 0 # ON residue NUMber -- keeps track of which residue we're extracting data from cress = [] # Complex RESidueS -- array that will contain all printed decomp residues in complex # All arrays below are arrays of 2-element arrays. First element is a running sum for the average, second # element is a running sum of squares for the standard deviation ctint = [] # Complex Total INTernal decomposition ctvdw = [] # Complex Total VDW decomposition cteel = [] # Complex Total EEL (electrostatic) decomposition ctpol = [] # Complex Total POLar solvation decomposition ctsas = [] # Complex Total non-polar solvation decomposition csint = [] # Complex Sidechain INTernal decomposition csvdw = [] # Complex Sidechain VDW decomposition cseel = [] # Complex Sidechain EEL (elecrostatic) decomposition cspol = [] # Complex Sidechain POLar solvation decomposition cssas = [] # Complex Sidechain non-polar solvation decomposition cbint = [] # Complex Backbone INTernal decomposition cbvdw = [] # Complex Backbone VDW decomposition cbeel = [] # Complex Backbone EEL decomposition cbpol = [] # Complex Backbone POLar solvation decomposition cbsas = [] # Complex Backbone non-polar solvation decomposition for line in mdout: if not intot and not inbac and not insc: # If we are not currently analyzing any decomp output if not "PRINT DECOMP" in line: continue elif "TOTAL" in line: # we've found the TOTAL decomp data intot = True # we are now in total skip = True # we must skip the filler-lines between PRINT DECOMP and actual data elif "SIDECHAIN" in line: insc = True # same as for TOTAL skip = True elif "BACKBONE" in line: inbac = True # same as for TOTAL skip = True if intot: if skip: # if we are still skipping lines... (to distinguish from blank lines at end of data) if not "TDC" in line: continue # skip ahead to next line if we haven't started analyzing data else: skip = False # we've found the data, no more skipping. if not "TDC" in line: # if we've reached a blank line now, we've finished the data onnum = 0 # we're now on residue 0 again, as we're going to a new data block intot = False # we're no longer in total block continue # we've done what we need to -- skip to next line words = line.split() # now we're in to analyzing data -- split line at whitespace if first: # if we're on the first frame, we must build the arrays cress.append(int(words[1])) ctint.append([float(words[2]),float(words[2])**2]) ctvdw.append([float(words[3]),float(words[3])**2]) cteel.append([float(words[4]),float(words[4])**2]) ctpol.append([float(words[5]),float(words[5])**2]) ctsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: # otherwise just append to the array for the given residue ctint[onnum][0] += float(words[2]) ctvdw[onnum][0] += float(words[3]) cteel[onnum][0] += float(words[4]) ctpol[onnum][0] += float(words[5]) ctsas[onnum][0] += float(words[6]) * surften ctint[onnum][1] += float(words[2]) ** 2 ctvdw[onnum][1] += float(words[3]) ** 2 cteel[onnum][1] += float(words[4]) ** 2 ctpol[onnum][1] += float(words[5]) ** 2 ctsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if insc: # same logic as we applied to TOTAL if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: csint.append([float(words[2]),float(words[2])**2]) csvdw.append([float(words[3]),float(words[3])**2]) cseel.append([float(words[4]),float(words[4])**2]) cspol.append([float(words[5]),float(words[5])**2]) cssas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: csint[onnum][0] += float(words[2]) csvdw[onnum][0] += float(words[3]) cseel[onnum][0] += float(words[4]) cspol[onnum][0] += float(words[5]) cssas[onnum][0] += float(words[6]) * surften csint[onnum][1] += float(words[2]) ** 2 csvdw[onnum][1] += float(words[3]) ** 2 cseel[onnum][1] += float(words[4]) ** 2 cspol[onnum][1] += float(words[5]) ** 2 cssas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if inbac: # same logic as we applied to TOTAL with one change... if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False # backbone is the last block to be printed, we now cannot be in first frame anymore continue words = line.split() if first: cbint.append([float(words[2]),float(words[2])**2]) cbvdw.append([float(words[3]),float(words[3])**2]) cbeel.append([float(words[4]),float(words[4])**2]) cbpol.append([float(words[5]),float(words[5])**2]) cbsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: cbint[onnum][0] += float(words[2]) cbvdw[onnum][0] += float(words[3]) cbeel[onnum][0] += float(words[4]) cbpol[onnum][0] += float(words[5]) cbsas[onnum][0] += float(words[6]) * surften cbint[onnum][1] += float(words[2]) ** 2 cbvdw[onnum][1] += float(words[3]) ** 2 cbeel[onnum][1] += float(words[4]) ** 2 cbpol[onnum][1] += float(words[5]) ** 2 cbsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(cress)): # we now compute averages and standard deviations ctint[x][0] /= frames ctint[x][1] = math.sqrt(abs(ctint[x][1]/frames - ctint[x][0] ** 2)) ctvdw[x][0] /= frames ctvdw[x][1] = math.sqrt(abs(ctvdw[x][1]/frames - ctvdw[x][0] ** 2)) cteel[x][0] /= frames cteel[x][1] = math.sqrt(abs(cteel[x][1]/frames - cteel[x][0] ** 2)) ctpol[x][0] /= frames ctpol[x][1] = math.sqrt(abs(ctpol[x][1]/frames - ctpol[x][0] ** 2)) ctsas[x][0] /= frames ctsas[x][1] = math.sqrt(abs(ctsas[x][1]/frames - ctsas[x][0] ** 2)) csint[x][0] /= frames csint[x][1] = math.sqrt(abs(csint[x][1]/frames - csint[x][0] ** 2)) csvdw[x][0] /= frames csvdw[x][1] = math.sqrt(abs(csvdw[x][1]/frames - csvdw[x][0] ** 2)) cseel[x][0] /= frames cseel[x][1] = math.sqrt(abs(cseel[x][1]/frames - cseel[x][0] ** 2)) cspol[x][0] /= frames cspol[x][1] = math.sqrt(abs(cspol[x][1]/frames - cspol[x][0] ** 2)) cssas[x][0] /= frames cssas[x][1] = math.sqrt(abs(cssas[x][1]/frames - cssas[x][0] ** 2)) cbint[x][0] /= frames cbint[x][1] = math.sqrt(abs(cbint[x][1]/frames - cbint[x][0] ** 2)) cbvdw[x][0] /= frames cbvdw[x][1] = math.sqrt(abs(cbvdw[x][1]/frames - cbvdw[x][0] ** 2)) cbeel[x][0] /= frames cbeel[x][1] = math.sqrt(abs(cbeel[x][1]/frames - cbeel[x][0] ** 2)) cbpol[x][0] /= frames cbpol[x][1] = math.sqrt(abs(cbpol[x][1]/frames - cbpol[x][0] ** 2)) cbsas[x][0] /= frames cbsas[x][1] = math.sqrt(abs(cbsas[x][1]/frames - cbsas[x][0] ** 2)) # Done with complex. Do the same thing for the receptor, only name arrays starting with 'r' # instead of 'c' for Receptor. mdout = open("_MMPBSA_receptor_" + type + ".mdout",'r') first = True intot = False inbac = False insc = False skip = False onnum = 0 rress = [] rtint = [] rtvdw = [] rteel = [] rtpol = [] rtsas = [] rsint = [] rsvdw = [] rseel = [] rspol = [] rssas = [] rbint = [] rbvdw = [] rbeel = [] rbpol = [] rbsas = [] for line in mdout: if not intot and not inbac and not insc: if not "PRINT DECOMP" in line: continue elif "TOTAL" in line: intot = True skip = True elif "SIDECHAIN" in line: insc = True skip = True elif "BACKBONE" in line: inbac = True skip = True if intot: if skip: if not "TDC" in line: continue else: skip = False if not "TDC" in line: onnum = 0 intot = False continue words = line.split() if first: rress.append(int(words[1])) rtint.append([float(words[2]),float(words[2])**2]) rtvdw.append([float(words[3]),float(words[3])**2]) rteel.append([float(words[4]),float(words[4])**2]) rtpol.append([float(words[5]),float(words[5])**2]) rtsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: rtint[onnum][0] += float(words[2]) rtvdw[onnum][0] += float(words[3]) rteel[onnum][0] += float(words[4]) rtpol[onnum][0] += float(words[5]) rtsas[onnum][0] += float(words[6]) * surften rtint[onnum][1] += float(words[2]) ** 2 rtvdw[onnum][1] += float(words[3]) ** 2 rteel[onnum][1] += float(words[4]) ** 2 rtpol[onnum][1] += float(words[5]) ** 2 rtsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if insc: if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: rsint.append([float(words[2]),float(words[2])**2]) rsvdw.append([float(words[3]),float(words[3])**2]) rseel.append([float(words[4]),float(words[4])**2]) rspol.append([float(words[5]),float(words[5])**2]) rssas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: rsint[onnum][0] += float(words[2]) rsvdw[onnum][0] += float(words[3]) rseel[onnum][0] += float(words[4]) rspol[onnum][0] += float(words[5]) rssas[onnum][0] += float(words[6]) * surften rsint[onnum][1] += float(words[2]) ** 2 rsvdw[onnum][1] += float(words[3]) ** 2 rseel[onnum][1] += float(words[4]) ** 2 rspol[onnum][1] += float(words[5]) ** 2 rssas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if inbac: if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False continue words = line.split() if first: rbint.append([float(words[2]),float(words[2])**2]) rbvdw.append([float(words[3]),float(words[3])**2]) rbeel.append([float(words[4]),float(words[4])**2]) rbpol.append([float(words[5]),float(words[5])**2]) rbsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: rbint[onnum][0] += float(words[2]) rbvdw[onnum][0] += float(words[3]) rbeel[onnum][0] += float(words[4]) rbpol[onnum][0] += float(words[5]) rbsas[onnum][0] += float(words[6]) * surften rbint[onnum][1] += float(words[2]) ** 2 rbvdw[onnum][1] += float(words[3]) ** 2 rbeel[onnum][1] += float(words[4]) ** 2 rbpol[onnum][1] += float(words[5]) ** 2 rbsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(rress)): rtint[x][0] /= frames rtint[x][1] = math.sqrt(abs(rtint[x][1]/frames - rtint[x][0] ** 2)) rtvdw[x][0] /= frames rtvdw[x][1] = math.sqrt(abs(rtvdw[x][1]/frames - rtvdw[x][0] ** 2)) rteel[x][0] /= frames rteel[x][1] = math.sqrt(abs(rteel[x][1]/frames - rteel[x][0] ** 2)) rtpol[x][0] /= frames rtpol[x][1] = math.sqrt(abs(rtpol[x][1]/frames - rtpol[x][0] ** 2)) rtsas[x][0] /= frames rtsas[x][1] = math.sqrt(abs(rtsas[x][1]/frames - rtsas[x][0] ** 2)) rsint[x][0] /= frames rsint[x][1] = math.sqrt(abs(rsint[x][1]/frames - rsint[x][0] ** 2)) rsvdw[x][0] /= frames rsvdw[x][1] = math.sqrt(abs(rsvdw[x][1]/frames - rsvdw[x][0] ** 2)) rseel[x][0] /= frames rseel[x][1] = math.sqrt(abs(rseel[x][1]/frames - rseel[x][0] ** 2)) rspol[x][0] /= frames rspol[x][1] = math.sqrt(abs(rspol[x][1]/frames - rspol[x][0] ** 2)) rssas[x][0] /= frames rssas[x][1] = math.sqrt(abs(rssas[x][1]/frames - rssas[x][0] ** 2)) rbint[x][0] /= frames rbint[x][1] = math.sqrt(abs(rbint[x][1]/frames - rbint[x][0] ** 2)) rbvdw[x][0] /= frames rbvdw[x][1] = math.sqrt(abs(rbvdw[x][1]/frames - rbvdw[x][0] ** 2)) rbeel[x][0] /= frames rbeel[x][1] = math.sqrt(abs(rbeel[x][1]/frames - rbeel[x][0] ** 2)) rbpol[x][0] /= frames rbpol[x][1] = math.sqrt(abs(rbpol[x][1]/frames - rbpol[x][0] ** 2)) rbsas[x][0] /= frames rbsas[x][1] = math.sqrt(abs(rbsas[x][1]/frames - rbsas[x][0] ** 2)) # Same as receptor above mdout = open("_MMPBSA_ligand_" + type + ".mdout",'r') first = True intot = False inbac = False insc = False skip = False onnum = 0 lress = [] ltint = [] ltvdw = [] lteel = [] ltpol = [] ltsas = [] lsint = [] lsvdw = [] lseel = [] lspol = [] lssas = [] lbint = [] lbvdw = [] lbeel = [] lbpol = [] lbsas = [] for line in mdout: if not intot and not inbac and not insc: if not "PRINT DECOMP" in line: continue elif "TOTAL" in line: intot = True skip = True elif "SIDECHAIN" in line: insc = True skip = True elif "BACKBONE" in line: inbac = True skip = True if intot: if skip: if not "TDC" in line: continue else: skip = False if not "TDC" in line: onnum = 0 intot = False continue words = line.split() if first: lress.append(int(words[1])) ltint.append([float(words[2]),float(words[2])**2]) ltvdw.append([float(words[3]),float(words[3])**2]) lteel.append([float(words[4]),float(words[4])**2]) ltpol.append([float(words[5]),float(words[5])**2]) ltsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: ltint[onnum][0] += float(words[2]) ltvdw[onnum][0] += float(words[3]) lteel[onnum][0] += float(words[4]) ltpol[onnum][0] += float(words[5]) ltsas[onnum][0] += float(words[6]) * surften ltint[onnum][1] += float(words[2]) ** 2 ltvdw[onnum][1] += float(words[3]) ** 2 lteel[onnum][1] += float(words[4]) ** 2 ltpol[onnum][1] += float(words[5]) ** 2 ltsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if insc: if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: lsint.append([float(words[2]),float(words[2])**2]) lsvdw.append([float(words[3]),float(words[3])**2]) lseel.append([float(words[4]),float(words[4])**2]) lspol.append([float(words[5]),float(words[5])**2]) lssas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: lsint[onnum][0] += float(words[2]) lsvdw[onnum][0] += float(words[3]) lseel[onnum][0] += float(words[4]) lspol[onnum][0] += float(words[5]) lssas[onnum][0] += float(words[6]) * surften lsint[onnum][1] += float(words[2]) ** 2 lsvdw[onnum][1] += float(words[3]) ** 2 lseel[onnum][1] += float(words[4]) ** 2 lspol[onnum][1] += float(words[5]) ** 2 lssas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 if inbac: if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False continue words = line.split() if first: lbint.append([float(words[2]),float(words[2])**2]) lbvdw.append([float(words[3]),float(words[3])**2]) lbeel.append([float(words[4]),float(words[4])**2]) lbpol.append([float(words[5]),float(words[5])**2]) lbsas.append([float(words[6])*surften,(float(words[6])*surften)**2]) else: lbint[onnum][0] += float(words[2]) lbvdw[onnum][0] += float(words[3]) lbeel[onnum][0] += float(words[4]) lbpol[onnum][0] += float(words[5]) lbsas[onnum][0] += float(words[6]) * surften lbint[onnum][1] += float(words[2]) ** 2 lbvdw[onnum][1] += float(words[3]) ** 2 lbeel[onnum][1] += float(words[4]) ** 2 lbpol[onnum][1] += float(words[5]) ** 2 lbsas[onnum][1] += float(words[6]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(lress)): ltint[x][0] /= frames ltint[x][1] = math.sqrt(abs(ltint[x][1]/frames - ltint[x][0] ** 2)) ltvdw[x][0] /= frames ltvdw[x][1] = math.sqrt(abs(ltvdw[x][1]/frames - ltvdw[x][0] ** 2)) lteel[x][0] /= frames lteel[x][1] = math.sqrt(abs(lteel[x][1]/frames - lteel[x][0] ** 2)) ltpol[x][0] /= frames ltpol[x][1] = math.sqrt(abs(ltpol[x][1]/frames - ltpol[x][0] ** 2)) ltsas[x][0] /= frames ltsas[x][1] = math.sqrt(abs(ltsas[x][1]/frames - ltsas[x][0] ** 2)) lsint[x][0] /= frames lsint[x][1] = math.sqrt(abs(lsint[x][1]/frames - lsint[x][0] ** 2)) lsvdw[x][0] /= frames lsvdw[x][1] = math.sqrt(abs(lsvdw[x][1]/frames - lsvdw[x][0] ** 2)) lseel[x][0] /= frames lseel[x][1] = math.sqrt(abs(lseel[x][1]/frames - lseel[x][0] ** 2)) lspol[x][0] /= frames lspol[x][1] = math.sqrt(abs(lspol[x][1]/frames - lspol[x][0] ** 2)) lssas[x][0] /= frames lssas[x][1] = math.sqrt(abs(lssas[x][1]/frames - lssas[x][0] ** 2)) lbint[x][0] /= frames lbint[x][1] = math.sqrt(abs(lbint[x][1]/frames - lbint[x][0] ** 2)) lbvdw[x][0] /= frames lbvdw[x][1] = math.sqrt(abs(lbvdw[x][1]/frames - lbvdw[x][0] ** 2)) lbeel[x][0] /= frames lbeel[x][1] = math.sqrt(abs(lbeel[x][1]/frames - lbeel[x][0] ** 2)) lbpol[x][0] /= frames lbpol[x][1] = math.sqrt(abs(lbpol[x][1]/frames - lbpol[x][0] ** 2)) lbsas[x][0] /= frames lbsas[x][1] = math.sqrt(abs(lbsas[x][1]/frames - lbsas[x][0] ** 2)) # get the list of all residue names in the complex, receptor, and ligand so we can label the printed residues in the output com_resnames = getallresinfo(comprm,"RESIDUE_LABEL") rec_resnames = getallresinfo(recprm,"RESIDUE_LABEL") lig_resnames = getallresinfo(ligprm,"RESIDUE_LABEL") # print output if decomp == 1: output.write("idecomp = 1: Decomposition per-residue adding 1-4 interactions added to Internal.\n") if decomp == 2: output.write("idecomp = 2: Decomposition per-residue adding 1-4 interactions added to EEL and VDW.\n") if ligstart == -1 and verbose <= 2: output.write("DELTAs can only be printed with the default mask assignments\n\n") verbose += 2 if type == 'gb': output.write("Energy Decomposition Analysis (All units kcal/mol): Generalized Born solvent\n\n") elif type == 'pb': output.write("Energy Decomposition Analysis (All units kcal/mol): Poisson Boltzmann solvent\n\n") if verbose == 2 or verbose == 3: output.write("Complex:" + """ Total Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = ctint[x][0] + ctvdw[x][0] + cteel[x][0] + ctpol[x][0] + ctsas[x][0] total_stdev = math.sqrt(ctint[x][1]**2 + ctvdw[x][1]**2 + cteel[x][1]**2 + ctpol[x][1]**2 + ctsas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + Format(ctint[x][0],5,3) + ' +/-' + Format(ctint[x][1],3,3) + ' |' + Format(ctvdw[x][0],5,3) + ' +/-' + Format(ctvdw[x][1],3,3) + ' |' + Format(cteel[x][0],5,3) + ' +/-' + Format(cteel[x][1],3,3) + ' |' + Format(ctpol[x][0],5,3) + ' +/-' + Format(ctpol[x][1],3,3) + ' |' + Format(ctsas[x][0],5,3) + ' +/-' + Format(ctsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3: output.write(""" Sidechain Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = csint[x][0] + csvdw[x][0] + cseel[x][0] + cspol[x][0] + cssas[x][0] total_stdev = math.sqrt(csint[x][1]**2 + csvdw[x][1]**2 + cseel[x][1]**2 + cspol[x][1]**2 + cssas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + Format(csint[x][0],5,3) + ' +/-' + Format(csint[x][1],3,3) + ' |' + Format(csvdw[x][0],5,3) + ' +/-' + Format(csvdw[x][1],3,3) + ' |' + Format(cseel[x][0],5,3) + ' +/-' + Format(cseel[x][1],3,3) + ' |' + Format(cspol[x][0],5,3) + ' +/-' + Format(cspol[x][1],3,3) + ' |' + Format(cssas[x][0],5,3) + ' +/-' + Format(cssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') output.write(""" Backbone Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = cbint[x][0] + cbvdw[x][0] + cbeel[x][0] + cbpol[x][0] + cbsas[x][0] total_stdev = math.sqrt(cbint[x][1]**2 + cbvdw[x][1]**2 + cbeel[x][1]**2 + cbpol[x][1]**2 + cbsas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + Format(cbint[x][0],5,3) + ' +/-' + Format(cbint[x][1],3,3) + ' |' + Format(cbvdw[x][0],5,3) + ' +/-' + Format(cbvdw[x][1],3,3) + ' |' + Format(cbeel[x][0],5,3) + ' +/-' + Format(cbeel[x][1],3,3) + ' |' + Format(cbpol[x][0],5,3) + ' +/-' + Format(cbpol[x][1],3,3) + ' |' + Format(cbsas[x][0],5,3) + ' +/-' + Format(cbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 2 or verbose == 3: output.write("\nReceptor:" + """ Total Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(rress)): total_avg = rtint[x][0] + rtvdw[x][0] + rteel[x][0] + rtpol[x][0] + rtsas[x][0] total_stdev = math.sqrt(rtint[x][1]**2 + rtvdw[x][1]**2 + rteel[x][1]**2 + rtpol[x][1]**2 + rtsas[x][1]**2) output.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' |' + Format(rtint[x][0],5,3) + ' +/-' + Format(rtint[x][1],3,3) + ' |' + Format(rtvdw[x][0],5,3) + ' +/-' + Format(rtvdw[x][1],3,3) + ' |' + Format(rteel[x][0],5,3) + ' +/-' + Format(rteel[x][1],3,3) + ' |' + Format(rtpol[x][0],5,3) + ' +/-' + Format(rtpol[x][1],3,3) + ' |' + Format(rtsas[x][0],5,3) + ' +/-' + Format(rtsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3: output.write(""" Sidechain Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(rress)): total_avg = rsint[x][0] + rsvdw[x][0] + rseel[x][0] + rspol[x][0] + rssas[x][0] total_stdev = math.sqrt(rsint[x][1]**2 + rsvdw[x][1]**2 + rseel[x][1]**2 + rspol[x][1]**2 + rssas[x][1]**2) output.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' |' + Format(rsint[x][0],5,3) + ' +/-' + Format(rsint[x][1],3,3) + ' |' + Format(rsvdw[x][0],5,3) + ' +/-' + Format(rsvdw[x][1],3,3) + ' |' + Format(rseel[x][0],5,3) + ' +/-' + Format(rseel[x][1],3,3) + ' |' + Format(rspol[x][0],5,3) + ' +/-' + Format(rspol[x][1],3,3) + ' |' + Format(rssas[x][0],5,3) + ' +/-' + Format(rssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') output.write(""" Backbone Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(rress)): total_avg = rbint[x][0] + rbvdw[x][0] + rbeel[x][0] + rbpol[x][0] + rbsas[x][0] total_stdev = math.sqrt(rbint[x][1]**2 + rbvdw[x][1]**2 + rbeel[x][1]**2 + rbpol[x][1]**2 + rbsas[x][1]**2) output.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' |' + Format(rbint[x][0],5,3) + ' +/-' + Format(rbint[x][1],3,3) + ' |' + Format(rbvdw[x][0],5,3) + ' +/-' + Format(rbvdw[x][1],3,3) + ' |' + Format(rbeel[x][0],5,3) + ' +/-' + Format(rbeel[x][1],3,3) + ' |' + Format(rbpol[x][0],5,3) + ' +/-' + Format(rbpol[x][1],3,3) + ' |' + Format(rbsas[x][0],5,3) + ' +/-' + Format(rbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 2 or verbose == 3: output.write("\nLigand:" + """ Total Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(lress)): total_avg = ltint[x][0] + ltvdw[x][0] + lteel[x][0] + ltpol[x][0] + ltsas[x][0] total_stdev = math.sqrt(ltint[x][1]**2 + ltvdw[x][1]**2 + lteel[x][1]**2 + ltpol[x][1]**2 + ltsas[x][1]**2) output.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' |' + Format(ltint[x][0],5,3) + ' +/-' + Format(ltint[x][1],3,3) + ' |' + Format(ltvdw[x][0],5,3) + ' +/-' + Format(ltvdw[x][1],3,3) + ' |' + Format(lteel[x][0],5,3) + ' +/-' + Format(lteel[x][1],3,3) + ' |' + Format(ltpol[x][0],5,3) + ' +/-' + Format(ltpol[x][1],3,3) + ' |' + Format(ltsas[x][0],5,3) + ' +/-' + Format(ltsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3: output.write(""" Sidechain Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(lress)): total_avg = lsint[x][0] + lsvdw[x][0] + lseel[x][0] + lspol[x][0] + lssas[x][0] total_stdev = math.sqrt(lsint[x][1]**2 + lsvdw[x][1]**2 + lseel[x][1]**2 + lspol[x][1]**2 + lssas[x][1]**2) output.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' |' + Format(lsint[x][0],5,3) + ' +/-' + Format(lsint[x][1],3,3) + ' |' + Format(lsvdw[x][0],5,3) + ' +/-' + Format(lsvdw[x][1],3,3) + ' |' + Format(lseel[x][0],5,3) + ' +/-' + Format(lseel[x][1],3,3) + ' |' + Format(lspol[x][0],5,3) + ' +/-' + Format(lspol[x][1],3,3) + ' |' + Format(lssas[x][0],5,3) + ' +/-' + Format(lssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') output.write(""" Backbone Energy Decomposition: Residue | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(lress)): total_avg = lbint[x][0] + lbvdw[x][0] + lbeel[x][0] + lbpol[x][0] + lbsas[x][0] total_stdev = math.sqrt(lbint[x][1]**2 + lbvdw[x][1]**2 + lbeel[x][1]**2 + lbpol[x][1]**2 + lbsas[x][1]**2) output.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' |' + Format(lbint[x][0],5,3) + ' +/-' + Format(lbint[x][1],3,3) + ' |' + Format(lbvdw[x][0],5,3) + ' +/-' + Format(lbvdw[x][1],3,3) + ' |' + Format(lbeel[x][0],5,3) + ' +/-' + Format(lbeel[x][1],3,3) + ' |' + Format(lbpol[x][0],5,3) + ' +/-' + Format(lbpol[x][1],3,3) + ' |' + Format(lbsas[x][0],5,3) + ' +/-' + Format(lbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # Time for the fun part: DELTAs. if ligstart != -1: dtint = [] # Delta Total INTernal decomp dtvdw = [] # Delta Total VDW decomp dteel = [] # Delta Total EEL decomp dtpol = [] # Delta Total POLar solvation decomp dtsas = [] # Delta Total non-polar solvation decomp dsint = [] # Delta Sidechain INTernal decomp dsvdw = [] # Delta Sidechain VDW decomp dseel = [] # Delta Sidechain EEL decomp dspol = [] # Delta Sidechain POLar solvation decomp dssas = [] # Delta Sidechain non-polar solvation decomp dbint = [] # Delta Backbone INTernal decomp dbvdw = [] # Delta Backbone VDW decomp dbeel = [] # Delta Backbone EEL decomp dbpol = [] # Delta Backbone POLar solvation decomp dbsas = [] # Delta Backbone non-polar solvation decomp loc = [] # Array that holds strings for locations of residues in ligand/receptor recdone = 0 # how many receptor residues have been used in differences -- used as an index ligdone = 0 # same as above, but for ligand for x in range(len(cress)): # accumulate the deltas # Figure out location of the complex residue if cress[x] < ligstart or cress[x] >= ligstart + len(lig_resnames): loc.append(' R ' + rec_resnames[rress[recdone]-1] + ' ' + str(rress[recdone]).rjust(3) + ' |') dtint.append([ctint[x][0] - rtint[recdone][0], math.sqrt(ctint[x][1]**2+rtint[recdone][1]**2)]) dtvdw.append([ctvdw[x][0] - rtvdw[recdone][0], math.sqrt(ctvdw[x][1]**2+rtvdw[recdone][1]**2)]) dteel.append([cteel[x][0] - rteel[recdone][0], math.sqrt(cteel[x][1]**2+rteel[recdone][1]**2)]) dtpol.append([ctpol[x][0] - rtpol[recdone][0], math.sqrt(ctpol[x][1]**2+rtpol[recdone][1]**2)]) dtsas.append([ctsas[x][0] - rtsas[recdone][0], math.sqrt(ctsas[x][1]**2+rtsas[recdone][1]**2)]) dsint.append([csint[x][0] - rsint[recdone][0], math.sqrt(csint[x][1]**2+rsint[recdone][1]**2)]) dsvdw.append([csvdw[x][0] - rsvdw[recdone][0], math.sqrt(csvdw[x][1]**2+rsvdw[recdone][1]**2)]) dseel.append([cseel[x][0] - rseel[recdone][0], math.sqrt(cseel[x][1]**2+rseel[recdone][1]**2)]) dspol.append([cspol[x][0] - rspol[recdone][0], math.sqrt(cspol[x][1]**2+rspol[recdone][1]**2)]) dssas.append([cssas[x][0] - rssas[recdone][0], math.sqrt(cssas[x][1]**2+rssas[recdone][1]**2)]) dbint.append([cbint[x][0] - rbint[recdone][0], math.sqrt(cbint[x][1]**2+rbint[recdone][1]**2)]) dbvdw.append([cbvdw[x][0] - rbvdw[recdone][0], math.sqrt(cbvdw[x][1]**2+rbvdw[recdone][1]**2)]) dbeel.append([cbeel[x][0] - rbeel[recdone][0], math.sqrt(cbeel[x][1]**2+rbeel[recdone][1]**2)]) dbpol.append([cbpol[x][0] - rbpol[recdone][0], math.sqrt(cbpol[x][1]**2+rbpol[recdone][1]**2)]) dbsas.append([cbsas[x][0] - rbsas[recdone][0], math.sqrt(cbsas[x][1]**2+rbsas[recdone][1]**2)]) recdone += 1 # elif cress[x] >= ligstart and cress[x] < ligstart + len(lig_resnames): # if it's in ligand else: loc.append(' L ' + lig_resnames[lress[ligdone]-1] + ' ' + str(lress[ligdone]).rjust(3) + ' |') dtint.append([ctint[x][0] - ltint[ligdone][0], math.sqrt(ctint[x][1]**2+ltint[ligdone][1]**2)]) dtvdw.append([ctvdw[x][0] - ltvdw[ligdone][0], math.sqrt(ctvdw[x][1]**2+ltvdw[ligdone][1]**2)]) dteel.append([cteel[x][0] - lteel[ligdone][0], math.sqrt(cteel[x][1]**2+lteel[ligdone][1]**2)]) dtpol.append([ctpol[x][0] - ltpol[ligdone][0], math.sqrt(ctpol[x][1]**2+ltpol[ligdone][1]**2)]) dtsas.append([ctsas[x][0] - ltsas[ligdone][0], math.sqrt(ctsas[x][1]**2+ltsas[ligdone][1]**2)]) dsint.append([csint[x][0] - lsint[ligdone][0], math.sqrt(csint[x][1]**2+lsint[ligdone][1]**2)]) dsvdw.append([csvdw[x][0] - lsvdw[ligdone][0], math.sqrt(csvdw[x][1]**2+lsvdw[ligdone][1]**2)]) dseel.append([cseel[x][0] - lseel[ligdone][0], math.sqrt(cseel[x][1]**2+lseel[ligdone][1]**2)]) dspol.append([cspol[x][0] - lspol[ligdone][0], math.sqrt(cspol[x][1]**2+lspol[ligdone][1]**2)]) dssas.append([cssas[x][0] - lssas[ligdone][0], math.sqrt(cssas[x][1]**2+lssas[ligdone][1]**2)]) dbint.append([cbint[x][0] - lbint[ligdone][0], math.sqrt(cbint[x][1]**2+lbint[ligdone][1]**2)]) dbvdw.append([cbvdw[x][0] - lbvdw[ligdone][0], math.sqrt(cbvdw[x][1]**2+lbvdw[ligdone][1]**2)]) dbeel.append([cbeel[x][0] - lbeel[ligdone][0], math.sqrt(cbeel[x][1]**2+lbeel[ligdone][1]**2)]) dbpol.append([cbpol[x][0] - lbpol[ligdone][0], math.sqrt(cbpol[x][1]**2+lbpol[ligdone][1]**2)]) dbsas.append([cbsas[x][0] - lbsas[ligdone][0], math.sqrt(cbsas[x][1]**2+lbsas[ligdone][1]**2)]) ligdone += 1 if verbose == 0 or verbose == 1 or verbose == 2 or verbose == 3: output.write(""" DELTAS: Total Energy Decomposition: Residue | Location | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = dtint[x][0] + dtvdw[x][0] + dteel[x][0] + dtpol[x][0] + dtsas[x][0] total_stdev = math.sqrt(dtint[x][1]**2 + dtvdw[x][1]**2 + dteel[x][1]**2 + dtpol[x][1]**2 + dtsas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + loc[x] + Format(dtint[x][0],5,3) + ' +/-' + Format(dtint[x][1],3,3) + ' |' + Format(dtvdw[x][0],5,3) + ' +/-' + Format(dtvdw[x][1],3,3) + ' |' + Format(dteel[x][0],5,3) + ' +/-' + Format(dteel[x][1],3,3) + ' |' + Format(dtpol[x][0],5,3) + ' +/-' + Format(dtpol[x][1],3,3) + ' |' + Format(dtsas[x][0],5,3) + ' +/-' + Format(dtsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 1 or verbose == 3: output.write(""" Sidechain Energy Decomposition: Residue | Location | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = dsint[x][0] + dsvdw[x][0] + dseel[x][0] + dspol[x][0] + dssas[x][0] total_stdev = math.sqrt(dsint[x][1]**2 + dsvdw[x][1]**2 + dseel[x][1]**2 + dspol[x][1]**2 + dssas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + loc[x] + Format(dsint[x][0],5,3) + ' +/-' + Format(dsint[x][1],3,3) + ' |' + Format(dsvdw[x][0],5,3) + ' +/-' + Format(dsvdw[x][1],3,3) + ' |' + Format(dseel[x][0],5,3) + ' +/-' + Format(dseel[x][1],3,3) + ' |' + Format(dspol[x][0],5,3) + ' +/-' + Format(dspol[x][1],3,3) + ' |' + Format(dssas[x][0],5,3) + ' +/-' + Format(dssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') output.write(""" Backbone Energy Decomposition: Residue | Location | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ------------------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): total_avg = dbint[x][0] + dbvdw[x][0] + dbeel[x][0] + dbpol[x][0] + dbsas[x][0] total_stdev = math.sqrt(dbint[x][1]**2 + dbvdw[x][1]**2 + dbeel[x][1]**2 + dbpol[x][1]**2 + dbsas[x][1]**2) output.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' |' + loc[x] + Format(dbint[x][0],5,3) + ' +/-' + Format(dbint[x][1],3,3) + ' |' + Format(dbvdw[x][0],5,3) + ' +/-' + Format(dbvdw[x][1],3,3) + ' |' + Format(dbeel[x][0],5,3) + ' +/-' + Format(dbeel[x][1],3,3) + ' |' + Format(dbpol[x][0],5,3) + ' +/-' + Format(dbpol[x][1],3,3) + ' |' + Format(dbsas[x][0],5,3) + ' +/-' + Format(dbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') output.write('\n\n') output.close() return 0 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def decoutput34(output, debug, verbose, frames, comprm, recprm, ligprm, ligstart, type, decomp, surften, outname): import sys, math, string, os frames = float(frames) mdout = open("_MMPBSA_complex_" + type + ".mdout",'r') first = True # Logical variable -- is this the first frame we've analyzed? intot = False # if we are currently analyzing the total energy decomp inbac = False # if we are currently analyzing the backbone energy decomp insc = False # if we are currently analyzing the sidechain energy decomp skip = False # if we've just seen "PRINT" and must now skip lines until we reach data onnum = 0 # ON residue NUMber -- keeps track of which residue we're extracting data from cress = [] # Complex RESidueS -- array that will contain all printed decomp residues in complex cress2 = [] # Complex RESidueS 2 -- array that will contain all printed decomp residues in complex # that Complex RESidueS interacts with # set up arrays and variables for the deltas, but only if we know where the ligand is if ligstart != -1: dtint = [] # Delta Total INTernal decomp dtvdw = [] # Delta Total VDW decomp dteel = [] # Delta Total EEL decomp dtpol = [] # Delta Total POLar solvation decomp dtsas = [] # Delta Total non-polar solvation decomp dsint = [] # Delta Sidechain INTernal decomp dsvdw = [] # Delta Sidechain VDW decomp dseel = [] # Delta Sidechain EEL decomp dspol = [] # Delta Sidechain POLar solvation decomp dssas = [] # Delta Sidechain non-polar solvation decomp dbint = [] # Delta Backbone INTernal decomp dbvdw = [] # Delta Backbone VDW decomp dbeel = [] # Delta Backbone EEL decomp dbpol = [] # Delta Backbone POLar solvation decomp dbsas = [] # Delta Backbone non-polar solvation decomp recdone = 0 # how many receptor residues have been used in differences -- used as an index ligdone = 0 # same as above, but for ligand # All arrays below are arrays of 2-element arrays. First element is a running sum for the average, second # element is a running sum of squares for the standard deviation ctint = [] # Complex Total INTernal decomposition ctvdw = [] # Complex Total VDW decomposition cteel = [] # Complex Total EEL (electrostatic) decomposition ctpol = [] # Complex Total POLar solvation decomposition ctsas = [] # Complex Total non-polar solvation decomposition csint = [] # Complex Sidechain INTernal decomposition csvdw = [] # Complex Sidechain VDW decomposition cseel = [] # Complex Sidechain EEL (elecrostatic) decomposition cspol = [] # Complex Sidechain POLar solvation decomposition cssas = [] # Complex Sidechain non-polar solvation decomposition cbint = [] # Complex Backbone INTernal decomposition cbvdw = [] # Complex Backbone VDW decomposition cbeel = [] # Complex Backbone EEL decomposition cbpol = [] # Complex Backbone POLar solvation decomposition cbsas = [] # Complex Backbone non-polar solvation decomposition for line in mdout: if not intot and not inbac and not insc: # If we are not currently analyzing any decomp output if not "PRINT PAIR" in line: continue elif "TOTAL" in line: # we've found the TOTAL decomp data intot = True # we are now in total skip = True # we must skip the filler-lines between PRINT DECOMP and actual data elif "SIDECHAIN" in line: insc = True # same as for TOTAL skip = True elif "BACKBONE" in line: inbac = True # same as for TOTAL skip = True if intot: if skip: # if we are still skipping lines... (to distinguish from blank lines at end of data) if not "TDC" in line: continue # skip ahead to next line if we haven't started analyzing data else: skip = False # we've found the data, no more skipping. if not "TDC" in line: # if we've reached a blank line now, we've finished the data onnum = 0 # we're now on residue 0 again, as we're going to a new data block intot = False # we're no longer in total block continue # we've done what we need to -- skip to next line words = line.split() # now we're in to analyzing data -- split line at whitespace if first: # if we're on the first frame, we must build the arrays cress.append(int(words[1].strip('->'))) cress2.append(int(words[2])) ctint.append([float(words[3]),float(words[3])**2]) ctvdw.append([float(words[4]),float(words[4])**2]) cteel.append([float(words[5]),float(words[5])**2]) ctpol.append([float(words[6]),float(words[6])**2]) ctsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: # otherwise add to the running sums/sums of squares ctint[onnum][0] += float(words[3]) ctvdw[onnum][0] += float(words[4]) cteel[onnum][0] += float(words[5]) ctpol[onnum][0] += float(words[6]) ctsas[onnum][0] += float(words[7]) * surften ctint[onnum][1] += float(words[3]) ** 2 ctvdw[onnum][1] += float(words[4]) ** 2 cteel[onnum][1] += float(words[5]) ** 2 ctpol[onnum][1] += float(words[6]) ** 2 ctsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if insc: # same logic as we applied to TOTAL if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: csint.append([float(words[3]),float(words[3])**2]) csvdw.append([float(words[4]),float(words[4])**2]) cseel.append([float(words[5]),float(words[5])**2]) cspol.append([float(words[6]),float(words[6])**2]) cssas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: csint[onnum][0] += float(words[3]) csvdw[onnum][0] += float(words[4]) cseel[onnum][0] += float(words[5]) cspol[onnum][0] += float(words[6]) cssas[onnum][0] += float(words[7]) * surften csint[onnum][1] += float(words[3]) ** 2 csvdw[onnum][1] += float(words[4]) ** 2 cseel[onnum][1] += float(words[5]) ** 2 cspol[onnum][1] += float(words[6]) ** 2 cssas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if inbac: # same logic as we applied to TOTAL with one change... if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False # backbone is the last block to be printed, we now cannot be in first frame anymore continue words = line.split() if first: cbint.append([float(words[3]),float(words[3])**2]) cbvdw.append([float(words[4]),float(words[4])**2]) cbeel.append([float(words[5]),float(words[5])**2]) cbpol.append([float(words[6]),float(words[6])**2]) cbsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: cbint[onnum][0] += float(words[3]) cbvdw[onnum][0] += float(words[4]) cbeel[onnum][0] += float(words[5]) cbpol[onnum][0] += float(words[6]) cbsas[onnum][0] += float(words[7]) * surften cbint[onnum][1] += float(words[3]) ** 2 cbvdw[onnum][1] += float(words[4]) ** 2 cbeel[onnum][1] += float(words[5]) ** 2 cbpol[onnum][1] += float(words[6]) ** 2 cbsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(cress)): # we now compute averages and standard deviations ctint[x][0] /= frames ctint[x][1] = math.sqrt(abs(ctint[x][1]/frames - ctint[x][0] ** 2)) ctvdw[x][0] /= frames ctvdw[x][1] = math.sqrt(abs(ctvdw[x][1]/frames - ctvdw[x][0] ** 2)) cteel[x][0] /= frames cteel[x][1] = math.sqrt(abs(cteel[x][1]/frames - cteel[x][0] ** 2)) ctpol[x][0] /= frames ctpol[x][1] = math.sqrt(abs(ctpol[x][1]/frames - ctpol[x][0] ** 2)) ctsas[x][0] /= frames ctsas[x][1] = math.sqrt(abs(ctsas[x][1]/frames - ctsas[x][0] ** 2)) csint[x][0] /= frames csint[x][1] = math.sqrt(abs(csint[x][1]/frames - csint[x][0] ** 2)) csvdw[x][0] /= frames csvdw[x][1] = math.sqrt(abs(csvdw[x][1]/frames - csvdw[x][0] ** 2)) cseel[x][0] /= frames cseel[x][1] = math.sqrt(abs(cseel[x][1]/frames - cseel[x][0] ** 2)) cspol[x][0] /= frames cspol[x][1] = math.sqrt(abs(cspol[x][1]/frames - cspol[x][0] ** 2)) cssas[x][0] /= frames cssas[x][1] = math.sqrt(abs(cssas[x][1]/frames - cssas[x][0] ** 2)) cbint[x][0] /= frames cbint[x][1] = math.sqrt(abs(cbint[x][1]/frames - cbint[x][0] ** 2)) cbvdw[x][0] /= frames cbvdw[x][1] = math.sqrt(abs(cbvdw[x][1]/frames - cbvdw[x][0] ** 2)) cbeel[x][0] /= frames cbeel[x][1] = math.sqrt(abs(cbeel[x][1]/frames - cbeel[x][0] ** 2)) cbpol[x][0] /= frames cbpol[x][1] = math.sqrt(abs(cbpol[x][1]/frames - cbpol[x][0] ** 2)) cbsas[x][0] /= frames cbsas[x][1] = math.sqrt(abs(cbsas[x][1]/frames - cbsas[x][0] ** 2)) # Done with complex. Do the same thing for the receptor, only name arrays starting with 'r' # instead of 'c' for Receptor. mdout = open("_MMPBSA_receptor_" + type + ".mdout",'r') first = True intot = False inbac = False insc = False skip = False onnum = 0 rress = [] rress2 = [] rtint = [] rtvdw = [] rteel = [] rtpol = [] rtsas = [] rsint = [] rsvdw = [] rseel = [] rspol = [] rssas = [] rbint = [] rbvdw = [] rbeel = [] rbpol = [] rbsas = [] for line in mdout: if not intot and not inbac and not insc: if not "PRINT PAIR" in line: continue elif "TOTAL" in line: intot = True skip = True elif "SIDECHAIN" in line: insc = True skip = True elif "BACKBONE" in line: inbac = True skip = True if intot: if skip: if not "TDC" in line: continue else: skip = False if not "TDC" in line: onnum = 0 intot = False continue words = line.split() if first: rress.append(int(words[1].strip('->'))) rress2.append(int(words[2])) rtint.append([float(words[3]),float(words[3])**2]) rtvdw.append([float(words[4]),float(words[4])**2]) rteel.append([float(words[5]),float(words[5])**2]) rtpol.append([float(words[6]),float(words[6])**2]) rtsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: rtint[onnum][0] += float(words[3]) rtvdw[onnum][0] += float(words[4]) rteel[onnum][0] += float(words[5]) rtpol[onnum][0] += float(words[6]) rtsas[onnum][0] += float(words[7]) * surften rtint[onnum][1] += float(words[3]) ** 2 rtvdw[onnum][1] += float(words[4]) ** 2 rteel[onnum][1] += float(words[5]) ** 2 rtpol[onnum][1] += float(words[6]) ** 2 rtsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if insc: if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: rsint.append([float(words[3]),float(words[3])**2]) rsvdw.append([float(words[4]),float(words[4])**2]) rseel.append([float(words[5]),float(words[5])**2]) rspol.append([float(words[6]),float(words[6])**2]) rssas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: rsint[onnum][0] += float(words[3]) rsvdw[onnum][0] += float(words[4]) rseel[onnum][0] += float(words[5]) rspol[onnum][0] += float(words[6]) rssas[onnum][0] += float(words[7]) * surften rsint[onnum][1] += float(words[3]) ** 2 rsvdw[onnum][1] += float(words[4]) ** 2 rseel[onnum][1] += float(words[5]) ** 2 rspol[onnum][1] += float(words[6]) ** 2 rssas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if inbac: if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False continue words = line.split() if first: rbint.append([float(words[3]),float(words[3])**2]) rbvdw.append([float(words[4]),float(words[4])**2]) rbeel.append([float(words[5]),float(words[5])**2]) rbpol.append([float(words[6]),float(words[6])**2]) rbsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: rbint[onnum][0] += float(words[3]) rbvdw[onnum][0] += float(words[4]) rbeel[onnum][0] += float(words[5]) rbpol[onnum][0] += float(words[6]) rbsas[onnum][0] += float(words[7]) * surften rbint[onnum][1] += float(words[3]) ** 2 rbvdw[onnum][1] += float(words[4]) ** 2 rbeel[onnum][1] += float(words[5]) ** 2 rbpol[onnum][1] += float(words[6]) ** 2 rbsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(rress)): rtint[x][0] /= frames rtint[x][1] = math.sqrt(abs(rtint[x][1]/frames - rtint[x][0] ** 2)) rtvdw[x][0] /= frames rtvdw[x][1] = math.sqrt(abs(rtvdw[x][1]/frames - rtvdw[x][0] ** 2)) rteel[x][0] /= frames rteel[x][1] = math.sqrt(abs(rteel[x][1]/frames - rteel[x][0] ** 2)) rtpol[x][0] /= frames rtpol[x][1] = math.sqrt(abs(rtpol[x][1]/frames - rtpol[x][0] ** 2)) rtsas[x][0] /= frames rtsas[x][1] = math.sqrt(abs(rtsas[x][1]/frames - rtsas[x][0] ** 2)) rsint[x][0] /= frames rsint[x][1] = math.sqrt(abs(rsint[x][1]/frames - rsint[x][0] ** 2)) rsvdw[x][0] /= frames rsvdw[x][1] = math.sqrt(abs(rsvdw[x][1]/frames - rsvdw[x][0] ** 2)) rseel[x][0] /= frames rseel[x][1] = math.sqrt(abs(rseel[x][1]/frames - rseel[x][0] ** 2)) rspol[x][0] /= frames rspol[x][1] = math.sqrt(abs(rspol[x][1]/frames - rspol[x][0] ** 2)) rssas[x][0] /= frames rssas[x][1] = math.sqrt(abs(rssas[x][1]/frames - rssas[x][0] ** 2)) rbint[x][0] /= frames rbint[x][1] = math.sqrt(abs(rbint[x][1]/frames - rbint[x][0] ** 2)) rbvdw[x][0] /= frames rbvdw[x][1] = math.sqrt(abs(rbvdw[x][1]/frames - rbvdw[x][0] ** 2)) rbeel[x][0] /= frames rbeel[x][1] = math.sqrt(abs(rbeel[x][1]/frames - rbeel[x][0] ** 2)) rbpol[x][0] /= frames rbpol[x][1] = math.sqrt(abs(rbpol[x][1]/frames - rbpol[x][0] ** 2)) rbsas[x][0] /= frames rbsas[x][1] = math.sqrt(abs(rbsas[x][1]/frames - rbsas[x][0] ** 2)) # Same as receptor above mdout = open("_MMPBSA_ligand_" + type + ".mdout",'r') first = True intot = False inbac = False insc = False skip = False onnum = 0 lress = [] lress2 = [] ltint = [] ltvdw = [] lteel = [] ltpol = [] ltsas = [] lsint = [] lsvdw = [] lseel = [] lspol = [] lssas = [] lbint = [] lbvdw = [] lbeel = [] lbpol = [] lbsas = [] for line in mdout: if not intot and not inbac and not insc: if not "PRINT PAIR" in line: continue elif "TOTAL" in line: intot = True skip = True elif "SIDECHAIN" in line: insc = True skip = True elif "BACKBONE" in line: inbac = True skip = True if intot: if skip: if not "TDC" in line: continue else: skip = False if not "TDC" in line: onnum = 0 intot = False continue words = line.split() if first: lress.append(int(words[1].strip('->'))) lress2.append(int(words[2])) ltint.append([float(words[3]),float(words[3])**2]) ltvdw.append([float(words[4]),float(words[4])**2]) lteel.append([float(words[5]),float(words[5])**2]) ltpol.append([float(words[6]),float(words[6])**2]) ltsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: ltint[onnum][0] += float(words[3]) ltvdw[onnum][0] += float(words[4]) lteel[onnum][0] += float(words[5]) ltpol[onnum][0] += float(words[6]) ltsas[onnum][0] += float(words[7]) * surften ltint[onnum][1] += float(words[3]) ** 2 ltvdw[onnum][1] += float(words[4]) ** 2 lteel[onnum][1] += float(words[5]) ** 2 ltpol[onnum][1] += float(words[6]) ** 2 ltsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if insc: if skip: if not "SDC" in line: continue else: skip = False if not "SDC" in line: onnum = 0 insc = False continue words = line.split() if first: lsint.append([float(words[3]),float(words[3])**2]) lsvdw.append([float(words[4]),float(words[4])**2]) lseel.append([float(words[5]),float(words[5])**2]) lspol.append([float(words[6]),float(words[6])**2]) lssas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: lsint[onnum][0] += float(words[3]) lsvdw[onnum][0] += float(words[4]) lseel[onnum][0] += float(words[5]) lspol[onnum][0] += float(words[6]) lssas[onnum][0] += float(words[7]) * surften lsint[onnum][1] += float(words[3]) ** 2 lsvdw[onnum][1] += float(words[4]) ** 2 lseel[onnum][1] += float(words[5]) ** 2 lspol[onnum][1] += float(words[6]) ** 2 lssas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 if inbac: if skip: if not "BDC" in line: continue else: skip = False if not "BDC" in line: onnum = 0 inbac = False first = False continue words = line.split() if first: lbint.append([float(words[3]),float(words[3])**2]) lbvdw.append([float(words[4]),float(words[4])**2]) lbeel.append([float(words[5]),float(words[5])**2]) lbpol.append([float(words[6]),float(words[6])**2]) lbsas.append([float(words[7])*surften,(float(words[7])*surften)**2]) else: lbint[onnum][0] += float(words[3]) lbvdw[onnum][0] += float(words[4]) lbeel[onnum][0] += float(words[5]) lbpol[onnum][0] += float(words[6]) lbsas[onnum][0] += float(words[7]) * surften lbint[onnum][1] += float(words[3]) ** 2 lbvdw[onnum][1] += float(words[4]) ** 2 lbeel[onnum][1] += float(words[5]) ** 2 lbpol[onnum][1] += float(words[6]) ** 2 lbsas[onnum][1] += float(words[7]) ** 2 * surften ** 2 onnum += 1 mdout.close() for x in range(len(lress)): ltint[x][0] /= frames ltint[x][1] = math.sqrt(abs(ltint[x][1]/frames - ltint[x][0] ** 2)) ltvdw[x][0] /= frames ltvdw[x][1] = math.sqrt(abs(ltvdw[x][1]/frames - ltvdw[x][0] ** 2)) lteel[x][0] /= frames lteel[x][1] = math.sqrt(abs(lteel[x][1]/frames - lteel[x][0] ** 2)) ltpol[x][0] /= frames ltpol[x][1] = math.sqrt(abs(ltpol[x][1]/frames - ltpol[x][0] ** 2)) ltsas[x][0] /= frames ltsas[x][1] = math.sqrt(abs(ltsas[x][1]/frames - ltsas[x][0] ** 2)) lsint[x][0] /= frames lsint[x][1] = math.sqrt(abs(lsint[x][1]/frames - lsint[x][0] ** 2)) lsvdw[x][0] /= frames lsvdw[x][1] = math.sqrt(abs(lsvdw[x][1]/frames - lsvdw[x][0] ** 2)) lseel[x][0] /= frames lseel[x][1] = math.sqrt(abs(lseel[x][1]/frames - lseel[x][0] ** 2)) lspol[x][0] /= frames lspol[x][1] = math.sqrt(abs(lspol[x][1]/frames - lspol[x][0] ** 2)) lssas[x][0] /= frames lssas[x][1] = math.sqrt(abs(lssas[x][1]/frames - lssas[x][0] ** 2)) lbint[x][0] /= frames lbint[x][1] = math.sqrt(abs(lbint[x][1]/frames - lbint[x][0] ** 2)) lbvdw[x][0] /= frames lbvdw[x][1] = math.sqrt(abs(lbvdw[x][1]/frames - lbvdw[x][0] ** 2)) lbeel[x][0] /= frames lbeel[x][1] = math.sqrt(abs(lbeel[x][1]/frames - lbeel[x][0] ** 2)) lbpol[x][0] /= frames lbpol[x][1] = math.sqrt(abs(lbpol[x][1]/frames - lbpol[x][0] ** 2)) lbsas[x][0] /= frames lbsas[x][1] = math.sqrt(abs(lbsas[x][1]/frames - lbsas[x][0] ** 2)) # get the list of all residue names in the complex, receptor, and ligand so we can label the printed residues in the output com_resnames = getallresinfo(comprm,"RESIDUE_LABEL") rec_resnames = getallresinfo(recprm,"RESIDUE_LABEL") lig_resnames = getallresinfo(ligprm,"RESIDUE_LABEL") # print output if decomp == 3: output.write("idecomp = 3: Pairwise decomposition adding 1-4 interactions added to Internal.\n") if decomp == 4: output.write("idecomp = 4: Pairwise decomposition adding 1-4 interactions added to EEL and VDW.\n") if ligstart == -1 and verbose <= 2: output.write("DELTAs can only be printed with the default mask assignments\n\n") verbose += 2 if type == 'gb': output.write("Pairwise Energy Decomposition Analysis (All units kcal/mol): Generalized Born solvent\n\n") elif type == 'pb': output.write("Pairwise Energy Decomposition Analysis (All units kcal/mol): Poisson Boltzmann solvent\n\n") # Everything else will hereforth be added to the final output file using command-line >>. Close it here. output.close() # To avoid having to loop through several times, we use several temporary files. # tcom : file that holds Total COMplex : Opened for DEC_VERBOSE 2,3 # scom : file that holds Sidechain COMplex : Opened for DEC_VERBOSE 3 # bcom : file that holds Backbone COMplex : Opened for DEC_VERBOSE 3 # trec : file that holds Total RECeptor : Opened for DEC_VERBOSE 2,3 # srec : file that holds Sidechain RECeptor : Opened for DEC_VERBOSE 3 # brec : file that holds Backbone RECeptor : Opened for DEC_VERBOSE 3 # tlig : file that holds Total LIGand : Opened for DEC_VERBOSE 2,3 # slig : file that holds Sidechain LIGand : Opened for DEC_VERBOSE 3 # blig : file that holds Backbone LIGand : Opened for DEC_VERBOSE 3 # tdel : file that holds Total DELta : Opened for DEC_VERBOSE 0,1,2,3 # sdel : file that holds Sidechain DELta : Opened for DEC_VERBOSE 1,3 # bdel : file that holds Backbone DELta : Opened for DEC_VERBOSE 1,3 tdel = open('_MMPBSA_tdel.txt','w') if verbose == 2 or verbose == 3: tcom = open('_MMPBSA_tcom.txt','w') trec = open('_MMPBSA_trec.txt','w') tlig = open('_MMPBSA_tlig.txt','w') if verbose == 3: scom = open('_MMPBSA_scom.txt','w') bcom = open('_MMPBSA_bcom.txt','w') srec = open('_MMPBSA_srec.txt','w') brec = open('_MMPBSA_brec.txt','w') slig = open('_MMPBSA_slig.txt','w') blig = open('_MMPBSA_blig.txt','w') if verbose == 1 or verbose == 3: sdel = open('_MMPBSA_sdel.txt','w') bdel = open('_MMPBSA_bdel.txt','w') tdel.write(""" DELTAS: Total Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) if verbose == 2 or verbose == 3: tcom.write("Complex:" + """ Total Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) trec.write("\nReceptor:" + """ Total Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) tlig.write("\nLigand:" + """ Total Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) if verbose == 3: scom.write(""" Sidechain Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) bcom.write(""" Backbone Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) srec.write(""" Sidechain Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) brec.write(""" Backbone Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) slig.write(""" Sidechain Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) blig.write(""" Backbone Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) if verbose == 1 or verbose == 3: sdel.write(""" Sidechain Energy Decomposition: Residue | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) bdel.write(""" Backbone Energy Decomposition: Resid 1 | Resid 2 | Internal | van der Waals | Electrostatic | Polar Solvation | Non-Polar Solv. | TOTAL ----------------------------------------------------------------------------------------------------------------------------------------------------- """) for x in range(len(cress)): # master loop -- loop through once (maybe twice, again for deltas) if verbose == 2 or verbose == 3: # total complex total_avg = ctint[x][0] + ctvdw[x][0] + cteel[x][0] + ctpol[x][0] + ctsas[x][0] total_stdev = math.sqrt(ctint[x][1]**2 + ctvdw[x][1]**2 + cteel[x][1]**2 + ctpol[x][1]**2 + ctsas[x][1]**2) tcom.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(ctint[x][0],5,3) + ' +/-' + Format(ctint[x][1],3,3) + ' |' + Format(ctvdw[x][0],5,3) + ' +/-' + Format(ctvdw[x][1],3,3) + ' |' + Format(cteel[x][0],5,3) + ' +/-' + Format(cteel[x][1],3,3) + ' |' + Format(ctpol[x][0],5,3) + ' +/-' + Format(ctpol[x][1],3,3) + ' |' + Format(ctsas[x][0],5,3) + ' +/-' + Format(ctsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3: # sidechain complex total_avg = csint[x][0] + csvdw[x][0] + cseel[x][0] + cspol[x][0] + cssas[x][0] total_stdev = math.sqrt(csint[x][1]**2 + csvdw[x][1]**2 + cseel[x][1]**2 + cspol[x][1]**2 + cssas[x][1]**2) scom.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(csint[x][0],5,3) + ' +/-' + Format(csint[x][1],3,3) + ' |' + Format(csvdw[x][0],5,3) + ' +/-' + Format(csvdw[x][1],3,3) + ' |' + Format(cseel[x][0],5,3) + ' +/-' + Format(cseel[x][1],3,3) + ' |' + Format(cspol[x][0],5,3) + ' +/-' + Format(cspol[x][1],3,3) + ' |' + Format(cssas[x][0],5,3) + ' +/-' + Format(cssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # backbone complex total_avg = cbint[x][0] + cbvdw[x][0] + cbeel[x][0] + cbpol[x][0] + cbsas[x][0] total_stdev = math.sqrt(cbint[x][1]**2 + cbvdw[x][1]**2 + cbeel[x][1]**2 + cbpol[x][1]**2 + cbsas[x][1]**2) bcom.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(cbint[x][0],5,3) + ' +/-' + Format(cbint[x][1],3,3) + ' |' + Format(cbvdw[x][0],5,3) + ' +/-' + Format(cbvdw[x][1],3,3) + ' |' + Format(cbeel[x][0],5,3) + ' +/-' + Format(cbeel[x][1],3,3) + ' |' + Format(cbpol[x][0],5,3) + ' +/-' + Format(cbpol[x][1],3,3) + ' |' + Format(cbsas[x][0],5,3) + ' +/-' + Format(cbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 2 or verbose == 3 and x < len(rress): # total receptor total_avg = rtint[x][0] + rtvdw[x][0] + rteel[x][0] + rtpol[x][0] + rtsas[x][0] total_stdev = math.sqrt(rtint[x][1]**2 + rtvdw[x][1]**2 + rteel[x][1]**2 + rtpol[x][1]**2 + rtsas[x][1]**2) trec.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' | ' + rec_resnames[rress2[x]-1] + str(rress2[x]).rjust(4) + ' |' + Format(rtint[x][0],5,3) + ' +/-' + Format(rtint[x][1],3,3) + ' |' + Format(rtvdw[x][0],5,3) + ' +/-' + Format(rtvdw[x][1],3,3) + ' |' + Format(rteel[x][0],5,3) + ' +/-' + Format(rteel[x][1],3,3) + ' |' + Format(rtpol[x][0],5,3) + ' +/-' + Format(rtpol[x][1],3,3) + ' |' + Format(rtsas[x][0],5,3) + ' +/-' + Format(rtsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3 and x < len(rress): # sidechain receptor total_avg = rsint[x][0] + rsvdw[x][0] + rseel[x][0] + rspol[x][0] + rssas[x][0] total_stdev = math.sqrt(rsint[x][1]**2 + rsvdw[x][1]**2 + rseel[x][1]**2 + rspol[x][1]**2 + rssas[x][1]**2) srec.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' | ' + rec_resnames[rress2[x]-1] + str(rress2[x]).rjust(4) + ' |' + Format(rsint[x][0],5,3) + ' +/-' + Format(rsint[x][1],3,3) + ' |' + Format(rsvdw[x][0],5,3) + ' +/-' + Format(rsvdw[x][1],3,3) + ' |' + Format(rseel[x][0],5,3) + ' +/-' + Format(rseel[x][1],3,3) + ' |' + Format(rspol[x][0],5,3) + ' +/-' + Format(rspol[x][1],3,3) + ' |' + Format(rssas[x][0],5,3) + ' +/-' + Format(rssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # backbone receptor total_avg = rbint[x][0] + rbvdw[x][0] + rbeel[x][0] + rbpol[x][0] + rbsas[x][0] total_stdev = math.sqrt(rbint[x][1]**2 + rbvdw[x][1]**2 + rbeel[x][1]**2 + rbpol[x][1]**2 + rbsas[x][1]**2) brec.write( rec_resnames[rress[x]-1] + str(rress[x]).rjust(4) + ' | ' + rec_resnames[rress2[x]-1] + str(rress2[x]).rjust(4) + ' |' + Format(rbint[x][0],5,3) + ' +/-' + Format(rbint[x][1],3,3) + ' |' + Format(rbvdw[x][0],5,3) + ' +/-' + Format(rbvdw[x][1],3,3) + ' |' + Format(rbeel[x][0],5,3) + ' +/-' + Format(rbeel[x][1],3,3) + ' |' + Format(rbpol[x][0],5,3) + ' +/-' + Format(rbpol[x][1],3,3) + ' |' + Format(rbsas[x][0],5,3) + ' +/-' + Format(rbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 2 or verbose == 3 and x < len(lress): # total ligand total_avg = ltint[x][0] + ltvdw[x][0] + lteel[x][0] + ltpol[x][0] + ltsas[x][0] total_stdev = math.sqrt(ltint[x][1]**2 + ltvdw[x][1]**2 + lteel[x][1]**2 + ltpol[x][1]**2 + ltsas[x][1]**2) tlig.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' | ' + lig_resnames[lress2[x]-1] + str(lress2[x]).rjust(4) + ' |' + Format(ltint[x][0],5,3) + ' +/-' + Format(ltint[x][1],3,3) + ' |' + Format(ltvdw[x][0],5,3) + ' +/-' + Format(ltvdw[x][1],3,3) + ' |' + Format(lteel[x][0],5,3) + ' +/-' + Format(lteel[x][1],3,3) + ' |' + Format(ltpol[x][0],5,3) + ' +/-' + Format(ltpol[x][1],3,3) + ' |' + Format(ltsas[x][0],5,3) + ' +/-' + Format(ltsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 3 and x < len(lress): # sidechain ligand total_avg = lsint[x][0] + lsvdw[x][0] + lseel[x][0] + lspol[x][0] + lssas[x][0] total_stdev = math.sqrt(lsint[x][1]**2 + lsvdw[x][1]**2 + lseel[x][1]**2 + lspol[x][1]**2 + lssas[x][1]**2) slig.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' | ' + lig_resnames[lress2[x]-1] + str(lress2[x]).rjust(4) + ' |' + Format(lsint[x][0],5,3) + ' +/-' + Format(lsint[x][1],3,3) + ' |' + Format(lsvdw[x][0],5,3) + ' +/-' + Format(lsvdw[x][1],3,3) + ' |' + Format(lseel[x][0],5,3) + ' +/-' + Format(lseel[x][1],3,3) + ' |' + Format(lspol[x][0],5,3) + ' +/-' + Format(lspol[x][1],3,3) + ' |' + Format(lssas[x][0],5,3) + ' +/-' + Format(lssas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # backbone ligand total_avg = lbint[x][0] + lbvdw[x][0] + lbeel[x][0] + lbpol[x][0] + lbsas[x][0] total_stdev = math.sqrt(lbint[x][1]**2 + lbvdw[x][1]**2 + lbeel[x][1]**2 + lbpol[x][1]**2 + lbsas[x][1]**2) blig.write( lig_resnames[lress[x]-1] + str(lress[x]).rjust(4) + ' | ' + lig_resnames[lress2[x]-1] + str(lress2[x]).rjust(4) + ' |' + Format(lbint[x][0],5,3) + ' +/-' + Format(lbint[x][1],3,3) + ' |' + Format(lbvdw[x][0],5,3) + ' +/-' + Format(lbvdw[x][1],3,3) + ' |' + Format(lbeel[x][0],5,3) + ' +/-' + Format(lbeel[x][1],3,3) + ' |' + Format(lbpol[x][0],5,3) + ' +/-' + Format(lbpol[x][1],3,3) + ' |' + Format(lbsas[x][0],5,3) + ' +/-' + Format(lbsas[x][1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # Time for the fun part: DELTAs. # for x in range(len(cress)): # accumulate the deltas if ligstart != -1: # Figure out location of the complex residue if (cress[x] < ligstart or cress[x] >= ligstart + len(lig_resnames)) and \ (cress2[x] < ligstart or cress2[x] >= ligstart + len(lig_resnames)): # if both are in receptor dtint = [ctint[x][0] - rtint[recdone][0], math.sqrt(ctint[x][1]**2+rtint[recdone][1]**2)] dtvdw = [ctvdw[x][0] - rtvdw[recdone][0], math.sqrt(ctvdw[x][1]**2+rtvdw[recdone][1]**2)] dteel = [cteel[x][0] - rteel[recdone][0], math.sqrt(cteel[x][1]**2+rteel[recdone][1]**2)] dtpol = [ctpol[x][0] - rtpol[recdone][0], math.sqrt(ctpol[x][1]**2+rtpol[recdone][1]**2)] dtsas = [ctsas[x][0] - rtsas[recdone][0], math.sqrt(ctsas[x][1]**2+rtsas[recdone][1]**2)] dsint = [csint[x][0] - rsint[recdone][0], math.sqrt(csint[x][1]**2+rsint[recdone][1]**2)] dsvdw = [csvdw[x][0] - rsvdw[recdone][0], math.sqrt(csvdw[x][1]**2+rsvdw[recdone][1]**2)] dseel = [cseel[x][0] - rseel[recdone][0], math.sqrt(cseel[x][1]**2+rseel[recdone][1]**2)] dspol = [cspol[x][0] - rspol[recdone][0], math.sqrt(cspol[x][1]**2+rspol[recdone][1]**2)] dssas = [cssas[x][0] - rssas[recdone][0], math.sqrt(cssas[x][1]**2+rssas[recdone][1]**2)] dbint = [cbint[x][0] - rbint[recdone][0], math.sqrt(cbint[x][1]**2+rbint[recdone][1]**2)] dbvdw = [cbvdw[x][0] - rbvdw[recdone][0], math.sqrt(cbvdw[x][1]**2+rbvdw[recdone][1]**2)] dbeel = [cbeel[x][0] - rbeel[recdone][0], math.sqrt(cbeel[x][1]**2+rbeel[recdone][1]**2)] dbpol = [cbpol[x][0] - rbpol[recdone][0], math.sqrt(cbpol[x][1]**2+rbpol[recdone][1]**2)] dbsas = [cbsas[x][0] - rbsas[recdone][0], math.sqrt(cbsas[x][1]**2+rbsas[recdone][1]**2)] recdone += 1 elif (cress[x] >= ligstart and cress[x] < ligstart + len(lig_resnames)) and \ (cress2[x] >= ligstart and cress2[x] < ligstart + len(lig_resnames)): # if both are in ligand dtint = [ctint[x][0] - ltint[ligdone][0], math.sqrt(ctint[x][1]**2+ltint[ligdone][1]**2)] dtvdw = [ctvdw[x][0] - ltvdw[ligdone][0], math.sqrt(ctvdw[x][1]**2+ltvdw[ligdone][1]**2)] dteel = [cteel[x][0] - lteel[ligdone][0], math.sqrt(cteel[x][1]**2+lteel[ligdone][1]**2)] dtpol = [ctpol[x][0] - ltpol[ligdone][0], math.sqrt(ctpol[x][1]**2+ltpol[ligdone][1]**2)] dtsas = [ctsas[x][0] - ltsas[ligdone][0], math.sqrt(ctsas[x][1]**2+ltsas[ligdone][1]**2)] dsint = [csint[x][0] - lsint[ligdone][0], math.sqrt(csint[x][1]**2+lsint[ligdone][1]**2)] dsvdw = [csvdw[x][0] - lsvdw[ligdone][0], math.sqrt(csvdw[x][1]**2+lsvdw[ligdone][1]**2)] dseel = [cseel[x][0] - lseel[ligdone][0], math.sqrt(cseel[x][1]**2+lseel[ligdone][1]**2)] dspol = [cspol[x][0] - lspol[ligdone][0], math.sqrt(cspol[x][1]**2+lspol[ligdone][1]**2)] dssas = [cssas[x][0] - lssas[ligdone][0], math.sqrt(cssas[x][1]**2+lssas[ligdone][1]**2)] dbint = [cbint[x][0] - lbint[ligdone][0], math.sqrt(cbint[x][1]**2+lbint[ligdone][1]**2)] dbvdw = [cbvdw[x][0] - lbvdw[ligdone][0], math.sqrt(cbvdw[x][1]**2+lbvdw[ligdone][1]**2)] dbeel = [cbeel[x][0] - lbeel[ligdone][0], math.sqrt(cbeel[x][1]**2+lbeel[ligdone][1]**2)] dbpol = [cbpol[x][0] - lbpol[ligdone][0], math.sqrt(cbpol[x][1]**2+lbpol[ligdone][1]**2)] dbsas = [cbsas[x][0] - lbsas[ligdone][0], math.sqrt(cbsas[x][1]**2+lbsas[ligdone][1]**2)] ligdone += 1 else: # if one's in each, it's just in complex dtint = [ctint[x][0], ctint[x][1]] dtvdw = [ctvdw[x][0], ctvdw[x][1]] dteel = [cteel[x][0], cteel[x][1]] dtpol = [ctpol[x][0], ctpol[x][1]] dtsas = [ctsas[x][0], ctsas[x][1]] dsint = [csint[x][0], csint[x][1]] dsvdw = [csvdw[x][0], csvdw[x][1]] dseel = [cseel[x][0], cseel[x][1]] dspol = [cspol[x][0], cspol[x][1]] dssas = [cssas[x][0], cssas[x][1]] dbint = [cbint[x][0], cbint[x][1]] dbvdw = [cbvdw[x][0], cbvdw[x][1]] dbeel = [cbeel[x][0], cbeel[x][1]] dbpol = [cbpol[x][0], cbpol[x][1]] dbsas = [cbsas[x][0], cbsas[x][1]] total_avg = dtint[0] + dtvdw[0] + dteel[0] + dtpol[0] + dtsas[0] total_stdev = math.sqrt(dtint[1]**2 + dtvdw[1]**2 + dteel[1]**2 + dtpol[1]**2 + dtsas[1]**2) tdel.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(dtint[0],5,3) + ' +/-' + Format(dtint[1],3,3) + ' |' + Format(dtvdw[0],5,3) + ' +/-' + Format(dtvdw[1],3,3) + ' |' + Format(dteel[0],5,3) + ' +/-' + Format(dteel[1],3,3) + ' |' + Format(dtpol[0],5,3) + ' +/-' + Format(dtpol[1],3,3) + ' |' + Format(dtsas[0],5,3) + ' +/-' + Format(dtsas[1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') if verbose == 1 or verbose == 3: total_avg = dsint[0] + dsvdw[0] + dseel[0] + dspol[0] + dssas[0] total_stdev = math.sqrt(dsint[1]**2 + dsvdw[1]**2 + dseel[1]**2 + dspol[1]**2 + dssas[1]**2) sdel.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(dsint[0],5,3) + ' +/-' + Format(dsint[1],3,3) + ' |' + Format(dsvdw[0],5,3) + ' +/-' + Format(dsvdw[1],3,3) + ' |' + Format(dseel[0],5,3) + ' +/-' + Format(dseel[1],3,3) + ' |' + Format(dspol[0],5,3) + ' +/-' + Format(dspol[1],3,3) + ' |' + Format(dssas[0],5,3) + ' +/-' + Format(dssas[1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') total_avg = dbint[0] + dbvdw[0] + dbeel[0] + dbpol[0] + dbsas[0] total_stdev = math.sqrt(dbint[1]**2 + dbvdw[1]**2 + dbeel[1]**2 + dbpol[1]**2 + dbsas[1]**2) bdel.write( com_resnames[cress[x]-1] + str(cress[x]).rjust(4) + ' | ' + com_resnames[cress2[x]-1] + str(cress2[x]).rjust(4) + ' |' + Format(dbint[0],5,3) + ' +/-' + Format(dbint[1],3,3) + ' |' + Format(dbvdw[0],5,3) + ' +/-' + Format(dbvdw[1],3,3) + ' |' + Format(dbeel[0],5,3) + ' +/-' + Format(dbeel[1],3,3) + ' |' + Format(dbpol[0],5,3) + ' +/-' + Format(dbpol[1],3,3) + ' |' + Format(dbsas[0],5,3) + ' +/-' + Format(dbsas[1],3,3) + ' |' + Format(total_avg,5,3) + ' +/-' + Format(total_stdev,3,3) + '\n') # close all files tdel.close() if verbose == 2 or verbose == 3: tcom.close() trec.close() tlig.close() if verbose == 3: scom.close() bcom.close() srec.close() brec.close() slig.close() blig.close() if verbose == 1 or verbose == 3: sdel.close() bdel.close() if verbose == 2 or verbose == 3: os.system('cat _MMPBSA_tcom.txt >> ' + outname + '; rm _MMPBSA_tcom.txt') if verbose == 3: os.system('cat _MMPBSA_scom.txt >> ' + outname + '; rm _MMPBSA_scom.txt') os.system('cat _MMPBSA_bcom.txt >> ' + outname + '; rm _MMPBSA_bcom.txt') if verbose == 2 or verbose == 3: os.system('cat _MMPBSA_trec.txt >> ' + outname + '; rm _MMPBSA_trec.txt') if verbose == 3: os.system('cat _MMPBSA_srec.txt >> ' + outname + '; rm _MMPBSA_srec.txt') os.system('cat _MMPBSA_brec.txt >> ' + outname + '; rm _MMPBSA_brec.txt') if verbose == 2 or verbose == 3: os.system('cat _MMPBSA_tlig.txt >> ' + outname + '; rm _MMPBSA_tlig.txt') if verbose == 3: os.system('cat _MMPBSA_slig.txt >> ' + outname + '; rm _MMPBSA_slig.txt') os.system('cat _MMPBSA_blig.txt >> ' + outname + '; rm _MMPBSA_blig.txt') os.system('cat _MMPBSA_tdel.txt >> ' + outname + '; rm _MMPBSA_tdel.txt') if verbose == 1 or verbose == 3: os.system('cat _MMPBSA_sdel.txt >> ' + outname + '; rm _MMPBSA_sdel.txt') os.system('cat _MMPBSA_bdel.txt >> ' + outname + '; rm _MMPBSA_bdel.txt') os.system('cat >> {0} < EOF \n EOF'.format(outname)) return 0 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def separate(finaloutput): finaloutput.write('-------------------------------------------------------------------------------\n') finaloutput.write('-------------------------------------------------------------------------------\n') #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def GetPtrajEntropy(filename): import sys try: file = open(filename,'r') except IOError: print >> sys.stderr, 'Error: Entropy file {0} was not created!'.format(filename) return -1 results = [] # array to hold results: [total, translational, rotational, vibrational] rotationals = 0 # only the third place "rotational" occurs is the correct value -- only take the 3rd for line in file: # loop through all lines in the file words = line.split() # split along whitespace if len(words) > 0: # only analyze non-blank lines if words[0] == 'Total': results.append(words[3]) # the total contribution is the 4th word elif words[0] == 'translational': results.append(words[3]) # the translational contribution is the 4th word elif words[0] == 'rotational': # the rotational contribution is the 4th word rotationals = rotationals + 1 if rotationals % 3 == 0: # but only take the 3rd one results.append(words[3]) elif words[0] == 'vibrational': # get the vibrational entropy results.append(words[3]) if len(results) != 12: # if it didn't get trans, rot, vib, total from com, rec, lig: return error code print >> sys.stderr, 'Warning: Only {0} terms found in quasi-harmonic entropy output file!'.format(len(results)) file.close() # close file before return return -1 file.close() # close file before return return results #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def GetNmodeEntropy(warnings, ligand_prmtop, prefix, temp): import math, sys unconverged = False # check to see if any frames didn't converge warned = False # check to see if we've already warned about it, so we don't warn many times results = [] # results array comtotal = [] # array to hold total complex entropy for each frame rectotal = [] # array to hold total receptor entropy for each frame ligtotal = [] # array to hold total ligand entropy for each frame comtrans = [] # array to hold translational complex entropy for each frame rectrans = [] # array to hold translational receptor entropy for each frame ligtrans = [] # array to hold translational ligand entropy for each frame comrot = [] # array to hold rotational complex entropy for each frame recrot = [] # array to hold rotational receptor entropy for each frame ligrot = [] # array to hold rotational ligand entropy for each frame comvib = [] # array to hold vibrational complex entropy for each frame recvib = [] # array to hold vibrational receptor entropy for each frame ligvib = [] # array to hold vibrational ligand entropy for each frame try: mdout1 = open(prefix + 'complex_nm.out','r') # open all 3 entropy files mdout2 = open(prefix + 'receptor_nm.out','r') mdout3 = open(prefix + 'ligand_nm.out','r') except IOError: # if any fail, return error code and print warning message print >> sys.stderr, 'Error: _MMPBSA_******_nm.out files not created! Check output files.' return -1 for line in mdout1: if 'Total:' in line: words = line.split() comtotal.append(float(words[3])*temp/1000) # convert to kcal/mol elif 'translational:' in line: words = line.split() comtrans.append(float(words[3])*temp/1000) # convert to kcal/mol elif 'rotational:' in line: words = line.split() comrot.append(float(words[3])*temp/1000) # convert to kcal/mol elif 'vibrational:' in line: words = line.split() comvib.append(float(words[3])*temp/1000) # convert to kcal/mol elif '|---- Entropy not Calculated---|' in line and not warned: # warn if we haven't before print >> sys.stderr, 'Warning: Not all of your snapshots minimized within tolerable limits for nmode' warnings.append('Not all of your snapshots minimized within tolerable limits for nmode') unconverged = True # we know we haven't converged everyone warned = True # and that we've already been warned for line in mdout2: # do the same for the 2nd mdout if 'Total:' in line: words = line.split() rectotal.append(float(words[3])*temp/1000) elif 'translational:' in line: words = line.split() rectrans.append(float(words[3])*temp/1000) elif 'rotational:' in line: words = line.split() recrot.append(float(words[3])*temp/1000) elif 'vibrational:' in line: words = line.split() recvib.append(float(words[3])*temp/1000) elif '|---- Entropy not Calculated---|' in line and not warned: print >> sys.stderr, 'Warning: Not all of your snapshots minimized within tolerable limits for nmode' warnings.append('Not all of your snapshots minimized within tolerable limits for nmode') unconverged = True warned = True for line in mdout3: if 'Total:' in line: words = line.split() ligtotal.append(float(words[3])*temp/1000) elif 'translational:' in line: words = line.split() ligtrans.append(float(words[3])*temp/1000) elif 'rotational:' in line: words = line.split() ligrot.append(float(words[3])*temp/1000) elif 'vibrational:' in line: words = line.split() ligvib.append(float(words[3])*temp/1000) elif '|---- Entropy not Calculated---|' in line and not warned: print >> sys.stderr, 'Warning: Not all of your snapshots minimized within tolerable limits for nmode' warnings.append('Not all of your snapshots minimized within tolerable limits for nmode') unconverged = True warned = True # It appears that nmode does not correctly handle a ligand that is a # single atom. Thus, we calculate the entropy for this simple case here # Note that for a monatomic, there are no rotational/vibrational # degrees of freedom, so svib = srot = 0, and stot = strans. This can be # found from typical stat thermo text (i.e. McQuarrie) and will simply # be applied here. strans is volume (or concentration) dependent, so I # adopt the convention used in nmode to be consistent. The results have # been checked. if len(ligtotal) == 0: ligmass = float(getresinfo(1,ligand_prmtop,'MASS')) / 1000 / 6.0221415e23 transentropy = 1.987 * math.log(7.18678e45 * math.sqrt(ligmass) * ligmass) # now in cal/molK transentropy = transentropy * temp/1000 # to kcal/mol, use 300 K as T for x in range(len(comtotal)): ligrot.append(0) # there are no rotations for a single atom ligvib.append(0) ligtrans.append(transentropy) # translational entropy is total entropy, doesn't change ligtotal.append(transentropy) if not unconverged: # we can only subtract without index errors if everything converged deltas = [] for x in range(len(comtotal)): deltas.append(float(comtotal[x])-float(rectotal[x])-float(ligtotal[x])) results.append(average(comtotal)) # append all 12 averages (com/rec/lig by tot/trans/rot/vib) results.append(average(comtrans)) results.append(average(comrot)) results.append(average(comvib)) results.append(average(rectotal)) results.append(average(rectrans)) results.append(average(recrot)) results.append(average(recvib)) results.append(average(ligtotal)) results.append(average(ligtrans)) results.append(average(ligrot)) results.append(average(ligvib)) results.append(stdev(comtotal,results[0])) # append all 12 standard dev. (com/rec/lig by tot/trans/rot/vib) results.append(stdev(comtrans,results[1])) results.append(stdev(comrot,results[2])) results.append(stdev(comvib,results[3])) results.append(stdev(rectotal,results[4])) results.append(stdev(rectrans,results[5])) results.append(stdev(recrot,results[6])) results.append(stdev(recvib,results[7])) results.append(stdev(ligtotal,results[8])) results.append(stdev(ligtrans,results[9])) results.append(stdev(ligrot,results[10])) results.append(stdev(ligvib,results[11])) if not unconverged: results.append(average(deltas)) # if we have a list of deltas, get their averages results.append(stdev(deltas,results[24])) # and standard deviations else: results.append(results[0]-results[4]-results[8]) # otherwise subtract the averages and propagate stdevs. results.append(math.sqrt(results[12] * results[12] + results[16] * results[16] + results[20] * results[20])) return results #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def PrintPtrajEnt(file, finaloutput, string, temp): data = GetPtrajEntropy(file) # get ptraj entropy data if data == -1: # if error code was returned, return error code to calling function return -1 data.append(float(data[0])-float(data[4])-float(data[8])) data.append(float(data[1])-float(data[5])-float(data[9])) data.append(float(data[2])-float(data[6])-float(data[10])) data.append(float(data[3])-float(data[7])-float(data[11])) for x in range(len(data)): data[x] = float(data[x]) * temp / 1000 finaloutput.write(string + 'ENTROPY RESULTS (QUASI-HARMONIC APPROXIMATION) CALCULATED WITH PTRAJ:') finaloutput.write('\n\n Translational Rotational Vibrational Total') finaloutput.write('\nComplex: ' + Format(data[1],8,4) + Format(data[2],11,4) + Format(data[3],12,4) + Format(data[0],11,4)) finaloutput.write('\nReceptor: ' + Format(data[5],8,4) + Format(data[6],11,4) + Format(data[7],12,4) + Format(data[4],11,4)) finaloutput.write('\nLigand: ' + Format(data[9],8,4) + Format(data[10],11,4) + Format(data[11],12,4) + Format(data[8],11,4)) finaloutput.write('\n\nDELTA S: ' + Format(data[13],8,4) + Format(data[14],11,4) + Format(data[15],12,4) + Format(data[12],11,4)) finaloutput.write('\n') return data[12] #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def PrintNmodeEnt(finaloutput, warnings, ligand_prmtop, prefix, string, temp): import sys data = GetNmodeEntropy(warnings, ligand_prmtop, prefix, temp) for x in range(len(data)): # check that there is entropy data to print out if data[x] == -1: # if not, print error and return no entropy print >> sys.stderr, "Warning: No snapshots for nmode minimized within tolerable limits! Entropy not calculated." finaloutput.write("NO NMODE ENTROPY CALCULATED. SNAPSHOTS NOT MINIMIZED!\n\n") return [0,0] finaloutput.write(string + 'ENTROPY RESULTS (HARMONIC APPROXIMATION) CALCULATED WITH NMODE:\n') finaloutput.write('\nComplex:\nEntropy Term Average Std. Dev.') finaloutput.write('\n-----------------------------------------------------------') finaloutput.write('\nTranslational:' + Format(data[1],10,4) + Format(data[13],12,4)) finaloutput.write('\nRotational: ' + Format(data[2],10,4) + Format(data[14],12,4)) finaloutput.write('\nVibrational: ' + Format(data[3],10,4) + Format(data[15],12,4)) finaloutput.write('\nTotal: ' + Format(data[0],10,4) + Format(data[12],12,4) + '\n\n') finaloutput.write('\nReceptor:\nEntropy Term Average Std. Dev.') finaloutput.write('\n-----------------------------------------------------------') finaloutput.write('\nTranslational:' + Format(data[5],10,4) + Format(data[17],12,4)) finaloutput.write('\nRotational: ' + Format(data[6],10,4) + Format(data[18],12,4)) finaloutput.write('\nVibrational: ' + Format(data[7],10,4) + Format(data[19],12,4)) finaloutput.write('\nTotal: ' + Format(data[4],10,4) + Format(data[16],12,4) + '\n\n') finaloutput.write('\nLigand:\nEntropy Term Average Std. Dev.') finaloutput.write('\n-----------------------------------------------------------') finaloutput.write('\nTranslational:' + Format(data[9],10,4) + Format(data[21],12,4)) finaloutput.write('\nRotational: ' + Format(data[10],10,4) + Format(data[22],12,4)) finaloutput.write('\nVibrational: ' + Format(data[11],10,4) + Format(data[23],12,4)) finaloutput.write('\nTotal: ' + Format(data[8],10,4) + Format(data[20],12,4) + '\n\n') finaloutput.write('\nDELTA S total=' + Format(data[24],10,4) + ' +/- ' + Format(data[25],7,4)) finaloutput.write('\n') return [data[24],data[25]] # return total entropy average and st. dev. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def PrintFinalResults(output_file, inputfile_name, initial_traj, solvated_prmtop, complex_prmtop, receptor_mdcrd, solvated_receptor_prmtop, ligand_mdcrd, solvated_ligand_prmtop, receptor_prmtop, ligand_prmtop, mdcrd, mutant_complex_prmtop, mutant_receptor_prmtop, mutant_ligand_prmtop, maskholder, numres_lig, numframes, nmoderun, numframesnmode, pbrun, sanderpb, entropy, mutant_only, mutstring, warnings, keep_files, ala_entropy, alarun, gbrun, debug, one_trajectory, verbose, sander_apbs, decompout, idecomp, dec_verbose, ligstart, decomprun, surften, cavity_surften, temp): import os, sys, math os.system('echo "| Run on `date`" > ' + output_file) # Open FINAL_RESULTS_MMPBSA.dat and append to it, since the previous # line overwrites the file with the date in the first line finaloutput = open(output_file,'a') # Print the input file to the final output finaloutput.write('\n|Input file:\n') if inputfile_name == 'none': finaloutput.write('|No input file specified. Default values used for all variables.\n') else: finaloutput.write('|--------------------------------------------------------------\n') inpfile = open(inputfile_name,'r') for line in inpfile: finaloutput.write('|' + line) inpfile.close() finaloutput.write('|--------------------------------------------------------------\n') # Print out the user-input files for the system if initial_traj != 1: finaloutput.write('|Solvated complex topology file: ' + solvated_prmtop + '\n') finaloutput.write('|Complex topology file: ' + complex_prmtop + '\n') if len(receptor_mdcrd) != 0 and initial_traj != 1: finaloutput.write('|Solvated Receptor topology file: ' + solvated_receptor_prmtop + '\n') if len(ligand_mdcrd) != 0 and initial_traj != 1: finaloutput.write('|Solvated Ligand topology file: ' + solvated_ligand_prmtop + '\n') finaloutput.write('|Receptor topology file: ' + receptor_prmtop + '\n') finaloutput.write('|Ligand topology file: ' + ligand_prmtop) finaloutput.write('\n|Initial mdcrd(s): ') for x in range(len(mdcrd)): # print out every mdcrd that was used if x != 0: finaloutput.write('\n| ') finaloutput.write(mdcrd[x]) finaloutput.write('\n') if len(receptor_mdcrd) != 0: # print out every receptor_mdcrd that was used finaloutput.write('|Initial receptor mdcrd(s): ') for x in range(len(receptor_mdcrd)): if x != 0: finaloutput.write('\n| ') finaloutput.write(receptor_mdcrd[x]) finaloutput.write('\n') if len(ligand_mdcrd) != 0: # print out every ligand_mdcrd that was used finaloutput.write('|Initial ligand mdcrd(s): ') for x in range(len(ligand_mdcrd)): if x != 0: finaloutput.write('\n ') finaloutput.write(ligand_mdcrd[x]) finaloutput.write('\n') if alarun: # print out mutant prmtops used finaloutput.write('|Mutant complex topology file: ' + mutant_complex_prmtop + '\n') finaloutput.write('|Mutant receptor topology file: ' + mutant_receptor_prmtop + '\n') finaloutput.write('|Mutant ligand topology file: ' + mutant_ligand_prmtop + '\n') if maskholder != ['','']: # if we guessed the masks, print what we guessed finaloutput.write('|\n|Best guess for receptor mask: "' + maskholder[1] + '"\n') finaloutput.write('|Best guess for ligand mask: "' + maskholder[0] + '"\n') if numres_lig == 1: # if ligand is a single residue, print its name finaloutput.write('|Ligand residue name is "' + getresinfo(1,ligand_prmtop, "RESIDUE_LABEL") + '"\n|') # print out the number of frames used in the calculations (normal and nmode) finaloutput.write('\n|Calculations performed using {0} frames.\n'.format(numframes)) if nmoderun: finaloutput.write('|NMODE calculations performed using {0} frames.\n'.format(numframesnmode)) # print which pb solver was used, APBS or internal PBSA if pbrun: if sanderpb.endswith('APBS'): finaloutput.write('|Poisson Boltzmann calculations performed using iAPBS interface to sander (sander.APBS)\n') else: finaloutput.write('|Poisson Boltzmann calculations performed using internal PBSA solver in sander.\n') # state all units are kcal/mol finaloutput.write('|\n|All units are reported in kcal/mole.\n') if entropy == 1 or nmoderun: finaloutput.write('|All entropy results have units kcal/mole (Temperature is ' + str(temp) + ' K).\n') separate(finaloutput) ######################################################################## # # END write some details about the calculation to the final output file # ######################################################################## ######################################################################## # # BEGIN print out calculation results to the final output file # ######################################################################## # print out entropy results first. indices are ptraj_norm ptraj_mutant # nmode_avg_norm nmode_stdev_norm nmode_avg_mutant nmode_stdev_mutant # we do this so that we have the entropy values to find total DELTA G's later entropyvalues = [0,0,0,0,0,0] if entropy == 1: # if we do quasi-harmonic entropy with ptraj if mutant_only != 1: # if we aren't just doing mutant... entropy_holder = PrintPtrajEnt('_MMPBSA_ptraj_entropy.out',finaloutput,'',temp) # print normal ptraj entropy results if entropy_holder == -1: # if returned error code, quit remove(debug,0) entropyvalues[0] = entropy_holder # put the entropy holder into the first location in entropyvalues separate(finaloutput) if ala_entropy == 1 and alarun: # print out mutant entropy results if desired entropy_holder = PrintPtrajEnt('_MMPBSA_mutant_ptraj_entropy.out',finaloutput,mutstring + ' MUTANT ',temp) if entropy_holder == -1: remove(debug,0) entropyvalues[1] = entropy_holder separate(finaloutput) if mutant_only != 1: # if we're not just doing mutant, (but we did do mutant), show delta delta S binding. finaloutput.write('\nRESULT OF ALANINE SCANNING: ({0}) DELTA DELTA S binding = '.format(mutstring) + Format(entropyvalues[1]-entropyvalues[0],4,4) + '\n') separate(finaloutput) if nmoderun: # if we ran nmode calculation, print these next if mutant_only != 1: # if we're not just doing the mutant... isinerr = PrintNmodeEnt(finaloutput,warnings,ligand_prmtop,"_MMPBSA_",'', temp) if isinerr == -1: # non-fatal error print >> sys.stderr, 'Warning: nmode calculations failed! Check outputs.' warnings.append('nmode calculations failed! Check outputs.') keep_files = 2 # keep all files else: separate(finaloutput) entropyvalues[2] = isinerr[0] # put nmode entropy avg/stdev in positions 3,4 in array entropyvalues[3] = isinerr[1] if ala_entropy == 1 and alarun: # if we did nmode on the alanine mutant... finaloutput.write("\nMutant Entropy Results: ({0})\n".format(mutstring)) isinerr = PrintNmodeEnt(finaloutput,warnings,ligand_prmtop,"_MMPBSA_mutant_",mutstring + ' MUTANT ', temp) if isinerr == -1: # if there were problems with the calculation, non-fatal error print >> sys.stderr, 'Warning: nmode calculations failed! Check outputs.' warnings.append('nmode calculations failed! Check outputs.') keep_files = 2 # keep fall files if problems else: separate(finaloutput) entropyvalues[4] = isinerr[0] # put nmode entropy avg/stdev in positions 5,6 in array entropyvalues[5] = isinerr[1] if mutant_only != 1: # if we didn't just do mutant, find delta delta S binding finaloutput.write('\nRESULT OF ALANINE SCANNING: ({0}) DELTA DELTA S binding = '.format(mutstring) + Format(entropyvalues[4]-entropyvalues[2],4,4) + ' +/- ' + Format(math.sqrt(entropyvalues[3]**2 + entropyvalues[5]**2),4,4) + '\n') separate(finaloutput) # print out gb results if gbrun: if not alarun or mutant_only != 1: noalascan = gboutput('_MMPBSA_complex_gb.mdout','_MMPBSA_receptor_gb.mdout','_MMPBSA_ligand_gb.mdout', '',finaloutput, debug, numframes, one_trajectory,verbose) # noalascan is [avg,stdev] for total DELTA G if noalascan == -1: # if error code, quit remove(debug,0) separate(finaloutput) if entropy == 1: # if we did quasi-harmonic entropy approximation, print DELTA G with entropy finaloutput.write("Using Quasi-harmonic Entropy Approximation: DELTA G binding = " + Format((noalascan[0] - entropyvalues[0]),4,4) + '\n') separate(finaloutput) if nmoderun: # if we did nmode, print DELTA G with this entropy finaloutput.write("Using Normal Mode Entropy Approximation: DELTA G binding = " + Format((noalascan[0]- entropyvalues[2]),4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2+entropyvalues[3]**2),2,4) + '\n') separate(finaloutput) if alarun: # analyze alanine scanning: alascan is [avg,stdev] for total DELTA G for mutant alascan = gboutput('_MMPBSA_mutant_complex_gb.mdout','_MMPBSA_mutant_receptor_gb.mdout', '_MMPBSA_mutant_ligand_gb.mdout',mutstring,finaloutput,debug,numframes,one_trajectory,verbose) if alascan == -1: # if we returned error code, quit remove(debug,0) separate(finaloutput) if entropy == 1: # print DELTA G combined with ptraj entropy finaloutput.write("Using Quasi-harmonic Entropy Approximation: DELTA G binding = " + Format((alascan[0] - entropyvalues[1]),4,4) + '\n') separate(finaloutput) if nmoderun: # print DELTA G combined with nmode entropy finaloutput.write("Using Normal Mode Entropy Approximation: DELTA G binding = " + Format((alascan[0] - entropyvalues[4]),4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2 + entropyvalues[5]**2),2,4) + '\n') separate(finaloutput) if mutant_only != 1: # if we are not only analyzing mutant... finaloutput.write('\nRESULT OF ALANINE SCANNING: ({0}) DELTA DELTA G binding = '.format(mutstring) + Format( alascan[0]-noalascan[0],4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2 + alascan[1]**2),4,4) + '\n') if entropy == 1 and ala_entropy == 1: # print out DELTA DELTA G binding for quasi-harmonic approximation finaloutput.write('\n (quasi-harmonic entropy) ({0}) DELTA DELTA G binding = '.format(mutstring) + Format(alascan[0]-noalascan[0] - entropyvalues[1]+entropyvalues[0],4,4) + '\n') if nmoderun and ala_entropy == 1: # print out DELTA DELTA G binding for nmode approximation finaloutput.write('\n (normal mode entropy) ({0}) DELTA DELTA G binding = '.format(mutstring) + Format(alascan[0]-noalascan[0]-entropyvalues[1]+entropyvalues[0],4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2 + entropyvalues[5]**2 + alascan[1]**2 + entropyvalues[3]**2),4,4) + '\n') separate(finaloutput) # print out pb results if pbrun: if not alarun or mutant_only != 1: # print out normal pb results if not only analyzing mutant noalascan = pboutput('_MMPBSA_complex_pb.mdout','_MMPBSA_receptor_pb.mdout','_MMPBSA_ligand_pb.mdout', '',finaloutput,debug,numframes,sander_apbs,one_trajectory,verbose) if noalascan == -1: # if it returns error code, quit remove(debug,0) separate(finaloutput) if entropy == 1: # if we did quasi-harmonic approx, print DELTA G binding corrected with this finaloutput.write("Using Quasi-harmonic Entropy Approximation: DELTA G binding = " + Format((noalascan[0] - entropyvalues[0]),4,4) + '\n') separate(finaloutput) if nmoderun: # if we did nmode, print DELTA G binding corrected with this finaloutput.write("Using Normal Mode Entropy Approximation: DELTA G binding = " + Format((noalascan[0] - entropyvalues[2]),4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2+entropyvalues[3]**2),2,4) + '\n') separate(finaloutput) if alarun: # if we did alanine scanning, analyze these pb outputs alascan = pboutput('_MMPBSA_mutant_complex_pb.mdout','_MMPBSA_mutant_receptor_pb.mdout','_MMPBSA_mutant_ligand_pb.mdout', mutstring,finaloutput,mutstring,numframes,sander_apbs,one_trajectory,verbose) if alascan == -1: # if it returned an error code, quit remove(debug,0) separate(finaloutput) if ala_entropy == 1: # add entropy term from quasi-harmonic approximation if done finaloutput.write("Using Quasi-harmonic Entropy Approximation: DELTA G binding = " + Format((alascan[0]-entropyvalues[1]),4,4) + '\n') separate(finaloutput) if nmoderun: # add nmode entropy term if done finaloutput.write("Using Normal Mode Entropy Approximation: DELTA G binding = " + Format((alascan[0] - entropyvalues[4]),4,4) + ' +/- ' + Format(math.sqrt(noalascan[1]**2 + entropyvalues[5]**2),2,4) + '\n') separate(finaloutput) if not alarun or mutant_only != 1: # if we did both mutant and regular, print DELTA DELTA G binding. finaloutput.write('\nRESULT OF ALANINE SCANNING: ({0}) DELTA DELTA G binding = '.format(mutstring) + Format(alascan[0]-noalascan[0],4,4) +' +/- '+ Format(math.sqrt((noalascan[1]**2+alascan[1]**2)),4,4) + '\n') if entropy == 1 and ala_entropy == 1: # add entropy correction terms if done (quasi-harmonic) finaloutput.write('\n (quasi-harmonic entropy) ({0}) DELTA DELTA G binding = '.format(mutstring) + Format(alascan[0]-noalascan[0]-entropyvalues[1]+entropyvalues[0],4,4) + '\n') if nmoderun and ala_entropy == 1: # nmode entropy correction finaloutput.write('\n (normal mode entropy) ({0}) DELTA DELTA G binding = '.format(mutstring) + Format(alascan[0]-noalascan[0]-entropyvalues[1]+entropyvalues[0],4,4)+' +/- ' + Format(math.sqrt(noalascan[1]**2 +entropyvalues[5]**2 + alascan[1]**2 +entropyvalues[3]**2),4,4) + '\n') separate(finaloutput) # print out all warnings, now that they will all have been generated if len(warnings) > 0: finaloutput.write('\nWARNINGS:\n') for x in range(len(warnings)): finaloutput.write(warnings[x] + '\n') finaloutput.close() # Output decomp results if decomprun: os.system('echo "| Run on `date`" > ' + decompout) if gbrun: success = decomp_out(decompout, idecomp, debug, dec_verbose, 'gb', numframes, complex_prmtop, receptor_prmtop, ligand_prmtop, ligstart, surften) if pbrun: success = decomp_out(decompout, idecomp, debug, dec_verbose, 'pb', numframes, complex_prmtop, receptor_prmtop, ligand_prmtop, ligstart, cavity_surften) if success == -1: print >> sys.stderr, 'All files kept for debugging purposes' keep_files = 2 return keep_files #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++