surface

Isosurfaces for molecules in crystals or isolation.

promolecule_density_isosurface(promol, isovalue=0.002, sep=0.2, props=True, extra_props=None, smoothing='laplacian')

Calculate the promolecule density isosurface.

Parameters:

Name	Type	Description	Default
isovalue	float	level set value for the isosurface (default=0.002) in au.	0.002
separation	float	separation between density grid used in the surface calculation (default 0.2) in Angstroms.	required
props	bool	calculate surface properties	True
extra_props	dict	dictionary of property names and functions of vertex positions for the calculation of additional properties on the isosurface	None

Returns: IsosurfaceMesh: a namedtuple of vertices, faces, face normals and vertex properties

Source code in chmpy/surface.py

def promolecule_density_isosurface(
    promol, isovalue=0.002, sep=0.2, props=True, extra_props=None, smoothing="laplacian"
):
    """Calculate the promolecule density isosurface.

    Args:
        isovalue (float, optional): level set value for the isosurface
            (default=0.002) in au.
        separation (float, optional):  separation between density grid
                used in the surface calculation (default 0.2) in Angstroms.
        props (bool, optional): calculate surface properties
        extra_props(dict, optional): dictionary of property names and functions
            of vertex positions for the calculation of additional properties
            on the isosurface
    Returns:
        IsosurfaceMesh: a namedtuple of vertices, faces,
            face normals and vertex properties
    """
    t1 = time.time()
    extras = {}
    if SURFACE_BACKEND == "occpy":
        verts, faces, normals, extras = _promolecule_mesh_occpy(promol, isovalue, sep)
    else:
        l, u = promol.bb()
        x_grid = np.arange(l[0], u[0], sep, dtype=np.float32)
        y_grid = np.arange(l[1], u[1], sep, dtype=np.float32)
        z_grid = np.arange(l[2], u[2], sep, dtype=np.float32)
        x, y, z = np.meshgrid(x_grid, y_grid, z_grid)
        shape = x.shape
        pts = np.c_[x.ravel(), y.ravel(), z.ravel()]
        d = promol.rho(pts).reshape(shape)
        verts, faces, normals, _ = marching_cubes(
            d, isovalue, spacing=(sep, sep, sep), gradient_direction="descent"
        )
        verts = np.c_[verts[:, 1], verts[:, 0], verts[:, 2]] + l
    LOG.debug("Surface centroid: %s", np.mean(verts, axis=0))
    LOG.debug("Mol centroid: %s", np.mean(promol.positions, axis=0))
    LOG.debug("Max (x,y,z): %s", np.max(verts, axis=0))
    LOG.debug("Min (x,y,z): %s", np.min(verts, axis=0))

    if smoothing == "laplacian":
        nverts_pre = len(verts)
        verts, faces = smooth_laplacian(verts, faces)
        if len(verts) != nverts_pre:
            extras = {}

    vertex_props = dict(extras)

    if props:
        if extra_props is not None:
            for k, func in extra_props.items():
                LOG.debug("Calculating additional surface property: %s", k)
                vertex_props[k] = func(verts)
        d_i, d_norm_i, vecs = promol.d_norm(verts)
        vertex_props["d_i"] = d_i
        vertex_props["d_norm_i"] = d_norm_i
        LOG.debug("d_i (min, max): (%.2f, %.2f)", np.min(d_i), np.max(d_i))
    t2 = time.time()
    LOG.info(
        "promolecule surface took %.3fs (%s, %d verts)",
        t2 - t1, SURFACE_BACKEND, len(verts),
    )
    return IsosurfaceMesh(verts, faces, normals, vertex_props)

smooth_laplacian(vertices, faces, **kwargs)

Smooth vertices and faces using a Laplacian filter

Source code in chmpy/surface.py

def smooth_laplacian(vertices, faces, **kwargs):
    """Smooth vertices and faces using a Laplacian filter"""
    Trimesh = require("trimesh", "smoothing an isosurface").Trimesh
    filter_humphrey = require(
        "trimesh.smoothing", "smoothing an isosurface"
    ).filter_humphrey

    kwargs.setdefault("iterations", 2)
    mesh = Trimesh(vertices, faces)
    filter_humphrey(mesh, **kwargs)
    return mesh.vertices, mesh.faces

