In [1]:
execfile('desolvationGrid.py')
self.receptor_SAS_points()


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 SAS points
 in 6.18 s

In [2]:
self.receptor_MS()
self.full_MS()


Determining whether grid points are high or low dielectric
 in 156.68 s
 in 0.31 s

In [3]:
# TODO: Numerical integral of r^{-4}

In [4]:
del set_inside_sphere_to

In [5]:
%load_ext Cython

In [13]:
%%cython

import cython

import numpy as np
cimport numpy as np

from cpython cimport bool

ctypedef np.float_t float_t
ctypedef np.int_t int_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cpdef set_inside_sphere_to(\
    int[:,:,:] grid, \
    double[:] spacing, \
    long[:] counts, \
    double[:] point, \
    double r, \
    int val):

  cdef np.ndarray[np.int_t] lower_lim, upper_lim
  cdef int d, i, j, k
  cdef double r2
  
  lower_lim = np.zeros((3,), dtype=int)
  upper_lim = np.zeros((3,), dtype=int)
  
  for d in [0,1,2]:
    lower_lim[d] = max(int(np.floor((point[d]-r)/spacing[d])),0)
    upper_lim[d] = min(int(np.ceil((point[d]+r)/spacing[d])),counts[d])
    
  r2 = r*r
  for i in range(lower_lim[0],upper_lim[0]):
    for j in range(lower_lim[1],upper_lim[1]):
      for k in range(lower_lim[2],upper_lim[2]):
        if (np.square(point[0]-i*spacing[0]) + \
            np.square(point[1]-j*spacing[1]) + \
            np.square(point[2]-k*spacing[2]))<r2:
          grid[i,j,k]=val

# This is the original python code
# def set_inside_sphere_to(grid, spacing, counts, point, r, val):
#   lower_lim = [max(int(np.floor((point[d]-r)/spacing[d])),0) \
#     for d in range(3)]
#   upper_lim = [min(int(np.ceil((point[d]+r)/spacing[d])),counts[d]) \
#     for d in range(3)]
#   r2 = r*r
#   for i in range(lower_lim[0],upper_lim[0]):
#     for j in range(lower_lim[1],upper_lim[1]):
#       for k in range(lower_lim[2],upper_lim[2]):
#         if (np.square(point[0]-i*spacing[0]) + \
#             np.square(point[1]-j*spacing[1]) + \
#             np.square(point[2]-k*spacing[2]))<r2:
#           grid[i,j,k]=val


Error compiling Cython file:
------------------------------------------------------------
...

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cpdef set_inside_sphere_to(\
    np.bool_[:,:,:] grid, \
   ^
------------------------------------------------------------

/Users/dminh/.ipython/cython/_cython_magic_39736fd3ab6823e4b259c85630e79712.pyx:16:4: 'bool_' is not a type identifier

Error compiling Cython file:
------------------------------------------------------------
...
    np.bool_[:,:,:] grid, \
    double[:] spacing, \
    long[:] counts, \
    double[:] point, \
    double r, \
    np.bool_ val):
   ^
------------------------------------------------------------

/Users/dminh/.ipython/cython/_cython_magic_39736fd3ab6823e4b259c85630e79712.pyx:21:4: 'bool_' is not a type identifier

In [12]:
SAS_radii = self.LJ_radii + self.kwargs['probe_radius']
    self.receptor_MS_grid = np.ones(shape=tuple(self.kwargs['counts']), \
      dtype=bool)
    # Tentatively assign the grid inside the SAS to low dielectric
    for atom_index in range(5):
      set_inside_sphere_to(self.receptor_MS_grid, self.kwargs['spacing'], \
        self.kwargs['counts'], self.crd[atom_index,:], \
        SAS_radii[atom_index], False)


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-fca2a3941a69> in <module>()
      3 # Tentatively assign the grid inside the SAS to low dielectric
      4 for atom_index in range(5):
----> 5   set_inside_sphere_to(self.receptor_MS_grid, self.kwargs['spacing'],     self.kwargs['counts'], self.crd[atom_index,:],     SAS_radii[atom_index], False)

_cython_magic_c9f61556980a2cd58655eca08dc3bbb8.pyx in _cython_magic_c9f61556980a2cd58655eca08dc3bbb8.set_inside_sphere_to (/Users/dminh/.ipython/cython/_cython_magic_c9f61556980a2cd58655eca08dc3bbb8.c:2741)()

ValueError: Does not understand character buffer dtype format string ('?')

In [ ]: