2

我正在尝试对循环反应中的反应障碍进行图形可视化（参见柠檬酸循环）。

我有两条反应路径的障碍，通过不同的路径关闭相同的反应循环。对于这个用例，一个简单的 2D 绘图变得不直观，所以我给了 matplotlib 的 Axes3D 对象一个镜头：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import (absolute_import, division,
                        print_function, unicode_literals)

from math import sin, cos

import numpy as np

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm


def z_cycle(mesh, origin, cycle, r, angle_offset, ring_width=0.5):
    disc_radius = r/(len(cycle)+1)
    angles = np.linspace(0, 2*np.pi, len(cycle), endpoint=False)+angle_offset
    mesh_x, mesh_y = mesh
    xo, yo = origin
    z = np.zeros_like(mesh_x)
    for a, d in zip(angles, cycle):
        xc, yc = xo + r/2*cos(a), yo + r/2*sin(a)
        z += d*((mesh_x-xc)**2 + (mesh_y-yc)**2 < (0.2*r)**2)
    z0_mask = np.where(z == 0)
    ring = np.logical_and(
        (mesh_x-xo)**2+(mesh_y-yo)**2 < (r/2+ring_width/2*disc_radius)**2,
        (mesh_x-xo)**2+(mesh_y-yo)**2 > (r/2-ring_width/2*disc_radius)**2
    )
    z += min(cycle)/2*ring*(z == 0)
    return z


def main(data='1,3,2,4,3,4,2;1,3,1.5,2', r=1.0):
    """
    Plot energy barriers for reaction cycles in 3d

    Cycles as separated by semicolon, and individual
    barrier heights are separated by comma. Cycles need
    to share the value of their respective first barrier height.
    """

    max_barrier = None
    cycles = []
    max_ = lambda a, b: max(a, b) if b is not None else a
    for cycle_string in data.split(';'):
        cycle = []
        for barrier in map(float, cycle_string.split(',')):
            max_barrier = max_(barrier, max_barrier)
            cycle.append(barrier)
        cycles.append(cycle)

    assert all(cycle[0] == cycles[0][0] for cycle in cycles[1:])

    grid1d = np.linspace(-2*r, 2*r, 200)
    mesh_x, mesh_y = mesh = np.meshgrid(grid1d, grid1d)
    mesh_z = np.zeros_like(mesh_x)

    angles = np.linspace(0, 2*np.pi, len(cycles), endpoint=False)
    for angle, cycle in zip(angles, cycles):
        # cycle is an iterable of numerical values corresponding to
        # barrier heights
        disc_radius = r/(len(cycle)+1)
        centre = ((r/2+disc_radius)*cos(angle), (r/2+disc_radius)*sin(angle))
        mesh_z += z_cycle(mesh, centre, cycle, r, np.pi+angle)

    # Plot
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(mesh_x, mesh_y, mesh_z, rstride=1, cstride=1,
                    cmap=cm.coolwarm, linewidths=0)

    plt.savefig('barriers3d.png')
    plt.show()

if __name__ == '__main__':
    main()

这给了我：

循环圆形 3dbar 图

这种方法的问题包括：

  • 由于需要非常密集的网格而导致性能不佳（旋转变得缓慢）
  • 图形人工制品

有谁知道解决这些问题的好方法？

我尝试使用手动生成表面，mpl_toolkits.mplot3d.art3d.Poly3DCollection但它的文档很少。

4

0 回答 0