Source Code for Module sff.amber

  1  # 
  2  # Copyright_notice 
  3  # 
  4   
  5  import _prmlib as prmlib 
  6  import re 
  7  import os 
  8  import Numeric 
  9  from types import StringType 
 10   
 11   
 12  realType = Numeric.Float 
 13  intType = Numeric.Int 
 14  pyArrayInt = prmlib.PyArray_INT 
 15  pyArrayDouble = prmlib.PyArray_DOUBLE 
 16   
 17  getpat = re.compile( 'parmstruct_(\w+)_get') 
 18  parmattrs = {} 
 19   
 20   
 21  for x in dir( prmlib): 
 22      match = getpat.match( x) 
 23      if match: 
 24          parmattrs[ match.group( 1) ] = None 
 25               
 26  parmbuffers = { 
 27      'AtomNames': (StringType, lambda x: x.Natom * 4 + 81, None), 
 28      'Charges': (realType, lambda x: x.Natom, pyArrayDouble ), 
 29      'Masses': (realType, lambda x: x.Natom, pyArrayDouble), 
 30      'Iac': (intType, lambda x: x.Natom, pyArrayInt), 
 31      'Iblo': (intType, lambda x: x.Natom, pyArrayInt), 
 32      'Cno': (intType, lambda x: x.Ntype2d, pyArrayInt), 
 33      'ResNames': (StringType, lambda x: x.Nres * 4 + 81, None), 
 34      'Ipres': (intType, lambda x: x.Nres + 1, pyArrayInt), 
 35      'Rk': (realType, lambda x: x.Numbnd, pyArrayDouble), 
 36      'Req': (realType, lambda x: x.Numbnd, pyArrayDouble), 
 37      'Tk': (realType, lambda x: x.Numang, pyArrayDouble), 
 38      'Teq': (realType, lambda x: x.Numang, pyArrayDouble), 
 39      'Pk': (realType, lambda x: x.Nptra, pyArrayDouble), 
 40      'Pn': (realType, lambda x: x.Nptra, pyArrayDouble), 
 41      'Phase': (realType, lambda x: x.Nptra, pyArrayDouble), 
 42      'Solty': (realType, lambda x: x.Natyp, pyArrayDouble), 
 43      'Cn1': (realType, lambda x: x.Nttyp, pyArrayDouble), 
 44      'Cn2': (realType, lambda x: x.Nttyp, pyArrayDouble), 
 45      'Boundary': (intType, lambda x: x.Nspm, pyArrayInt), 
 46      'BondHAt1': (intType, lambda x: x.Nbonh, pyArrayInt), 
 47      'BondHAt2': (intType, lambda x: x.Nbonh, pyArrayInt), 
 48      'BondHNum': (intType, lambda x: x.Nbonh, pyArrayInt), 
 49      'BondAt1': (intType, lambda x: x.Nbona, pyArrayInt), 
 50      'BondAt2': (intType, lambda x: x.Nbona, pyArrayInt), 
 51      'BondNum': (intType, lambda x: x.Nbona, pyArrayInt), 
 52      'AngleHAt1': (intType, lambda x: x.Ntheth, pyArrayInt), 
 53      'AngleHAt2': (intType, lambda x: x.Ntheth, pyArrayInt), 
 54      'AngleHAt3': (intType, lambda x: x.Ntheth, pyArrayInt), 
 55      'AngleHNum': (intType, lambda x: x.Ntheth, pyArrayInt), 
 56      'AngleAt1': (intType, lambda x: x.Ntheta, pyArrayInt), 
 57      'AngleAt2': (intType, lambda x: x.Ntheta, pyArrayInt), 
 58      'AngleAt3': (intType, lambda x: x.Ntheta, pyArrayInt), 
 59      'AngleNum': (intType, lambda x: x.Ntheta, pyArrayInt), 
 60      'DihHAt1': (intType, lambda x: x.Nphih, pyArrayInt), 
 61      'DihHAt2': (intType, lambda x: x.Nphih, pyArrayInt), 
 62      'DihHAt3': (intType, lambda x: x.Nphih, pyArrayInt), 
 63      'DihHAt4': (intType, lambda x: x.Nphih, pyArrayInt), 
 64      'DihHNum': (intType, lambda x: x.Nphih, pyArrayInt), 
 65      'DihAt1': (intType, lambda x: x.Nphia, pyArrayInt), 
 66      'DihAt2': (intType, lambda x: x.Nphia, pyArrayInt), 
 67      'DihAt3': (intType, lambda x: x.Nphia, pyArrayInt), 
 68      'DihAt4': (intType, lambda x: x.Nphia, pyArrayInt), 
 69      'DihNum': (intType, lambda x: x.Nphia, pyArrayInt), 
 70      'ExclAt': (intType, lambda x: x.Nnb, pyArrayInt), 
 71      'HB12': (realType, lambda x: x.Nphb, pyArrayDouble), 
 72      'HB10': (realType, lambda x: x.Nphb, pyArrayDouble), 
 73      'Box': (realType, lambda x: 3, pyArrayDouble), 
 74      'AtomSym': (StringType, lambda x: x.Natom * 4 + 81, None), 
 75      'AtomTree': (StringType, lambda x: x.Natom * 4 + 81, None), 
 76      'TreeJoin': (intType, lambda x: x.Natom, pyArrayInt), 
 77      'AtomRes': (intType, lambda x: x.Natom, pyArrayInt), 
 78      'N14pairs': (intType, lambda x: x.Natom, pyArrayInt), 
 79      'N14pairlist': (intType, lambda x: 10*x.Natom, pyArrayInt), 
 80     } 
 81       
 82 -def checkbuffersize( parmobj, attr, value): 
 83      attrlen = parmbuffers[ attr][ 1]( parmobj) 
 84      if attr in ['AtomNames', 'AtomSym', 'AtomTree']: 
 85          attrlen = parmobj.Natom * 4 
 86      elif attr == 'ResNames': 
 87          attrlen = parmobj.Nres * 4 
 88      elif attr == 'Ipres': 
 89          attrlen = parmobj.Nres 
 90      elif attr == 'N14pairlist': 
 91          from operator import add 
 92          sum = reduce( add, parmobj.N14pairs ) 
 93          attrlen = sum 
 94          if sum!=len(value): 
 95              print 'WARNING: N14pairlist length' 
 96              attrlen = len(value) 
 97               
 98      if len( value) < attrlen: 
 99          raise ValueError( attr, attrlen, len( value)) 
