In [1]:

    

AstexDiv_Dir = '/Users/dminh/clusters/CCB/AstexDiv_xtal'
import argparse
parser = argparse.ArgumentParser(description='Calculates a desolvation grids')
parser.add_argument('--prmtop_FN', \
  default = AstexDiv_Dir + '/1-build/1tow/receptor.prmtop', \
  help='Input AMBER PRMTOP file')
parser.add_argument('--inpcrd_FN', \
  default = AstexDiv_Dir + '/3-grids/1tow/receptor.trans.inpcrd', \
  help='Input coordinates')
parser.add_argument('--header_FN', \
  default = AstexDiv_Dir + '/3-grids/1tow/header_coarse.dx', \
  help='Input grid header (optional)')
parser.add_argument('--grid_FN', \
  default='desolv.dx', \
  help='Output for desolvation grid')
parser.add_argument('--probe_radius', default=1.4, \
  help='Radius of the solvent probe, in A')
parser.add_argument('--ligand_atom_radius', default=1.4, \
  help='Radius of the ligand atom, in A')
parser.add_argument('--SAS_points', default=1000, \
  help='Number of points on solvent accessible surface per receptor atom')
parser.add_argument('--integration_cutoff', default=10,
  help='Numerical integration cutoff, in A')
parser.add_argument('--spacing', nargs=3, type=float, \
  help='Grid spacing (overrides header)')
parser.add_argument('--counts', nargs=3, type=int, \
  help='Number of point in each direction (overrides header)')
parser.add_argument('-f')
args = parser.parse_args()

import desolvationGrid
self = desolvationGrid.desolvationGridCalculation(**vars(args))
self.calc_receptor_SAS_points()
self.calc_receptor_MS()


    

    




Reading header from /Users/dminh/clusters/CCB/AstexDiv_xtal/3-grids/1tow/header_coarse.dx
*** Files and parameters ***
Input AMBER prmtop      :	/Users/dminh/clusters/CCB/AstexDiv_xtal/1-build/1tow/receptor.prmtop
Input AMBER inpcrd      :	/Users/dminh/clusters/CCB/AstexDiv_xtal/3-grids/1tow/receptor.trans.inpcrd
Input grid header file  :	/Users/dminh/clusters/CCB/AstexDiv_xtal/3-grids/1tow/header_coarse.dx
Output grid             :	desolv.dx
Grid spacing            :	[ 0.5  0.5  0.5]
Grid counts             :	[73 73 73]

Finding receptor SAS points
 in 6.35 s
Determining the number of SAS points marking each grid point
 in 3.78 s

In [14]:

    

#     self.desolvationGrid = calc_desolvationGrid(self.receptor_MS_grid, \
#       self.kwargs['spacing'], self.kwargs['counts'], \
#       self.receptor_SAS_points, self.crd, \
#       SAS_sphere_pts, self.LJ_r2, max(np.sqrt(self.LJ_r2)), \
#       self.kwargs['ligand_atom_radius'], \
#       self.kwargs['probe_radius'], self.kwargs['integration_cutoff'])

#           desolvationGrid[i,j,k] = fraction_r4inv_low_dielectric(grid_c, \
#             spacing, counts, grid_point_x, grid_point_y, grid_point_z, \
#             ligand_atom_radius, integration_cutoff)

# cpdef fraction_r4inv_low_dielectric(\
#     int_t[:,:,:] grid, \
#     float_t[:] spacing, \
#     int_t[:] counts, \
#     float_t point_x, \
#     float_t point_y, \
#     float_t point_z, \
#     float_t r_min, \
#     float_t r_max)

grid = self.receptor_MS_grid
spacing = self.kwargs['spacing']
counts = self.kwargs['counts']

r_min = self.kwargs['probe_radius']
r_max = self.kwargs['integration_cutoff']*2

point_x = 30*spacing[0]
point_y = 25*spacing[1]
point_z = 25*spacing[2]

I_low_dielectric = 0.
I_total = 0.

i_min = max(int((point_x-r_max)/spacing[0]),0)
i_max = min(int((point_x+r_max)/spacing[0])+1,counts[0])
j_min = max(int((point_y-r_max)/spacing[1]),0)
j_max = min(int((point_y+r_max)/spacing[1])+1,counts[1])
k_min = max(int((point_z-r_max)/spacing[2]),0)
k_max = min(int((point_z+r_max)/spacing[2])+1,counts[2])

r_min2 = r_min*r_min
r_max2 = r_max*r_max
i = i_min
while i<i_max:
  dx  = point_x-i*spacing[0]
  dx2 = dx*dx
  j = j_min
  while j<j_max:
    dy  = point_y-j*spacing[1]
    dy2 = dy*dy
    dx2dy2 = dx2 + dy2
    if dx2dy2 < r_max2:
      k = k_min
      while k<k_max:
        dz  = point_z-k*spacing[2]
        dz2 = dz*dz
        r2 = dx2dy2 + dz2
        if (r2 < r_max2) and (r2 > r_min2):
          r4inv = 1/(r2*r2)
          if grid[i,j,k]<1:
            I_low_dielectric += r4inv
          I_total += r4inv
        k += 1
    j += 1
  i += 1


    

In [15]:

    

print 'I_low_dielectric\t', I_low_dielectric
print 'I_total         \t', I_total

fD = I_low_dielectric/I_total
print 'Fraction Desolvated\t', fD


    

    




I_low_dielectric	3.55904423027
I_total         	68.6321151429
Fraction Desolvated	0.0518568344115

In [16]:

    

# Previous Igrid
fD * (1/r_min - 1/r_max)


    

    
Out[16]:





0.03704059600821101

In [20]:

    

# A simpler numerical integral
import math
dV = spacing[0]*spacing[1]*spacing[2]
I_low_dielectric*dV/(4*math.pi)


    

    
Out[20]:





0.035402467620679064

In [18]:

    

I_low_dielectric*dV


    

    
Out[18]:





0.44488052878430351

In [19]:

    

fD * (1/r_min - 1/r_max) * 4 * math.pi


    

    
Out[19]:





0.46546585721593248

In [ ]: