PyAutoDock.desolvation

1 ############################################################################ 2 # 3 # Authors: William Lindstrom, Ruth Huey 4 # 5 # Copyright: A. Olson TSRI 2004 6 # 7 ############################################################################# 8 9 # 10 # $Id: desolvation.py,v 1.11 2007/06/14 23:30:35 gillet Exp $ 11 # 12 13 14 import Numeric, math 15 from scorer import DistDepPairwiseScorer 16 17 18

19 -class DesolvationRefImpl(DistDepPairwiseScorer):

20 - def __init__(self, ms=None):

21 DistDepPairwiseScorer.__init__(self) 22 # add the required attributes of this subclass to the dict 23 self.required_attr_dictA.setdefault('autodock_element', False) 24 self.required_attr_dictA.setdefault('AtVol', False) #only rec needs AtVol 25 self.required_attr_dictB.setdefault('autodock_element', False) 26 if ms is not None: 27 self.set_molecular_system(ms) 28 29 # set up or get desolvation parameters for ligand 30 # access these dictionaries with defaults for missing keys 31 self.sol_par = { 32 'C': 0.004, 33 'A': 0.0006} # units of kcals/ Ang**2 34 35 # In Autogrid3, a constant is added to each map according 36 # to the sol_const dictionary below. In Autogrid3 the constant 37 # is applied to the map independent of the empirical weight. 38 # Here, we apply the constant before the weight and must divide 39 # the constant by the weight in order to recover the constant when 40 # the weight is applied. 41 # PROBLEM: our term now knows about its Autogrid3 weight. 42 dsolv_weight = 0.1711 # Autogrid3 weight 43 self.sol_const = { 44 'O': 0.236/dsolv_weight, 45 'H': 0.118/dsolv_weight, 46 }

47 48

49 - def get_ddsol(self, distance):

50 """return distance dependent solvation parameter 51 52 WHERE DID THIS COME FROM?? 53 =========================== 54 55 Stouten @@ Document this @@ 56 """ 57 sigma = 3.6 58 minus_inv_two_sigma_sqd = -1.0 / (2.0*(sigma*sigma)) 59 exponent = minus_inv_two_sigma_sqd*(distance*distance) 60 ddsol = math.e**exponent 61 return ddsol

62 63

64 - def _f(self, at_a, at_b, dist):

65 """Return desolvation energy in kcal/mole. 66 67 Atom Volumes are in units of kcal/mol*Ang**3 68 NB: pdbqs volumes are cal/mol*Ang**3 and must be converted!! 69 This is done by MolecularSystem.py for now. 70 71 In Autodock3.0.5, this term is valid only for C atoms 72 in the ligand and is used to descriminate aromatic carbons 73 from aliphatic carbons 74 """ 75 elm_a = at_a.autodock_element 76 #NB we hope we have checked for AtVol by this point 77 # in the MolecularSystem._check_interface 78 vol_a = at_a.AtVol 79 elm_b = at_b.autodock_element 80 energy = 0.0 81 sp = self.sol_par.get(elm_b, 0.0) # dict with default 82 if sp != 0.0: 83 energy = - 1.0 * sp * vol_a * self.get_ddsol(dist) 84 return energy # kcal/mol

86 - def post_process(self):

87 88 # ligand (by convention) 89 atoms_bx = self.ms.configuration[1] 90 atoms_b = self.ms.get_entities(atoms_bx) 91 92 # get the number of receptor atoms 93 atoms_ax = self.ms.configuration[0] 94 atoms_a = self.ms.get_entities(0) 95 num_R_atoms = len(atoms_a) 96 # Oxygen 97 # 1.3793 = 0.236/0.1711 (so these are unweighted values here!) 98 smeared_Oxygen_sol_const = 1.3793/ num_R_atoms 99 100 # Hydrogen 101 # 0.6897 = 0.110/0.1711 (so these are unweighted values here!) 102 smeared_Hydrogen_sol_const = 0.6897/ num_R_atoms 103 104 # these desolvation constants come from the gpf "constant" parameter 105 # where they are weighted by the desolvation weight (0.1711) 106 107 # get the score array from the derived class 108 # self.score_array = self.get_score_array() 109 110 for row in xrange(len(atoms_a)): 111 # we need to call get_entities to reset the iterator after firts loop 112 for i,ats in enumerate(self.ms.get_entities(atoms_bx)): 113 element = ats.element 114 if element == 'H': 115 self.array[row][i] += smeared_Hydrogen_sol_const 116 elif element =='O': 117 self.array[row][i] += smeared_Oxygen_sol_const

