Source Code for Module MolKit.pdb2pqr.src.psize

  1  """ psize class 
  2   
  3      Get dimensions and other information from a PQR file. 
  4   
  5      Originally written by Dave Sept 
  6      Additional APBS-specific features added by Nathan Baker 
  7      Ported to Python/Psize class by Todd Dolinsky and subsequently 
  8      hacked by Nathan Baker 
  9   
 10          ---------------------------- 
 11      
 12      PDB2PQR -- An automated pipeline for the setup, execution, and analysis of 
 13      Poisson-Boltzmann electrostatics calculations 
 14   
 15      Nathan A. Baker (baker@biochem.wustl.edu) 
 16      Todd Dolinsky (todd@ccb.wustl.edu) 
 17      Dept. of Biochemistry and Molecular Biophysics 
 18      Center for Computational Biology 
 19      Washington University in St. Louis 
 20   
 21      Jens Nielsen (Jens.Nielsen@ucd.ie) 
 22      University College Dublin 
 23   
 24      Additional contributing authors listed in documentation and supporting 
 25      package licenses. 
 26   
 27      Copyright (c) 2003-2007.  Washington University in St. Louis.   
 28      All Rights Reserved. 
 29   
 30      This file is part of PDB2PQR. 
 31   
 32      PDB2PQR is free software; you can redistribute it and/or modify 
 33      it under the terms of the GNU General Public License as published by 
 34      the Free Software Foundation; either version 2 of the License, or 
 35      (at your option) any later version. 
 36    
 37      PDB2PQR is distributed in the hope that it will be useful, 
 38      but WITHOUT ANY WARRANTY; without even the implied warranty of 
 39      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 40      GNU General Public License for more details. 
 41      
 42      You should have received a copy of the GNU General Public License 
 43      along with PDB2PQR; if not, write to the Free Software 
 44      Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA 
 45   
 46      ---------------------------- 
 47  """ 
 48   
 49  # User - Definable Variables: Default values 
 50   
 51  # CFAC = 1.7                  # Factor by which to expand mol dims to 
 52                                # get coarse grid dims 
 53  # FADD = 20                   # Amount to add to mol dims to get fine 
 54                                # grid dims 
 55  # SPACE = 0.50                # Desired fine mesh resolution 
 56  # GMEMFAC = 200               # Number of bytes per grid point required  
 57                                # for sequential MG calculation  
 58  # GMEMCEIL = 400              # Max MB allowed for sequential MG 
 59                                # calculation.  Adjust this to force the 
 60                                # script to perform faster calculations (which 
 61                                # require more parallelism). 
 62  # OFAC = 0.1                  # Overlap factor between mesh partitions 
 63  # REDFAC = 0.25               # The maximum factor by which a domain 
 64                                # dimension can be reduced during focusing 
 65  # TFAC_ALPHA = 9e-5           # Number of sec/unkown for setup/solve on 667 
 66                                # MHz EV67 Alpha CPU -- VERY ROUGH ESTIMATE 
 67  # TFAC_XEON  = 3e-4           # Number of sec/unkown for setup/solve on 500 
 68                                # MHz PIII Xeon CPU -- VERY ROUGH ESTIMATE 
 69  # TFAC_SPARC  = 5e-4          # Number of sec/unkown for setup/solve on 400 
 70                                # MHz UltraSPARC II CPU -- VERY ROUGH ESTIMATE 
 71   
 72  import string, sys 
 73  from sys import stdout 
 74  from math import log 
 75   
 76 -class Psize: 
 77 -    def __init__(self): 
 78          self.constants = {"CFAC":1.7, "FADD":20, "SPACE":0.50, "GMEMFAC":200, "GMEMCEIL":400, "OFAC":0.1, "REDFAC":0.25, "TFAC_ALPHA":9e-5, "TFAC_XEON":3e-4, "TFAC_SPARC": 5e-4} 
 79          self.minlen = [360.0, 360.0, 360.0] 
 80          self.maxlen = [0.0, 0.0, 0.0] 
 81          self.q = 0.0 
 82          self.gotatom = 0 
 83          self.gothet = 0 
 84          self.olen = [0.0, 0.0, 0.0] 
 85          self.cen = [0.0, 0.0, 0.0] 
 86          self.clen = [0.0, 0.0, 0.0] 
 87          self.flen = [0.0, 0.0, 0.0] 
 88          self.n = [0, 0, 0] 
 89          self.np = [0.0, 0.0, 0.0] 
 90          self.nsmall = [0,0,0] 
 91          self.nfocus = 0 
 92   
 93 -    def parseString(self, structure): 
 94          """ Parse the input structure as a string in PDB or PQR format """ 
 95          lines = string.split(structure, "\n") 
 96          self.parseLines(lines) 
 97   
 98 -    def parseInput(self, filename): 
 99          """ Parse input structure file in PDB or PQR format """ 