stockholder_weight_isosurface(s, isovalue=0.5, sep=0.2, props=True, extra_props=None, smoothing='laplacian')

Calculate stockholder weight (Hirshfeld) surface.

Parameters:

Name	Type	Description	Default
isovalue	float	level set value for the isosurface (default=0.5), dimensionless.	0.5
separation	float	separation between density grid used in the surface calculation (default 0.2) in Angstroms.	required
props	bool	calculate surface properties	True
extra_props	dict	dictionary of property names and functions of vertex positions for the calculation of additional properties on the isosurface	None

Returns: IsosurfaceMesh: a namedtuple of vertices, faces, face normals and vertex properties

Source code in chmpy/surface.py

def stockholder_weight_isosurface(
    s, isovalue=0.5, sep=0.2, props=True, extra_props=None, smoothing="laplacian"
):
    """Calculate stockholder weight (Hirshfeld) surface.

    Args:
        isovalue (float, optional): level set value for the isosurface
            (default=0.5), dimensionless.
        separation (float, optional):  separation between density grid
                used in the surface calculation (default 0.2) in Angstroms.
        props (bool, optional): calculate surface properties
        extra_props(dict, optional): dictionary of property names and functions
            of vertex positions for the calculation of additional properties
            on the isosurface
    Returns:
        IsosurfaceMesh: a namedtuple of vertices, faces,
            face normals and vertex properties
    """
    t1 = time.time()
    extras = {}
    if SURFACE_BACKEND == "occpy":
        verts, faces, normals, extras = _stockholder_mesh_occpy(s, isovalue, sep)
    else:
        l, u = s.bb()
        x_grid = np.arange(l[0], u[0], sep, dtype=np.float32)
        y_grid = np.arange(l[1], u[1], sep, dtype=np.float32)
        z_grid = np.arange(l[2], u[2], sep, dtype=np.float32)
        x, y, z = np.meshgrid(x_grid, y_grid, z_grid)
        separations = np.array((sep, sep, sep))
        shape = x.shape
        pts = np.c_[x.ravel(), y.ravel(), z.ravel()]
        pts = np.array(pts, dtype=np.float32)
        weights = s.weights(pts).reshape(shape)
        verts, faces, normals, _ = marching_cubes(
            weights, isovalue, spacing=separations, gradient_direction="descent"
        )
        verts = np.c_[verts[:, 1], verts[:, 0], verts[:, 2]] + l

    if smoothing == "laplacian":
        nverts_pre = len(verts)
        verts, faces = smooth_laplacian(verts, faces)
        if len(verts) != nverts_pre:
            extras = {}

    LOG.debug("Surface centroid: %s", np.mean(verts, axis=0))
    LOG.debug("Mol centroid: %s", np.mean(s.dens_a.positions, axis=0))
    LOG.debug("Max (x,y,z): %s", np.max(verts, axis=0))
    LOG.debug("Min (x,y,z): %s", np.min(verts, axis=0))
    vertex_props = dict(extras)
    if props:
        if extra_props is not None:
            for k, func in extra_props.items():
                LOG.debug("Calculating additional surface property: %s", k)
                vertex_props[k] = func(verts)
        d_i, d_e, d_norm_i, d_norm_e, dp, angles = s.d_norm(verts)
        vertex_props["d_i"] = d_i
        vertex_props["d_e"] = d_e
        vertex_props["d_norm_i"] = d_norm_i
        vertex_props["d_norm_e"] = d_norm_e
        d_norm = d_norm_i + d_norm_e
        vertex_props["d_norm"] = d_norm
        vertex_props["dp"] = dp
        vertex_props["angle"] = np.abs(angles)
        LOG.debug("d_norm (min, max): (%.2f, %.2f)", np.min(d_norm), np.max(d_norm))
    t2 = time.time()
    LOG.info(
        "stockholder weight surface took %.3fs (%s, %d verts)",
        t2 - t1, SURFACE_BACKEND, len(verts),
    )
    return IsosurfaceMesh(verts, faces, normals, vertex_props)