118 119 Desolvation = DesolvationRefImpl 120 # Desolvation 121 122 123 124

125 -class NewDesolvationRefImpl(DistDepPairwiseScorer):

126 127 #Vol = 4/3*pi*r**3 where r is Rij[XX]/2. 128 Vols = {} 129 Vols["C"] = Vols["C"] = Vols["A"] = 33.5103 130 Vols["N"] = Vols["NA"] = Vols["Na"] = Vols["n"] = 22.4493 131 Vols["O"] = Vols["OA"] = Vols["Oa"] = 17.1573 132 Vols["P"] = 38.7924 133 Vols["S"] = Vols["SA"] = Vols["Sa"] = 33.5103 134 Vols["H"] = Vols["HD"] = Vols["HS"] = Vols["Hd"] = Vols["Hs"] =0.00 135 Vols["F"] = 15.448 136 Vols["f"] = Vols["FE"] = Vols["Fe"] = 1.84 #Iron oldvalue==1.1503 137 Vols["I"] = 55.0585 138 Vols["M"] = Vols["MG"] = Vols["Mg"] = 1.56 # Magnesium oldvalue=1.1503 139 Vols["MN"]= 2.14 # Manganese 140 Vols["Z"] = Vols["ZN"] = Vols["Zn"] = 1.70 # Zinc oldvalue = 1.6974 141 Vols["L"] = Vols["CA"] = Vols["Ca"] = 2.77 # Calcium oldvalue = 4.0644 142 Vols["c"] = Vols["CL"] = Vols["Cl"] = 35.8235 # Chlorine 143 Vols["b"] = Vols["BR"] = Vols["Br"] = 42.5661 # Bromine 144 145 # set up or get desolvation parameters 146 # access these dictionaries with defaults for missing keys 147 Solpars = { #oldvalue * 0.10188 148 'C': -0.00143, #oldvalue = -0.00143 149 'A': -0.00052, #oldvalue = -0.00052 150 'N': -0.00162, #oldvalue = -0.00162 151 'NA':-0.00162, #oldvalue = -0.00162 152 'O': -0.00251, #oldvalue = -0.00251 153 'OA':-0.00251, #oldvalue = -0.00251 154 'H': 0.00051, #oldvalue = 0.00051 155 'HD': 0.00051, #oldvalue = 0.00051 156 'S': -0.00214, #oldvalue = -0.00214 157 'SA':-0.00214} #oldvalue = -0.00214 158 # units of kcals/ Ang**2 159 160

161 - def __init__(self, ms=None, solpar_q = 0.01097): #oldvalue = 0.01097

162 DistDepPairwiseScorer.__init__(self) 163 # add the required attributes of this subclass to the dict 164 self.required_attr_dictA = {} 165 self.required_attr_dictA.setdefault('autodock_element', False) 166 self.required_attr_dictB.setdefault('autodock_element', False) 167 if ms is not None: 168 self.set_molecular_system(ms) 169 170 #set up constant 171 self.solpar_q = solpar_q

172 173 174

175 - def get_ddsol(self, distance):

176 """return distance dependent solvation parameter 177 178 WHERE DID THIS COME FROM?? 179 =========================== 180 181 Stouten @@ Document this @@ 182 """ 183 sigma = 3.6 184 minus_inv_two_sigma_sqd = -1.0 / (2.0*(sigma*sigma)) 185 exponent = minus_inv_two_sigma_sqd*(distance*distance) 186 ddsol = math.e**exponent 187 return ddsol

188 189

190 - def _f(self, at_a, at_b, dist):

