Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
executable file 113 lines (94 sloc) 4.08 KB
# sasa.py computes the solvent accessible surface area
# Written by Chris Lockhart in Python3
from LLTK.structure import Structure
import numpy as np
from sklearn.metrics.pairwise import pairwise_distances
import sys
# Info: We are computing the solvent accessible surface area (SASA) for a
# Structure. Adapted from VMD code, an implementation of Shrake-Rupley 73
# Input: A Structure object, the list restrict that limits our output to a
# certain subset of the Structure. The probe radius = 1.4 A by default.
# The number of random points (npts) we use to check the SASA for an
# atom. The random number generator seed (seed). The boolean simple
# specifies if the same random points are used for all atoms.
# Output: The area, a float.
# To Do: I would like to make this code able to handle PBC
def sasa(structure,restrict=[],probe=1.4,npts=500,seed=-1,simple=True):
# Handle restrict
if isinstance(restrict,int): restrict=[restrict]
if isinstance(restrict,(list,tuple,range)): restrict=np.array(restrict)
# Check types
if not isinstance(structure,Structure):
print('Error: structure must be a Structure object.')
sys.exit(3)
if (not isinstance(restrict,np.ndarray) or
any(not isinstance(i,np.int32) for i in restrict)):
print('Error: restrict must be an array of integers.')
sys.exit(3)
if not isinstance(probe,float):
print('Error: must pass float (default: 1.4 A) to probe.')
sys.exit(3)
if not isinstance(npts,int):
print('Error: must pass integer (default: 500) to npts.')
sys.exit(3)
if not isinstance(seed,int):
print('Error: must pass integer to seed random number generator.')
sys.exit(3)
# If restrict is empty, it means we will compute for all atoms
if restrict.size == 0: restrict=np.array(range(structure.num_atom))
# Make sure that restrict is within bounds of Structure
if restrict.min() < 0 or restrict.max() >= structure.num_atom:
print('Error: index out of bounds.')
sys.exit(3)
# Make sure npts is positive
if npts <= 0:
print('Error: npts must be positive.')
sys.exit(3)
# Set seed for random number generator
np.random.seed() if seed < 0 else np.random.seed(seed)
# Get coordinates and radii (inflated by probe) for substructure
coord=structure.get_coord()
radius=structure.get('radius')+probe
# Build pairlist of atoms within 2*radius.max()
pairlist=pairwise_distances(coord)<=2.*radius.max()
# Fill pairlist diagonal with False so an atom doesn't interact with itself
np.fill_diagonal(pairlist,False)
# Define random points around arbitrary atom (same for all atoms)
if simple:
u1=np.random.rand(npts)
u2=np.random.rand(npts)
z=2.*u1-1.
phi=2.*np.pi*u2
r=np.sqrt(1.-z*z)
pts=np.zeros((npts,3))
pts[:,0]=r*np.cos(phi)
pts[:,1]=r*np.sin(phi)
pts[:,2]=z
# Finally, compute sasa
# For large number of atoms, this is embarrasingly slow. It needs to be
# optimized.
area=0.
for index in restrict:
# Obtain indices of those atoms paired to atom with index
pairs=pairlist[index,:].nonzero()[0]
# If ipairs is empty, we can skip further execution
if pairs.size == 0: continue
# Define random points around atom
if not simple:
u1=np.random.rand(npts)
u2=np.random.rand(npts)
z=2.*u1-1.
phi=2.*np.pi*u2
r=np.sqrt(1.-z*z)
pts=np.zeros((npts,3))
pts[:,0]=r*np.cos(phi)
pts[:,1]=r*np.sin(phi)
pts[:,2]=z
# Count points not impacted by paired atoms
r=pairwise_distances(coord[index,:]+radius[index]*pts,coord[pairs,:])
accessible=[j.min() for j in (r>radius[pairs])]
surfpts=np.sum(accessible)
# Increment area
area=area+surfpts*radius[index]**2
# return computed sasa
return 4.*np.pi*area/npts