Get void with OpenEye Spicoli

notes_on_OESpicoli_tests.md

I wanted to get the void volume with OpenEye Spicoli for use in OERocs.

But it did not create a usable surface.

This here is for my reference as I am missing something.

Standard reading:

# -----------------------------------
lig_filename = 'peptide.mol'
pdb_filename = 'template2.pdb'
surf_filename = 'void.oesrf'
max_neigh_distance = 4
# -----------------------------------

import os
import numpy as np
from pathlib import Path
# doubletap
os.environ['OE_LICENSE']=os.environ.get('OE_LICENSE', os.environ['HOME']+'/ASAP-oe_license.txt')
assert Path(os.environ['OE_LICENSE']).exists()
from openeye import oechem, oeomega, oespicoli, oegrid

# -----------------------------------

# read protein
pifs = oechem.oemolistream()
assert pifs.open(pdb_filename), f'Unable to open {pdb_filename} for reading as stream.'
prot = oechem.OEGraphMol()
assert oechem.OEReadMolecule(pifs, prot), f'Unable to read protein {pdb_filename}'
oechem.OEAddExplicitHydrogens(prot)
oechem.OEAssignBondiVdWRadii(prot)

# read ligand
lifs = oechem.oemolistream()
assert lifs.open(lig_filename), f'Unable to open {lig_filename} for reading as stream.'
lig = oechem.OEGraphMol()
assert oechem.OEReadMolecule(lifs, lig), f'Unable to read ligand {lig_filename}'
oechem.OEAddExplicitHydrogens(lig)
oechem.OEAssignBondiVdWRadii(lig)

But the problems happen.

Given a surf = oespicoli.OESurface() surface, I export to a MRC map for viewing in PyMol

def to_grid(surf: oespicoli.OESurface, grid_filename = 'void.ccp4')
  oespicoli.OEMakeGridFromSurface(grid, surf)
  oegrid.OEWriteGrid(grid_filename, grid)
  return grid.GetSize()

The most basic surface would be

surf = oespicoli.OESurface()
oespicoli.OEMakeMolecularSurface(surf, prot)
to_grid(surf)

Even with 773,054 points (~100 points per axis) it was very grainy and not the surface but the cavity-ish

This is nearly the same as

surf = oespicoli.OESurface()
oespicoli.OEMakeCavitySurfaces(prot, surf)
to_grid(surf)

or oddly

surf = oespicoli.OESurface()
oespicoli.OEMakeCavitySurfaces(prot, surf)
oespicoli.OEInvertSurface(surf)
to_grid(surf)

Now, some files on the internet do close variants of the follow in the surface form:

max_neigh_distance = 2

for i in range(surf.GetNumVertices()):
    vert = surf.GetVertex(i)
    for atom in lig.GetAtoms():
        dist = np.linalg.norm( np.array(lig.GetCoords(atom)) - np.array( vert))
        if dist <  max_neigh_distance:
            surf.SetVertexCliqueElement(i, 1)

oespicoli.OESurfaceCropToClique(surf, 1)

Having gone nowhere with this, I will run an OpenMM run for a water map and extract with MDAnalysis