100          file = open(filename, "r") 
101          self.parseLines(file.readlines()) 
102   
103 -    def parseLines(self, lines): 
104          """ Parse the lines """ 
105          for line in lines: 
106              if string.find(line,"ATOM") == 0: 
107                  subline = string.replace(line[30:], "-", " -") 
108                  words = string.split(subline) 
109                  if len(words) < 4:     
110                      #sys.stderr.write("Can't parse following line:\n") 
111                      #sys.stderr.write("%s\n" % line) 
112                      #sys.exit(2) 
113                      continue 
114                  self.gotatom = self.gotatom + 1 
115                  self.q = self.q + float(words[3]) 
116                  rad = float(words[4]) 
117                  center = [] 
118                  for word in words[0:3]: 
119                      center.append(float(word)) 
120                  for i in range(3): 
121                      self.minlen[i] = min(center[i]-rad, self.minlen[i]) 
122                      self.maxlen[i] = max(center[i]+rad, self.maxlen[i]) 
123              elif string.find(line, "HETATM") == 0: 
124                  self.gothet = self.gothet + 1 
125       
126 -    def setConstant(self, name, value): 
127          """ Set a constant to a value; returns 0 if constant not found """ 
128          try: 
129              self.constants[name] = value 
130              return 1 
131          except KeyError: 
132              return 0 
133   
134 -    def getConstant(self, name): 
135          """ Get a constant value; raises KeyError if constant not found """ 
136          return self.constants[name] 
137   
138 -    def setLength(self, maxlen, minlen): 
139          """ Compute molecule dimensions """ 
140          for i in range(3): 
141              self.olen[i] = maxlen[i] - minlen[i] 
142              if self.olen[i] < 0.1: 
143                  self.olen[i] = 0.1 
144          return self.olen 
145   
146   
147 -    def setCoarseGridDims(self, olen): 
148          """ Compute coarse mesh dimensions """ 
149          for i in range(3): 
150              self.clen[i] = self.constants["CFAC"] * olen[i] 
151          return self.clen 
152   
153 -    def setFineGridDims(self, olen, clen): 
154          """ Compute fine mesh dimensions """ 
155          for i in range(3): 
156              self.flen[i] = olen[i] + self.constants["FADD"] 
157              if self.flen[i] > clen[i]: 
158                  #str = "WARNING: Fine length (%.2f) cannot be larger than coarse length (%.2f)\n" % (self.flen[i], clen[i]) 
159                  #str = str + "         Setting fine grid length equal to coarse grid length\n" 
160                  #stdout.write(str) 
161                  self.flen[i] = clen[i] 
162          return self.flen 
163   
164 -    def setCenter(self, maxlen, minlen): 
165          """ Compute molecule center """ 
166          for i in range(3): 
167              self.cen[i] = (maxlen[i] + minlen[i]) / 2 
168          return self.cen 
169   
170   
171 -    def setFineGridPoints(self, flen): 
172          """ Compute mesh grid points, assuming 4 levels in MG hierarchy """ 
173          tn = [0,0,0] 
174          for i in range(3): 
175              tn[i] = int(flen[i]/self.constants["SPACE"] + 0.5) 
176              self.n[i] = 32*(int((tn[i] - 1) / 32 + 0.5)) + 1 
177              if self.n[i] < 33: 
178                  self.n[i] = 33 
179          return self.n 
180   
181 -    def setSmallest(self, n): 
182          """ Compute parallel division in case memory requirement above ceiling 
183          Find the smallest dimension and see if the number of grid points in 
184          that dimension will fit below the memory ceiling 
185          Reduce nsmall until an nsmall^3 domain will fit into memory """ 
186          nsmall = [] 
187          for i in range(3): 
188              nsmall.append(n[i]) 
189          while 1: 
190              nsmem = 200.0 * nsmall[0] * nsmall[1] * nsmall[2] / 1024 / 1024 
191              if nsmem < self.constants["GMEMCEIL"]: break 
192              else: 
193                  i = nsmall.index(max(nsmall)) 
194                  nsmall[i] = 32 * ((nsmall[i] - 1)/32 - 1) + 1 
195                  if nsmall <= 0: 
196                      stdout.write("You picked a memory ceiling that is too small\n") 
197                      sys.exit(0)         
198   
199          self.nsmall = nsmall 
200          return nsmall 
201   
202 -    def setProcGrid(self, n, nsmall): 
203          """ Calculate the number of processors required to span each  
204          dimension """ 
205   
206          zofac = 1 + 2 * self.constants["OFAC"] 
207          for i in range(3): 
208              self.np[i] = n[i]/float(nsmall[i]) 
209              if self.np[i] > 1: self.np[i] = int(zofac*n[1]/nsmall[i] + 1.0) 
210          return self.np 
211                                                   
212 -    def setFocus(self, flen, np, clen): 
213          """ Calculate the number of levels of focusing required for each 
214          processor subdomain """ 
215   
216          nfoc = [0,0,0] 
217          for i in range(3): 
218              nfoc[i] = int(log((flen[i]/np[i])/clen[i])/log(self.constants["REDFAC"]) + 1.0) 
219          nfocus = nfoc[0] 
220          if nfoc[1] > nfocus: nfocus = nfoc[1] 
221          if nfoc[2] > nfocus: nfocus = nfoc[2] 
222          if nfocus > 0: nfocus = nfocus + 1 
223          self.nfocus = nfocus 
224   
225 -    def setAll(self): 
226          """ Set up all of the things calculated individually above """ 
227          maxlen = self.getMax() 
228          minlen = self.getMin() 
229          self.setLength(maxlen, minlen) 
230          olen = self.getLength() 
231           
232          self.setCoarseGridDims(olen) 
233          clen = self.getCoarseGridDims()         
234           
235          self.setFineGridDims(olen, clen) 
236          flen = self.getFineGridDims() 
237           
238          self.setCenter(maxlen, minlen) 
239          cen = self.getCenter() 
240           
241          self.setFineGridPoints(flen) 
242          n = self.getFineGridPoints() 
243           
244          self.setSmallest(n) 
245          nsmall = self.getSmallest() 
246           
247          self.setProcGrid(n, nsmall) 
248          np = self.getProcGrid() 
249           
250          self.setFocus(flen, np, clen) 
251          nfocus = self.getFocus() 
252           
253 -    def getMax(self): return self.maxlen 
254 -    def getMin(self): return self.minlen 
255 -    def getCharge(self): return self.q 
256 -    def getLength(self): return self.olen 
257 -    def getCoarseGridDims(self): return self.clen 
258 -    def getFineGridDims(self): return self.flen 
259 -    def getCenter(self): return self.cen 
260 -    def getFineGridPoints(self): return self.n 
261 -    def getSmallest(self): return self.nsmall 
262 -    def getProcGrid(self): return self.np 
263 -    def getFocus(self): return self.nfocus 
264   
265 -    def runPsize(self, filename): 
266          """ Parse input PQR file and set parameters """ 
267          self.parseInput(filename) 
268          self.setAll() 
269       
270 -    def printResults(self): 
271          """ Return a string with the formatted results """ 
272   
273          str = "\n" 
274           
275          if self.gotatom > 0: 
276   
277              maxlen = self.getMax() 
278              minlen = self.getMin() 
279              q = self.getCharge() 
280              olen = self.getLength() 
281              clen = self.getCoarseGridDims()         
282              flen = self.getFineGridDims() 
283              cen = self.getCenter() 
284              n = self.getFineGridPoints() 
285              nsmall = self.getSmallest() 
286              np = self.getProcGrid() 
287              nfocus = self.getFocus() 
288   
289              # Compute memory requirements 
290   
291              nsmem = 200.0 * nsmall[0] * nsmall[1] * nsmall[2] / 1024 / 1024 
292              gmem = 200.0 * n[0] * n[1] * n[2] / 1024 / 1024 
293               
294              # Calculate VERY ROUGH wall clock times 
295   
296              tsolve_alpha = nfocus*nsmall[0]*nsmall[1]*nsmall[2]*self.constants["TFAC_ALPHA"]; 
297              tsolve_xeon = nfocus*nsmall[0]*nsmall[1]*nsmall[2]*self.constants["TFAC_XEON"]; 
298              tsolve_sparc = nfocus*nsmall[0]*nsmall[1]*nsmall[2]*self.constants["TFAC_SPARC"]; 
299   
300              # Print the calculated entries 
301              str = str + "################# MOLECULE INFO ####################\n" 
302              str = str + "Number of ATOM entries = %i\n" % self.gotatom 
303              str = str + "Number of HETATM entries (ignored) = %i\n" % self.gothet 
304              str = str + "Total charge = %.3f e\n" % q 
305              str = str + "Dimensions = %.3f x %.3f x %.3f A\n" % (olen[0], olen[1], olen[2]) 
306              str = str + "Center = %.3f x %.3f x %.3f A\n" % (cen[0], cen[1], cen[2]) 
307              str = str + "Lower corner = %.3f x %.3f x %.3f A\n" % (minlen[0], minlen[1], minlen[2]) 
308              str = str + "Upper corner = %.3f x %.3f x %.3f A\n" % (maxlen[0], maxlen[1], maxlen[2]) 
309   
310              str = str + "\n" 
311              str = str + "############## GENERAL CALCULATION INFO #############\n" 
312              str = str + "Coarse grid dims = %.3f x %.3f x %.3f A\n" % (clen[0], 
313  clen[1], clen[2]) 
314              str = str + "Fine grid dims = %.3f x %.3f x %.3f A\n" % (flen[0], flen[1], flen[2]) 
315              str = str + "Num. fine grid pts. = %i x %i x %i\n" % (n[0], n[1], n[2]) 
316   
317              if gmem > self.constants["GMEMCEIL"]: 
318                  str = str + "Parallel solve required (%.3f MB > %.3f MB)\n" % (gmem, self.constants["GMEMCEIL"]) 
319                  str = str + "Total processors required = %i\n" % (np[0]*np[1]*np[2]) 
320                  str = str + "Proc. grid = %i x %i x %i\n" % (np[0], np[1], np[2]) 
321                  str = str + "Grid pts. on each proc. = %i x %i x %i\n" % (nsmall[0], nsmall[1], nsmall[2]) 
322                  xglob = np[0]*round(nsmall[0]/(1 + 2*self.constants["OFAC"]) - .001) 
323                  yglob = np[1]*round(nsmall[1]/(1 + 2*self.constants["OFAC"]) - .001) 
324                  zglob = np[2]*round(nsmall[2]/(1 + 2*self.constants["OFAC"]) - .001) 
325                  if np[0] == 1: xglob = nsmall[0] 
326                  if np[1] == 1: yglob = nsmall[1] 
327                  if np[2] == 1: zglob = nsmall[2] 
328                  str = str + "Fine mesh spacing = %g x %g x %g A\n" % (flen[0]/(xglob-1), flen[1]/(yglob-1), flen[2]/(zglob-1)) 
329                  str = str + "Estimated mem. required for parallel solve = %.3f MB/proc.\n" % nsmem 
330                  ntot = nsmall[0]*nsmall[1]*nsmall[2] 
331   
332              else: 
333                  str = str + "Fine mesh spacing = %g x %g x %g A\n" % (flen[0]/(n[0]-1), flen[1]/(n[1]-1), flen[2]/(n[2]-1)) 
334                  str = str + "Estimated mem. required for sequential solve = %.3f MB\n" % gmem 
335                  ntot = n[0]*n[1]*n[2] 
336   
337              str = str + "Number of focusing operations = %i\n" % nfocus 
338   
339              str = str + "\n" 
340              str = str + "################# ESTIMATED REQUIREMENTS ####################\n" 
341              str = str + "Memory per processor                   = %.3f MB\n" % (200.0*ntot/1024/1024) 
342              str = str + "Grid storage requirements (ASCII)      = %.3f MB\n" % (8.0*12*np[0]*np[1]*np[2]*ntot/1024/1024) 
343              str = str + "Grid storage requirements (XDR)        = %.3f MB\n" % (8.0*ntot*np[0]*np[1]*np[2]/1024/1024) 
344              str = str + "Time to solve on 667 MHz EV67 Alpha    = %.3f sec\n" % tsolve_alpha 
345              str = str + "Time to solve on 500 MHz PIII Xeon     = %.3f sec\n" % tsolve_xeon 
346              str = str + "Time to solve on 400 MHz UltraSparc II = %.3f sec\n" % tsolve_sparc 
347              str = str + "\n" 
348   
349          else: 
350              str = str + "No ATOM entires in file!\n\n" 
351   
352          return str 
353   
354 -def usage(): 
355      psize = Psize() 
356      usage = "\n" 
357      usage = usage + "Psize script\n" 
358      usage = usage + "Usage: psize.py [opts] <filename>\n" 
359      usage = usage + "Optional Arguments:\n" 
360      usage = usage + "  --help               : Display this text\n" 
361      usage = usage + "  --CFAC=<value>       : Factor by which to expand mol dims to\n" 
362      usage = usage + "                         get coarse grid dims\n" 
363      usage = usage + "                         [default = %g]\n" % psize.getConstant("CFAC") 
364      usage = usage + "  --FADD=<value>       : Amount to add to mol dims to get fine\n" 
365      usage = usage + "                         grid dims\n" 
366      usage = usage + "                         [default = %g]\n" % psize.getConstant("FADD") 
367      usage = usage + "  --SPACE=<value>      : Desired fine mesh resolution\n" 
368      usage = usage + "                         [default = %g]\n" % psize.getConstant("SPACE") 
369      usage = usage + "  --GMEMFAC=<value>    : Number of bytes per grid point required\n" 
370      usage = usage + "                         for sequential MG calculation\n" 
371      usage = usage + "                         [default = %g]\n" % psize.getConstant("GMEMFAC") 
372      usage = usage + "  --GMEMCEIL=<value>   : Max MB allowed for sequential MG\n" 
373      usage = usage + "                         calculation.  Adjust this to force the\n" 
374      usage = usage + "                         script to perform faster calculations (which\n" 
375      usage = usage + "                         require more parallelism).\n" 
376      usage = usage + "                         [default = %g]\n" % psize.getConstant("GMEMCEIL") 
377      usage = usage + "  --OFAC=<value>       : Overlap factor between mesh partitions\n" 
378      usage = usage + "                         [default = %g]\n" % psize.getConstant("OFAC") 
379      usage = usage + "  --REDFAC=<value>     : The maximum factor by which a domain\n" 
380      usage = usage + "                         dimension can be reduced during focusing\n" 
381      usage = usage + "                         [default = %g]\n" % psize.getConstant("REDFAC") 
382      usage = usage + "  --TFAC_ALPHA=<value> : Number of sec/unknown for setup/solve on 667\n" 
383      usage = usage + "                         MHz EV67 Alpha CPU -- VERY ROUGH ESTIMATE\n" 
384      usage = usage + "                         [default = %g]\n" % psize.getConstant("TFAC_ALPHA") 
385      usage = usage + "  --TFAC_XEON=<value>  : Number of sec/unknown for setup/solve on 500\n" 
386      usage = usage + "                         MHz PIII Xeon CPU -- VERY ROUGH ESTIMATE\n" 
387      usage = usage + "                         [default = %g]\n" % psize.getConstant("TFAC_XEON") 
388      usage = usage + "  --TFAC_SPARC=<value> : Number of sec/unknown for setup/solve on 400\n" 
389      usage = usage + "                         MHz UltraSPARC II CPU -- VERY ROUGH ESTIMATE\n" 
390      usage = usage + "                         [default = %g]\n" % psize.getConstant("TFAC_SPARC") 
391   
392       
393      sys.stderr.write(usage) 
394      sys.exit(2) 
395   
396 -def main(): 
397      import getopt 
398      filename = "" 
399      shortOptList = "" 
400      longOptList = ["help", "CFAC=", "FADD=", "SPACE=", "GMEMFAC=", "GMEMCEIL=", "OFAC=", "REDFAC=", "TFAC_ALPHA=", "TFAC_XEON=", "TFAC_ALPHA="] 
401      try: 
402          opts, args = getopt.getopt(sys.argv[1:], shortOptList, longOptList) 
403      except getopt.GetoptError, details: 
404          sys.stderr.write("Option error (%s)!\n" % details) 
405          usage() 
406      if len(args) != 1:  
407          sys.stderr.write("Invalid argument list!\n") 
408          usage() 
409      else: 
410          filename = args[0] 
411   
412      psize = Psize()     
413   
414      for o, a in opts: 
415          if o == "--help": 
416              usage() 
417          if o == "--CFAC": 
418              psize.setConstant("CFAC", float(a)) 
419          if o == "--FADD": 
420              psize.setConstant("FADD", int(a)) 
421          if o == "--SPACE": 
422              psize.setConstant("SPACE", float(a)) 
423          if o == "--GMEMFAC": 
424              psize.setConstant("GMEMFAC", int(a)) 
425          if o == "--GMEMCEIL": 
426              psize.setConstant("GMEMCEIL",  int(a)) 
427          if o == "--OFAC": 
428              psize.setConstant("OFAC", float(a)) 
429          if o == "--REDFAC": 
430              psize.setConstant("REDFAC", float(a)) 
431          if o == "--TFAC_ALPHA": 
432              psize.setConstant("TFAC_ALPHA", float(a)) 
433          if o == "--TFAC_XEON": 
434              psize.setConstant("TFAC_XEON",  float(a)) 
435          if o == "--TFAC_SPARC":     
436              psize.setConstant("TFAC_SPARC",  float(a)) 
437   
438      psize.runPsize(filename) 
439      sys.stdout.write("Default constants used (./psize.py --help for more information):\n") 
440      sys.stdout.write("%s\n" % psize.constants) 
441      sys.stdout.write(psize.printResults()) 
442   
443       
444  if __name__ == "__main__": main() 
445