100   
101 -class AmberParm: 
102 -    def __init__( self, name, parmdict=None): 
103          """ 
104          name - string 
105          parmdict - map 
106          """ 
107          import types 
108   
109          self.name  = name 
110          if parmdict: 
111              parmptr = self._parmptr_ = prmlib.parmcalloc() 
112              for name in parmattrs.keys(): 
113                  value = parmdict[ name] 
114                  try: 
115                      bufdesc = parmbuffers[ name] 
116                  except KeyError: 
117                      pass 
118                  else: 
119                      if bufdesc[ 0] != StringType\ 
120                         and not isinstance( value, StringType): 
121                          value = Numeric.array( value).astype(bufdesc[ 0]) 
122                  self.__dict__[ name] = value 
123                  if name == 'Box': 
124                      self.Box[:] = value 
125                  else: 
126                      getattr( prmlib, 'parmstruct_%s_set' % name)( parmptr, value) 
127   
128          else: 
129              assert os.path.exists( name ) 
130              self._parmptr_ = parmptr = prmlib.readparm( name) 
131              for attr in filter( lambda x: not parmbuffers.has_key( x), 
132                                  parmattrs.keys()): 
133                  value = getattr( prmlib, 'parmstruct_%s_get' % attr)( parmptr) 
134                  self.__dict__[ attr] = value 
135                  
136              for attr in filter( lambda x: parmbuffers.has_key( x), 
137                                  parmattrs.keys()): 
138                  # these _get() functions  must not be called from anywhere else 
139                  #print "attr:", attr, 
140                  value = getattr( prmlib, 'parmstruct_%s_get' % attr)( parmptr) 
141                  #print "value: ", value 
142                  if value is None: 
143                      value = () 
144                  else: 
145                      bufdesc = parmbuffers[ attr] 
146                      if bufdesc[ 0] != StringType: 
147                          value = prmlib.createNumArr(value,  bufdesc[ 1]( self), bufdesc[2]) 
148   
149                  self.__dict__[ attr] = value 
150          if __debug__: 
151              for attr in parmbuffers.keys(): 
152                  val = getattr(self, attr) 
153                  if isinstance(val, Numeric.ArrayType) or isinstance(val, StringType): 
154                      checkbuffersize(self, attr, val) 
155               
156 -    def __setattr__( self, name, value): 
157          if parmattrs.has_key( name): 
158              raise AttributeError( 'constant parm attribute') 
159          self.__dict__[ name] = value 
160   
161 -    def __del__( self): 
162          prmlib.parmfree( self._parmptr_) 
163          delattr( self, '_parmptr_') 
164           
165   
166 -    def asDict(self): 
167          # return the content of the parm structure as a python dict 
168          d = {} 
169          for k in self.__dict__.keys(): 
170              if k[0]=='_': continue 
171              if parmbuffers.has_key(k): 
172                  value = list(getattr(self, k)) 
173              else: 
174                  value = getattr(self, k) 
175              d[k] = value 
176          return d 
177   
178   
179  import threading 
180  amberlock = threading.RLock() 
181   
182  import struct, Numeric 
183   
184 -class BinTrajectory: 
185       
186      ## WARNING nothing is done about byte order 
187 -    def __init__(self, filename): 
188          import os 
189          assert os.path.isfile(filename) 
190          self.filename = filename 
191          self.typecode = 'f' 
192          if prmlib.UseDouble: 
193              self.typecode = 'd' 
194          self.coordSize = struct.calcsize(self.typecode) 
195          self.intSize = struct.calcsize('i') 
196          self.fileHandle = None 
197           
198   
199 -    def getNumberOfAtoms(self, filename): 
200          f = open(filename, 'rb') 
201          lenstr = f.read(struct.calcsize('i')) 
202          f.close() 
203          return struct.unpack('i', lenstr)[0] 
204       
205   
206 -    def closeFile(self): 
207          self.fileHandle.close() 
208          self.fileHandle = None 
209   
210           
211 -    def getNextConFormation(self): 
212          # open file if necessary 
213          if self.fileHandle is None: 
214              self.fileHandle = open(self.filename) 
215   
216          # read the number of atoms as an integer 
217          lenstr = self.fileHandle.read(self.intSize) 
218          if len(lenstr) < self.intSize: #EOF reached 
219              self.closeFile() 
220              return None 
221          nba = struct.unpack('i', lenstr)[0] 
222          size = 3 * nba * self.coordSize 
223   
224          # read the coordinates for nba atoms 
225          crdstr = self.fileHandle.read( size ) 
226          if len(crdstr) != size: #EOF reached 
227              self.closeFile() 
228              return None 
229          c = Numeric.array( struct.unpack( '%dd'%3*nba, crdstr), 
230                             self.typecode ) 
231          c.shape = (-1, 3) 
232          return c 
233   
234           
235   
236 -class Amber94: 
237   
238      from MolKit import parm94_dat 
239   
240 -    def __init__(self, atoms, parmtop=None, prmfile=None, dataDict={}): 
241   
242          from MolKit.amberPrmTop import Parm 
243   
244          self.atoms = atoms 
245          self.parmtop = parmtop 
246          if prmfile: 
247              self.oprm = AmberParm( prmfile ) 
248          else: 
249              if self.parmtop is None: 
250                  # create parmtop info 
251                  if not len(dataDict): 
252                      self.parmtop = Parm() 
253                  else: 
254                      #dataDict is a dictionary with possible keys: 
255                      #allDictList, ntDictList, ctDictList 
256                      #whose values are lists of python files such as 
257                      #found in MolKit/data...which end in dat.py 
258                      #dataDict['allDictList']=[all_amino94_dat] 
259                      #if len(list)>1, the first is updated by the rest 
260                      self.parmtop = apply(Parm, (), dataDict) 
261                  self.parmtop.processAtoms(atoms, self.parm94_dat) 
262                  #this read is broken 
263                  #self.parmtop.loadFromFile(prmfile) 
264              else: 
265                  assert isinstance(parmtop, Parm) 
266   
267              # create the C-data structure 
268              self.oprm = AmberParm( 'test', self.parmtop.prmDict ) 
269   
270          from operator import add 
271          coords = self.atoms.coords[:] 
272          lcoords = reduce( add, coords) 
273          self.coords = Numeric.array( lcoords).astype(realType ) 
274           
275          # create Numeric array for frozen 
276          self.frozen = Numeric.zeros( self.oprm.Natom).astype(intType) 
277   
278          # create Numeric array for constrained 
279          self.constrained = Numeric.zeros( self.oprm.Natom). astype(intType) 
280   
281          # create Numeric array for anchor 
282          self.anchor = Numeric.zeros( 3*self.oprm.Natom).astype(realType)  
283   
284          # create Numeric array for minv (md) 
285          self.minv = Numeric.zeros( 3*self.oprm.Natom).astype(realType ) 
286           
287          # create Numeric array for v (md) 
288          self.mdv = Numeric.zeros( 3*self.oprm.Natom).astype(realType) 
289           
290          # create Numeric array for f (md) 
291          self.mdf = Numeric.zeros( 3*self.oprm.Natom).astype( realType ) 
292   
293          # is the number of variables 
294          self.nbVar = Numeric.array([3*self.oprm.Natom]).astype(intType) 
295   
296          # will contain the value of the objective function at the end 
297          self.objFunc = Numeric.zeros( 1).astype(realType ) 
298   
299          # return when sum of squares of gradient is less than dgrad 
300          drms = 0.1 
301          self.dgrad = Numeric.array([drms*3*self.oprm.Natom]).astype(realType) 
302   
303          # expected decrease in the function on the first iteration 
304          self.dfpred = Numeric.array( [10.0,]).astype( realType ) 
305   
306          #  
307          self.maxIter = Numeric.array([500,]).astype(intType) 
308   
309          # 
310          self.energies = Numeric.zeros(20).astype(realType ) 
311   
312          # filename used to save trajectory 
313          self.filename = None 
314          self.sff_opts = prmlib.init_sff_options() 
315   
316 -    def setMinimizeOptions(self, **kw): 
317          # WARNING when cut it set mme_init needs to be called to allocate a 
318          # list of non-bonded paires of the proper size 
319          for k,v in kw.items(): 
320              assert k in ['cut', 'nsnb', 'ntpr', 'scnb', 'scee', 
321                           'mme_init_first', 'dield', 'verbosemm', 
322                           'wcons'] 
323              #prmlib.mm_options(k, v) 
324              #setattr(prmlib.cvar, k, v) 
325              prmlib.mm_options(k, v, self.sff_opts) 
326           
327   
328 -    def setMdOptions(self, **kw): 
329          #nb: for the moment set verbosemm for verbosemd 
330          for k,v in kw.items(): 
331              assert k in ['t', 'dt', 'tautp', 'temp0', 'boltz2', 'verbosemd', 
332                            'ntwx','vlimit', 'ntpr_md', 'zerov', 'tempi', 'idum' ] 
333              #prmlib.md_options(k, v) 
334              #setattr(prmlib.cvar, k, v) 
335              prmlib.md_options(k, v, self.sff_opts) 
336   
337   
338 -    def setCallback(self, func, frequency): 
339          assert callable(func) 
340          prmlib.set_callback(func, frequency, 0) 
341   
342           
343 -    def freezeAtoms(self, atomIndices): 
344          assert len(atomIndices)==len(self.atoms), 'atomIndices wrong length' 
345          self.frozen = Numeric.array(atomIndices).astype(intType) 
346   
347   
348 -    def constrainAtoms(self, atomIndices, anchor): 
349          atlen = len(self.atoms) 
350          assert len(atomIndices)==atlen, 'atomIndices wrong length' 
351          #this is not right: 
352          #constNum = Numeric.add.reduce(atomIndices) 
353          #anchors have garbage for non-constrained atoms 
354          assert len(anchor)==atlen*3, 'anchor wrong length' 
355          self.constrained = Numeric.array(atomIndices).astype(intType) 
356          self.anchor = Numeric.array(anchor).astype(realType) 
357       
358 -    def minimize(self, drms=None, maxIter=None, dfpred=None): 
359   
360          if drms is not None: self.dgrad[0] = drms*3*self.oprm.Natom 
361          if maxIter is not None: self.maxIter[0] = maxIter 
362          if dfpred is not None: self.dfpred[0] = dfpred 
363               
364          prmlib.mme_init(self.frozen, self.constrained, self.anchor, 
365                          None, self.oprm._parmptr_, self.sff_opts) 
366   
367          amberlock.acquire() 
368          result_code = -6 
369          try: 
370              # add new thread here 
371              result_code = prmlib.conjgrad( self.coords, self.nbVar, self.objFunc, 
372                     prmlib.mme_fun, self.dgrad, self.dfpred, self.maxIter, 
373                     self.oprm._parmptr_, self.energies, self.sff_opts ) 
374          finally: 
375              amberlock.release() 
376   
377          return result_code 
378   
379       
380 -    def md(self, maxStep, filename=None): 
381   
382          self.filename = filename 
383          if filename is not None: 
384              f = open(filename, 'w') 
385          else: 
386              f = None 
387   
388          prmlib.mme_init(self.frozen, self.constrained, self.anchor, 
389                          f, self.oprm._parmptr_, self.sff_opts) 
390   
391          amberlock.acquire() 
392          result_code = -6 
393          try: 
394              # add new thread here 
395              result_code = prmlib.md( 3*self.oprm.Natom, maxStep, self.coords, 
396                                       self.minv, self.mdv, self.mdf, 
397                                       prmlib.mme_fun, self.energies, 
398                                       self.oprm._parmptr_ , self.sff_opts) 
399          finally: 
400              amberlock.release() 
401   
402          if filename is not None: 
403              f.close() 
404   
405          return result_code 
406