Source Code for Module PyAutoDock.electrostatics

  1  ############################################################################ 
  2  # 
  3  # Authors: William Lindstrom, Ruth Huey 
  4  # 
  5  # Copyright: A. Olson TSRI 2004 
  6  # 
  7  ############################################################################# 
  8   
  9  # 
 10  # $Id: electrostatics.py,v 1.3 2004/07/23 00:54:46 lindy Exp $ 
 11  # 
 12   
 13  import math 
 14  from scorer import DistDepPairwiseScorer 
 15   
 16   
 17 -class ElectrostaticsRefImpl(DistDepPairwiseScorer): 
 18      """reference implementation of the Electrostatics class 
 19   
 20      Notes 
 21      ====== 
 22      * For now, the MolecularSystem has the responsibility of 
 23        supplying charges. This could change in the future if the 
 24        chargeCalculators migrate away from MolKit. This requires 
 25        a reduced dependency on the Atoms. 
 26   
 27      * The electrostatics calculation depends on relative distance 
 28        not absolute coords. Hopefully, we won't need coords in this 
 29        class. 
 30   
 31      * A GridMap might be represencted as a 'subset' in this class 
 32        by defining psuedo-atoms at each grid point. 
 33   
 34      """ 
 35 -    def __init__(self, ms=None): 
 36          DistDepPairwiseScorer.__init__(self) 
 37          # add the required attributes of this subclass to the dict 
 38          self.required_attr_dictA.setdefault('charge', False) 
 39          self.required_attr_dictB.setdefault('charge', False) 
 40          if ms is not None: 
 41              self.set_molecular_system(ms) 
 42          self.use_non_bond_cutoff = False 
 43   
 44   
 45 -    def get_dddp(self, distance): 
 46          """return distance dependent dielectric permittivity 
 47   
 48          weight to be remove to WeightedMultiTermDistDepPairwiseScorer 
 49           
 50          Reference 
 51          ========== 
 52          Mehler & Solmajer (1991) Protein Engineering vol.4 no.8, pp. 903-910. 
 53   
 54          Basically, this describes a sigmoidal function from a low 
 55          dielectric (like organic solvent) at small distance (<5Ang) 
 56          asymptoticly to dielectric of water at large distance (30Ang). 
 57           
 58          epsilon(r) = A + (B/[1+ k*exp(-lambda*B*r)])     (equation 6) 
 59   
 60          where, A, lambda, and k are parameters supplied by the paper, 
 61                 B = epsilon0 - A (so B is also a parameter) 
 62                 epsilon0 is the dielectric constant of water at 25C (78.4) 
 63   
 64          Two sets of parameters are given in the paper: 
 65          {'A' : -8.5525, 'k' : 7.7839, 'lambda_' : 0.003627} [Conway] 
 66          {'A' : -20.929, 'k' : 3.4781, 'lambda_' : 0.001787} [Mehler&Eichele] 
 67          """ 
 68           
 69          epsilon0 = 78.4     # dielectric constant of water at 25C 
 70          lambda_ = 0.003627  # supplied parameter 
 71          A = -8.5525         # supplied parameter 
 72          k = 7.7839          # supplied parameter 
 73          B = epsilon0 - A 
 74   
 75          # epsilon is a unitless scaling factor  
 76          epsilon = A + B / (1.0 + k*math.e**(-lambda_*B*distance)) 
 77          return epsilon 
 78       
 79   
 80   
 81 -    def _f(self, at_a, at_b, dist): 
 82          """Return electrostatic potential in kcal/mole. 
 83   
 84          Here's how the 332. unit conversion factor is calculated: 
 85           
 86          q = 1.60217733E-19  # (C) charge on electron 
 87          eps0 = 8.854E-12    # (C^2/J*m) vacuum permittivity 
 88          J_per_cal = 4.18400 # Joules per calorie 
 89          avo = 6.02214199E23 # Avogadro's number 
 90   
 91          factor = (q*avo*q*1E10 )/(4.0*math.pi*eps0*J_per_cal*1000) 
 92   
 93          """ 
 94          q1 = at_a.charge; 
 95          q2 = at_b.charge; 
 96          epsilon = self.get_dddp(dist) 
 97          factor = 332. # 4*pi*esp0*units conversions 
 98          if dist != 0.0: 
 99              return (factor*q1*q2) / (epsilon*dist) # kcal/mole 
100          else: return 0. 
101   
102   
103  Electrostatics = ElectrostaticsRefImpl 
104  # ElectrostaticsRefImpl 
105   
106  if __name__ == '__main__': 
107      pass 
108      # test_scorer.run_test() 
109