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