191 """Return desolvation energy in kcal/mole. 192 193 Atom Volumes are in units of kcal/mol*Ang**3 194 NB: pdbqs volumes are cal/mol*Ang**3 and must be converted!! 195 This is done by MolecularSystem.py for now. 196 197 In Autodock3.0.5, this term is valid only for C atoms 198 in the ligand and is used to discriminate aromatic carbons 199 from aliphatic carbons 200 """ 201 elm_a = at_a.autodock_element 202 charge_a = at_a.charge 203 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 204 vol_a = self.Vols.get(elm_a, 33.5103) 205 #vol_a = at_a.AtVol 206 #-.0.00110*0.10188, the weight from the recalibration 207 sp_a = self.Solpars.get(elm_a, -0.00110) # dict with default 208 #sp_a = self.Solpars.get(elm_a, 0.0) # dict with default 209 #NB we hope we have checked for AtVol by this point 210 # in the MolecularSystem._check_interface 211 elm_b = at_b.autodock_element 212 charge_b = at_b.charge 213 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 214 vol_b = self.Vols.get(elm_b, 33.5103) 215 #vol_b = at_b.AtVol 216 #-.0.00110*0.10188, the weight from the recalibration 217 sp_b = self.Solpars.get(elm_b, -0.00110) # dict with default 218 219 lig_energy = 0.0 220 rec_energy = 0.0 221 222 sol_func = self.get_ddsol(dist) 223 #if sp != 0.0: 224 # energy = - sp * vol_a * self.get_ddsol(dist) 225 rec_at_energy = (sp_a + self.solpar_q*abs(charge_a))*vol_b * sol_func 226 lig_at_energy = (sp_b + self.solpar_q*abs(charge_b))*vol_a * sol_func 227 return rec_at_energy + lig_at_energy # kcal/mol

228 229 230 NewDesolvation = NewDesolvationRefImpl 231 #end NewDesolvation 232 233 234

235 -class NewDesolvationLigOnlyRefImpl(NewDesolvation):

236 237 #def __init__(self, ms=None, solpar_q=0.01097): 238 # NewDesolvation.__init__(self, ms=ms, solpar_q=solpar_q) 239 240

241 - def _f(self, at_a, at_b, dist):

242 """Return desolvation energy in kcal/mole. 243 244 Atom Volumes are in units of kcal/mol*Ang**3 245 NB: pdbqs volumes are cal/mol*Ang**3 and must be converted!! 246 This is done by MolecularSystem.py for now. 247 248 In Autodock3.0.5, this term is valid only for C atoms 249 in the ligand and is used to discriminate aromatic carbons 250 from aliphatic carbons 251 """ 252 elm_a = at_a.autodock_element 253 charge_a = at_a.charge 254 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 255 vol_a = self.Vols.get(elm_a, 33.5103) 256 #vol_a = at_a.AtVol 257 #-.0.00110*0.10188, the weight from the recalibration 258 sp_a = self.Solpars.get(elm_a, -0.00110) # dict with default 259 #sp_a = self.Solpars.get(elm_a, 0.0) # dict with default 260 #NB we hope we have checked for AtVol by this point 261 # in the MolecularSystem._check_interface 262 elm_b = at_b.autodock_element 263 charge_b = at_b.charge 264 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 265 vol_b = self.Vols.get(elm_b, 33.5103) 266 #vol_b = at_b.AtVol 267 #-.0.00110*0.10188, the weight from the recalibration 268 sp_b = self.Solpars.get(elm_b, -0.00110) # dict with default 269 270 lig_energy = 0.0 271 rec_energy = 0.0 272 273 sol_func = self.get_ddsol(dist) 274 #if sp != 0.0: 275 # energy = - sp * vol_a * self.get_ddsol(dist) 276 #rec_at_energy = (sp_a + self.solpar_q*abs(charge_a))*vol_b * sol_func 277 lig_at_energy = (sp_b + self.solpar_q*abs(charge_b))*vol_a * sol_func 278 #return rec_at_energy + lig_at_energy # kcal/mol 279 return lig_at_energy

280 281 NewDesolvationLigOnly = NewDesolvationLigOnlyRefImpl 282 #end NewDesolvationLigOnly 283 284 285 #NewDesolvationAtomMap

286 -class NewDesolvationAtomMapRefImpl(NewDesolvation):

287 288 #def __init__(self, ms=None, solpar_q = 0.01097): 289 # NewDesolvation.__init__(self, ms=ms, solpar_q=solpar_q) 290 291

292 - def _f(self, at_a, at_b, dist):

293 """Return desolvation energy in kcal/mole. 294 295 Atom Volumes are in units of kcal/mol*Ang**3 296 NB: pdbqs volumes are cal/mol*Ang**3 and must be converted!! 297 This is done by MolecularSystem.py for now. 298 299 In Autodock3.0.5, this term is valid only for C atoms 300 in the ligand and is used to discriminate aromatic carbons 301 from aliphatic carbons 302 """ 303 elm_a = at_a.autodock_element 304 charge_a = at_a.charge 305 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 306 vol_a = self.Vols.get(elm_a, 33.5103) 307 sp_a = self.Solpars.get(elm_a, -0.00110) # dict with default 308 #NB we hope we have checked for AtVol by this point 309 # in the MolecularSystem._check_interface 310 elm_b = at_b.autodock_element 311 charge_b = at_b.charge 312 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 313 vol_b = self.Vols.get(elm_b, 33.5103) 314 sp_b = self.Solpars.get(elm_b, -0.00110) # dict with default 315 316 lig_energy = 0.0 317 rec_energy = 0.0 318 319 sol_func = self.get_ddsol(dist) 320 # 321 rec_at_energy = (sp_a + self.solpar_q*abs(charge_a))*vol_b * sol_func 322 #NOTE THIS IS TRUNCATED: 323 #lig_at_energy = (sp_b + self.solpar_q*abs(charge_b))*vol_a * sol_func 324 #this portion goes into the atom map 325 lig_at_energy = (sp_b)*vol_a * sol_func 326 # 327 return rec_at_energy + lig_at_energy # kcal/mol

328 329 NewDesolvationAtomMap = NewDesolvationAtomMapRefImpl 330 331 332 #NewDesolvationDesolvMap

333 -class NewDesolvationDesolvMapRefImpl(NewDesolvation):

334 335 #def __init__(self, ms=None, solpar_q = 0.01097): 336 # NewDesolvation.__init__(self, ms=ms, solpar_q=solpar_q) 337 338

339 - def _f(self, at_a, at_b, dist):

340 """Return desolvation energy in kcal/mole. 341 342 Atom Volumes are in units of kcal/mol*Ang**3 343 NB: pdbqs volumes are cal/mol*Ang**3 and must be converted!! 344 This is done by MolecularSystem.py for now. 345 346 In Autodock3.0.5, this term is valid only for C atoms 347 in the ligand and is used to discriminate aromatic carbons 348 from aliphatic carbons 349 """ 350 elm_a = at_a.autodock_element 351 charge_a = at_a.charge 352 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 353 vol_a = self.Vols.get(elm_a, 33.5103) 354 sp_a = self.Solpars.get(elm_a, -0.00110) # dict with default 355 #NB we hope we have checked for AtVol by this point 356 # in the MolecularSystem._check_interface 357 elm_b = at_b.autodock_element 358 charge_b = at_b.charge 359 #WHAT SHOULD BE A DEFAULT VOLUME? CARBON? 360 vol_b = self.Vols.get(elm_b, 33.5103) 361 sp_b = self.Solpars.get(elm_b, -0.00110) # dict with default 362 363 lig_energy = 0.0 364 rec_energy = 0.0 365 366 sol_func = self.get_ddsol(dist) 367 # 368 #NOTE how following line is truncated: 369 #lig_at_energy = (sp_b + self.solpar_q*abs(charge_b))*vol_a * sol_func 370 #this portion goes into the atom map 371 #lig_at_energy = (sp_b)*vol_a * sol_func + 372 #rec_at_energy = (sp_a + self.solpar_q*abs(charge_a))*vol_b * sol_func 373 374 # for the .d.map: 375 lig_at_energy = self.solpar_q*abs(charge_b)*vol_a * sol_func 376 ####this was wrong but passed because atom had charge 1.0 377 ###lig_at_energy = self.solpar_q*abs(charge_b)*vol_a * sol_func 378 lig_at_energy = self.solpar_q * vol_a * sol_func 379 return lig_at_energy # kcal/mol

380 381 NewDesolvationDesolvMap = NewDesolvationDesolvMapRefImpl 382 383 #end NewDesolvation 384 385 if __name__ == '__main__': 386 print "run tests in Tests.test_desolvation" 387

Source Code for Module PyAutoDock.desolvation