SCHWED-2349: use faster call directly to fortran binary of naccess.

modules/bindings/pymod/naccess.py

@@ -30,16 +30,35 @@ def _GetExecutable(naccess_exe):
   """
   return settings.Locate('naccess', explicit_file_name=naccess_exe)
 
+def _CheckNaccessRoot(naccess_root):
+  """
+  :return: True, if given directory contains "accall" binary and files
+           "vdw.radii" and "standard.data".
+  :param naccess_root: Path to naccess folder to check.
+  """
+  accall_exe = os.path.join(naccess_root, "accall")
+  check = (os.path.exists(accall_exe) and os.access(accall_exe, os.X_OK) \
+           and os.path.exists(os.path.join(naccess_root, "vdw.radii")) \
+           and os.path.exists(os.path.join(naccess_root, "standard.data")))
+  if not check:
+    LogWarning("NACCESS: Could not find required files to launch accall " \
+               "directly in %s." % naccess_root)
+  return check
+
 def _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms):
   """
   Setup files for naccess calculation in temporary directory
 
-  :param entity:      EntityHandle or EntityView to calculate surface
-  :param selection:   Calculate surface for subset of entity
-  :param scratch_dir: Directory for temporary files (NACCESS is sensitive to "." in directory names
-  :param max_number_of_atoms: Max Number of atoms in the entity (i.e. is limited in the default NACCESS version to 50 000)
-  :returns:            array containing temporary directory, input filename for naccess and directory of the input file
-  :exception:         RuntimeError if selection is not valid
+  :param entity:              OST entity to calculate surface
+  :param selection:           Calculate surface for subset of entity
+  :param scratch_dir:         Scratch directory. A subfolder for temporary files
+                              is created in there. If not specified, a default
+                              directory is used (see :func:`tempfile.mkdtemp`).
+  :param max_number_of_atoms: Max Number of atoms in the entity (i.e. is limited
+                              in the default NACCESS version to 50 000)
+  :returns: Tuple containing temporary directory, input filename for naccess and
+            directory of the input file
+  :exception: RuntimeError if selection is not valid or too many atoms
   """
 
   # create temporary directory
@@ -103,7 +122,8 @@ def _ParseAsaFile(entity, file, asa_atom):
       if(a.IsValid()):
         a.SetFloatProp(asa_atom, float(asa))
       else:
-        print chain_id, resNum, atom_name
+        LogWarning("NACCESS: invalid asa entry %s %s %s" \
+                   % (chain_id, resNum, atom_name))
       
 def _ParseRsaFile(enti,file, asa_abs, asa_rel):
   """
@@ -161,33 +181,62 @@ def __CleanupFiles(dir_name):
   import shutil
   shutil.rmtree(dir_name)
 
+
+def _RunACCALL(command, temp_dir, query):
+  """
+  Fast method to run the Naccess surface calculation.
+
+  This method starts the accall binary directly and pipes in the input provided
+  in *query*. This is faster than calling the "naccess" script since the script
+  has a constant overhead of roughly 1.3s in each call.
+
+  :param command:  Command to execute
+  :param temp_dir: Command is executed with this working directory
+  :param query:    User input to pipe into *command*
+  :returns:        stdout of command
+  :exception:      CalledProcessError for non-zero return value
+  """
+  proc = subprocess.Popen(command, stdout=subprocess.PIPE,
+                          stderr=subprocess.PIPE, stdin=subprocess.PIPE,
+                          cwd=temp_dir)
+  stdout_value, stderr_value = proc.communicate(query)
+
+  # check for successful completion of naccess
+  if proc.returncode != 0:
+    LogWarning("WARNING: naccess error\n%s\n%s" % (stdout_value, stderr_value))
+    raise subprocess.CalledProcessError(proc.returncode, command)
+
+  return stdout_value
+
+
 def _RunNACCESS(command, temp_dir):
   """
-  Method to run the MSMS surface calculation
+  Method to run the Naccess surface calculation.
 
-  This method starts the external MSMS executable and returns the stdout of MSMS
+  This method starts the external Naccess executable and returns the stdout.
 
-  :param command:          Command to execute
-  :returns:                stdout of MSMS
-  :exception:              CalledProcessError for non-zero return value
+  :param command:  Command to execute
+  :param temp_dir: Command is executed with this working directory
+  :returns:        stdout of command
+  :exception:      CalledProcessError for non-zero return value
   """
   proc = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE,
                           cwd=temp_dir)
   stdout_value, stderr_value = proc.communicate()
 
   # check for successful completion of naccess
-  if proc.returncode!=0:
-    print "WARNING: naccess error\n", stdout_value
+  if proc.returncode != 0:
+    LogWarning("WARNING: naccess error\n%s" % stdout_value)
     raise subprocess.CalledProcessError(proc.returncode, command)
 
   return stdout_value
