Skip to content

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, optional

dictionary of property names and functions of vertex positions for the calculation of additional properties on the isosurface

 required

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

Source code in chmpy/surface.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
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()
    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)
    separations = np.array((sep, sep, sep))
    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"
    )
    LOG.debug("Separation (x,y,z): %s", separations)
    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))
    vertex_props = {}

    if smoothing == "laplacian":
        verts, faces = smooth_laplacian(verts, faces)

    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, %d pts", t2 - t1, len(pts))
    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 
12
13
14
15
16
17
18
19
20
def smooth_laplacian(vertices, faces, **kwargs):
    """Smooth vertices and faces using a Laplacian filter"""
    from trimesh.smoothing import filter_humphrey
    from trimesh import Trimesh

    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, optional

dictionary of property names and functions of vertex positions for the calculation of additional properties on the isosurface

 required

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

Source code in chmpy/surface.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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()
    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"
    )

    if smoothing == "laplacian":
        verts, faces = smooth_laplacian(verts, faces)

    LOG.debug("Separation (x,y,z): %s", separations)
    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(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 = {}
    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, %d pts", t2 - t1, len(pts))
    return IsosurfaceMesh(verts, faces, normals, vertex_props)