-  
+
 
 def CalculateSurfaceArea(entity,  radius=1.4,  
                          include_hydrogens=False, include_hetatm = False, 
                          include_water=False, selection="", naccess_exe=None,
-                         keep_files=False, asa_abs= "asaAbs", asa_rel="asaRel",
-                         asa_atom="asaAtom", scratch_dir=None,
+                         naccess_root=None, keep_files=False, asa_abs= "asaAbs",
+                         asa_rel="asaRel", asa_atom="asaAtom", scratch_dir=None,
                          max_number_of_atoms=50000):
   """
   Calculates analytical the solvent accessible surface area by using the
@@ -206,6 +255,9 @@ def CalculateSurfaceArea(entity,  radius=1.4,
   :param include_water:       Calculate surface including water
   :param selection:           Calculate surface for subset of entity
   :param naccess_exe:         naccess executable (full path to executable)
+  :param naccess_root:        Path to folder containing "accall" binary and
+                              files "vdw.radii" and "standard.data". This is the
+                              fastest way to call naccess!
   :param keep_files:          If True, do not delete temporary files
   :param asa_abs:             Attaches per residue absolute SASA to specified
                               FloatProp on residue level
@@ -223,38 +275,70 @@ def CalculateSurfaceArea(entity,  radius=1.4,
   """
   
   # check if naccess executable is specified
-  naccess_executable=_GetExecutable(naccess_exe)
+  if naccess_root and _CheckNaccessRoot(naccess_root):
+    # use faster, direct call to accall binary
+    fast_mode = True
+  else:
+    # get naccess executable
+    naccess_executable = _GetExecutable(naccess_exe)
+    # see if we can extract naccess_root from there (fallback to old mode)
+    naccess_root = os.path.dirname(naccess_executable)
+    fast_mode = _CheckNaccessRoot(naccess_root)
   
   # setup files for naccess
   (naccess_data_dir, naccess_data_file, naccess_data_base) \
     = _SetupFiles(entity, selection, scratch_dir, max_number_of_atoms)
 
-  # set command line
-  command = "%s %s -p %f " % \
-            (naccess_executable, naccess_data_file, radius)
-  if include_hydrogens:
-    command = "%s -w" % command
-  if include_water:
-    command = "%s -y" % command
-  if include_hetatm:
-    command = "%s -h" % command
-  # execute naccess
-  stdout_value = _RunNACCESS(command, naccess_data_dir)
-  
-  # parse outout
-  new_asa= os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base) 
-  _ParseAsaFile(entity, new_asa, asa_atom)
+  try:
+    # call naccess
+    if fast_mode:
+      # cook up stdin query (same logic as naccess script)
+      query = "PDBFILE %s\n" \
+              "VDWFILE %s\n" \
+              "STDFILE %s\n" \
+              "PROBE %f\n" \
+              "ZSLICE 0.05\n" \
+              % (naccess_data_file, os.path.join(naccess_root, "vdw.radii"),
+                 os.path.join(naccess_root, "standard.data"), radius)
+      if include_hydrogens:
+        query += "HYDROGENS\n"
+      if include_water:
+        query += "WATERS\n"
+      if include_hetatm:
+        query += "HETATOMS\n"
+      # call it
+      command = os.path.join(naccess_root, "accall")
+      _RunACCALL(command, naccess_data_dir, query)
+    else:
+      LogWarning("NACCESS: Falling back to slower call to %s." \
+                 % naccess_executable)
+      # set command line
+      command = "%s %s -p %f " % \
+                (naccess_executable, naccess_data_file, radius)
+      if include_hydrogens:
+        command = "%s -w" % command
+      if include_water:
+        command = "%s -y" % command
+      if include_hetatm:
+        command = "%s -h" % command
+      # execute naccess
+      _RunNACCESS(command, naccess_data_dir)
     
-  new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base) 
-  _ParseRsaFile(entity, new_rsa, asa_abs, asa_rel)
+    # parse outout
+    new_asa = os.path.join(naccess_data_dir, "%s.asa" % naccess_data_base) 
+    _ParseAsaFile(entity, new_asa, asa_atom)
+      
+    new_rsa = os.path.join(naccess_data_dir, "%s.rsa" % naccess_data_base) 
+    _ParseRsaFile(entity, new_rsa, asa_abs, asa_rel)
+
+  finally:
+    # clean up
+    if not keep_files:
+      __CleanupFiles(naccess_data_dir)
   
   # sum up Asa for all atoms
   sasa = 0.0 
   for a in entity.atoms:
     sasa += a.GetFloatProp(asa_atom, 0.0)
-    
-  # clean up
-  if not keep_files:
-    __CleanupFiles(naccess_data_dir)
 
   return sasa