From 836d644432663199c33b1f9a349657e7da5bf83b Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 16:36:42 +0100 Subject: [PATCH 01/18] Add function for plotting ball-and-stick models --- src/biotite/structure/graphics/atoms.py | 41 ++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/biotite/structure/graphics/atoms.py b/src/biotite/structure/graphics/atoms.py index 5129faa2c..bf790ea32 100644 --- a/src/biotite/structure/graphics/atoms.py +++ b/src/biotite/structure/graphics/atoms.py @@ -4,7 +4,7 @@ __name__ = "biotite.structure.graphics" __author__ = "Patrick Kunzmann" -__all__ = ["plot_atoms"] +__all__ = ["plot_atoms", "plot_ball_and_stick_model"] import numpy as np import matplotlib.pyplot as plt @@ -108,6 +108,45 @@ def plot_atoms(axes, atoms, colors, line_width=1.0, background_color=None, _set_box(axes, atoms.coord, center, size, zoom) +def plot_ball_and_stick_model(axes, atoms, colors, ball_size=200, + line_color="black", line_width=1.0, + background_color=None, center=None, + size=None, zoom=1.0): + if not isinstance(axes, Axes3D): + raise ValueError("The given axes mut be an 'Axes3D'") + if atoms.bonds is None: + raise ValueError("The atom array must have an associated bond list") + + # Calculating connections between atoms + line_coord = [ + (atoms.coord[index1], atoms.coord[index2]) + for index1, index2 in atoms.bonds.as_array()[:,:2] + ] + + # Plot sticks + # Use 'Line3DCollection' for higher efficiency + sticks = Line3DCollection( + line_coord, color=line_color, linewidths=line_width + ) + axes.add_collection(sticks) + + # Plot balls + axes.scatter( + *atoms.coord.T, s=ball_size, c=colors, linewidth=0, alpha=1 + ) + + # Set viewing angle + axes.azim = -90 + axes.elev = 90 + # Remove frame + axes.axis("off") + # Set background color + if background_color is not None: + axes.set_facecolor(background_color) + axes.get_figure().set_facecolor(background_color) + _set_box(axes, atoms.coord, center, size, zoom) + + def _set_box(axes, coord, center, size, zoom): """ This ensures an approximately equal aspect ratio in a 3D plot under From 52416ea4a6dee31ef4d83257c1a711999f971838 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 16:36:52 +0100 Subject: [PATCH 02/18] Small fixes --- src/biotite/structure/__init__.py | 2 +- src/biotite/structure/charges.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/biotite/structure/__init__.py b/src/biotite/structure/__init__.py index 71ff3b232..bb6b8f989 100644 --- a/src/biotite/structure/__init__.py +++ b/src/biotite/structure/__init__.py @@ -117,5 +117,5 @@ from .superimpose import * from .transform import * from .basepairs import * -from .charges.py import * +from .charges import * # util is used internally diff --git a/src/biotite/structure/charges.py b/src/biotite/structure/charges.py index 629ecb26f..3161ec246 100644 --- a/src/biotite/structure/charges.py +++ b/src/biotite/structure/charges.py @@ -13,7 +13,7 @@ __all__ = ["partial_charges"] import numpy as np -from .structure.info import residue +from .info import residue import warnings From 20ff411956dfec658981cc9b2d707cf3aa672de6 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 16:37:11 +0100 Subject: [PATCH 03/18] Remove 'pytest-cov' requirement --- environment.yml | 1 - pyproject.toml | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 0efd50526..282d75fb1 100644 --- a/environment.yml +++ b/environment.yml @@ -24,7 +24,6 @@ dependencies: - mdtraj >=1.9.3 - matplotlib >=3.3 - pytest >=5.2 - - pytest-cov >=2.5 - sphinx >=1.7 - sphinx-gallery =0.7.0 - numpydoc >=0.8 diff --git a/pyproject.toml b/pyproject.toml index bad552e68..c3fae7a1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,4 +5,4 @@ requires = [ "wheel >= 0.30", "numpy >= 1.13", "cython >= 0.29" -] \ No newline at end of file +] From 42e5216c98aba7fbecb6a9cb3dc57a081d0dc2fa Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 16:37:21 +0100 Subject: [PATCH 04/18] Remove comment --- doc/examples/scripts/structure/molecular_visualization.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/doc/examples/scripts/structure/molecular_visualization.py b/doc/examples/scripts/structure/molecular_visualization.py index 095d444d1..138515b61 100644 --- a/doc/examples/scripts/structure/molecular_visualization.py +++ b/doc/examples/scripts/structure/molecular_visualization.py @@ -44,8 +44,6 @@ colors[caffeine.element == "N"] = (0.0, 0.0, 0.8) # blue colors[caffeine.element == "O"] = (0.8, 0.0, 0.0) # red -# Create a quadratic figure to ensure a correct aspect ratio -# in the visualized molecule fig = plt.figure(figsize=(8.0, 8.0)) ax = fig.add_subplot(111, projection="3d") graphics.plot_atoms( From 2be794b86536e63afd2db144284bd7d75f60fe76 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 16:37:29 +0100 Subject: [PATCH 05/18] Add example script --- .../scripts/structure/peoe_visualization.py | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 doc/examples/scripts/structure/peoe_visualization.py diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py new file mode 100644 index 000000000..97fa7d0c9 --- /dev/null +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -0,0 +1,76 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.cm import ScalarMappable +from matplotlib.animation import FuncAnimation +import biotite.structure as struc +import biotite.structure.info as info +import biotite.structure.graphics as graphics + + +MOLECULE_NAME = "CFF" + +ITERATION_NUMBER = 6 +ELEMENT_FONT_SIZE = 10 +BALL_SCALE = 200 +RAY_SCALE = 3000 +CMAP_NAME = "bwr_r" + + +# Get an atom array for the selected molecule +# Caffeine has the PDB reside name 'CFF' +molecule = info.residue(MOLECULE_NAME) + +# TODO: Align molecule with PCA + +# Balls should be colored by partial charge +charges = struc.partial_charges(molecule, ITERATION_NUMBER) +# Later this variable stores values between 0 and 1 for use in color map +charge_in_color_range = charges.copy() +# Norm charge values to highest absolute value +charge_in_color_range /= np.max(np.abs(charge_in_color_range)) +# Transform range (-1, 1) to range (0,1) +charge_in_color_range = (charge_in_color_range + 1) / 2 +# Calculate colors +color_map = plt.get_cmap(CMAP_NAME) +colors = color_map(charge_in_color_range) + +# Ball size should be proportional to VdW radius of the respective atom +ball_sizes = np.array( + [info.vdw_radius_single(e) for e in molecule.element] +) * BALL_SCALE + + + +fig = plt.figure(figsize=(8.0, 16.0)) +ax = fig.add_subplot(111, projection="3d") +graphics.plot_ball_and_stick_model( + ax, molecule, colors, ball_size=ball_sizes, line_width=3, + line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 +) + +for atom in molecule: + ax.text( + *atom.coord, atom.element, + fontsize=ELEMENT_FONT_SIZE, color="black", + ha="center", va="center", zorder=100 + ) + +# Gradient of ray strength +ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE +ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 +ax.scatter( + *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 +) +ax.scatter( + *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 +) + +color_bar = fig.colorbar(ScalarMappable(cmap=color_map)) +color_bar.set_label("Partial charge", color="white") +color_bar.ax.yaxis.set_tick_params(color="white") +color_bar.outline.set_edgecolor("white") +for label in color_bar.ax.get_yticklabels(): + label.set_color("white") + +fig.tight_layout() +plt.show() \ No newline at end of file From b38e01bfc2affcce8ae45f0bd059970903f53b6e Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 17:53:16 +0100 Subject: [PATCH 06/18] Fix docstring --- src/biotite/structure/transform.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/biotite/structure/transform.py b/src/biotite/structure/transform.py index c6580c062..1e2213f8f 100644 --- a/src/biotite/structure/transform.py +++ b/src/biotite/structure/transform.py @@ -243,7 +243,7 @@ def align_vectors(atoms, origin_direction, target_direction, At first the transformation (translation and rotation), that is necessary to align the origin vector to the target vector is calculated. - This means means, that the application of the transformation on the + This means, that the application of the transformation on the origin vector would give the target vector. Then the same transformation is applied to the given atoms/coordinates. @@ -264,8 +264,8 @@ def align_vectors(atoms, origin_direction, target_direction, Returns ------- transformed : Atom or AtomArray or AtomArrayStack or ndarray, shape=(3,) or shape=(n,3) or shape=(m,n,3) - A copy of the input atoms or coordinates, rotated about the - given axis. + A copy of the input atoms or coordinates with the applied + transformation. See Also -------- From cda3020d025e0a2ff9cf1c1fe4f510c38f327a3d Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 18:01:10 +0100 Subject: [PATCH 07/18] Add scikit-learn as development dependency --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index 282d75fb1..82283854b 100644 --- a/environment.yml +++ b/environment.yml @@ -17,6 +17,7 @@ dependencies: - requests >=2.12 - numpy >=1.13 - scipy >=0.14 + - scikit-learn >= 0.18.0 - networkx >=2.0 - pydot >=1.4 - msgpack-python >=0.5.6 From e3f3f09516aa83442fa388295f37657ecc9e516d Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 18:01:43 +0100 Subject: [PATCH 08/18] Orient molecule --- doc/examples/scripts/structure/peoe_visualization.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 97fa7d0c9..8e50374fe 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -1,4 +1,5 @@ import numpy as np +from sklearn.decomposition import PCA import matplotlib.pyplot as plt from matplotlib.cm import ScalarMappable from matplotlib.animation import FuncAnimation @@ -20,7 +21,13 @@ # Caffeine has the PDB reside name 'CFF' molecule = info.residue(MOLECULE_NAME) -# TODO: Align molecule with PCA +# Align molecule with principal component analysis: +# The component with the least variance, i.e. the axis with the lowest +# number of atoms lying over each other, is aligned to the z-axis, +# which points into the plane of the figure +pca = PCA(n_components=3) +pca.fit(molecule.coord) +molecule = struc.align_vectors(molecule, pca.components_[-1], [0, 0, 1]) # Balls should be colored by partial charge charges = struc.partial_charges(molecule, ITERATION_NUMBER) From b18eb858395988be67296ec96aab432ad558a45f Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 18:15:37 +0100 Subject: [PATCH 09/18] Add proper colorbar --- doc/examples/scripts/structure/peoe_visualization.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 8e50374fe..dcdc73d6f 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -1,6 +1,7 @@ import numpy as np from sklearn.decomposition import PCA import matplotlib.pyplot as plt +from matplotlib.colors import Normalize from matplotlib.cm import ScalarMappable from matplotlib.animation import FuncAnimation import biotite.structure as struc @@ -34,7 +35,8 @@ # Later this variable stores values between 0 and 1 for use in color map charge_in_color_range = charges.copy() # Norm charge values to highest absolute value -charge_in_color_range /= np.max(np.abs(charge_in_color_range)) +max_charge = np.max(np.abs(charge_in_color_range)) +charge_in_color_range /= max_charge # Transform range (-1, 1) to range (0,1) charge_in_color_range = (charge_in_color_range + 1) / 2 # Calculate colors @@ -72,8 +74,11 @@ *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 ) -color_bar = fig.colorbar(ScalarMappable(cmap=color_map)) -color_bar.set_label("Partial charge", color="white") +color_bar = fig.colorbar(ScalarMappable( + norm=Normalize(vmin=-max_charge, vmax=max_charge), + cmap=color_map +)) +color_bar.set_label("Partial charge (e)", color="white") color_bar.ax.yaxis.set_tick_params(color="white") color_bar.outline.set_edgecolor("white") for label in color_bar.ax.get_yticklabels(): From e2d0af953b2e81f943955a560e94f2defb4ff84f Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 25 Nov 2020 18:43:05 +0100 Subject: [PATCH 10/18] Add animation --- .../scripts/structure/peoe_visualization.py | 107 ++++++++++-------- 1 file changed, 59 insertions(+), 48 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index dcdc73d6f..2e15fae8d 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -18,61 +18,66 @@ CMAP_NAME = "bwr_r" -# Get an atom array for the selected molecule -# Caffeine has the PDB reside name 'CFF' -molecule = info.residue(MOLECULE_NAME) +def draw_peoe_plot(fig, molecule, color_map, charges, max_charge): + # Align molecule with principal component analysis: + # The component with the least variance, i.e. the axis with the lowest + # number of atoms lying over each other, is aligned to the z-axis, + # which points into the plane of the figure + pca = PCA(n_components=3) + pca.fit(molecule.coord) + molecule = struc.align_vectors(molecule, pca.components_[-1], [0, 0, 1]) -# Align molecule with principal component analysis: -# The component with the least variance, i.e. the axis with the lowest -# number of atoms lying over each other, is aligned to the z-axis, -# which points into the plane of the figure -pca = PCA(n_components=3) -pca.fit(molecule.coord) -molecule = struc.align_vectors(molecule, pca.components_[-1], [0, 0, 1]) - -# Balls should be colored by partial charge -charges = struc.partial_charges(molecule, ITERATION_NUMBER) -# Later this variable stores values between 0 and 1 for use in color map -charge_in_color_range = charges.copy() -# Norm charge values to highest absolute value -max_charge = np.max(np.abs(charge_in_color_range)) -charge_in_color_range /= max_charge -# Transform range (-1, 1) to range (0,1) -charge_in_color_range = (charge_in_color_range + 1) / 2 -# Calculate colors -color_map = plt.get_cmap(CMAP_NAME) -colors = color_map(charge_in_color_range) + # Balls should be colored by partial charge + # Later this variable stores values between 0 and 1 for use in color map + charge_in_color_range = charges.copy() + # Norm charge values to highest absolute value + charge_in_color_range /= max_charge + # Transform range (-1, 1) to range (0,1) + charge_in_color_range = (charge_in_color_range + 1) / 2 + # Calculate colors + colors = color_map(charge_in_color_range) + + # Ball size should be proportional to VdW radius of the respective atom + ball_sizes = np.array( + [info.vdw_radius_single(e) for e in molecule.element] + ) * BALL_SCALE -# Ball size should be proportional to VdW radius of the respective atom -ball_sizes = np.array( - [info.vdw_radius_single(e) for e in molecule.element] -) * BALL_SCALE + ax = fig.add_subplot(111, projection="3d") + graphics.plot_ball_and_stick_model( + ax, molecule, colors, ball_size=ball_sizes, line_width=3, + line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 + ) -fig = plt.figure(figsize=(8.0, 16.0)) -ax = fig.add_subplot(111, projection="3d") -graphics.plot_ball_and_stick_model( - ax, molecule, colors, ball_size=ball_sizes, line_width=3, - line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 -) + for atom in molecule: + ax.text( + *atom.coord, atom.element, + fontsize=ELEMENT_FONT_SIZE, color="black", + ha="center", va="center", zorder=100 + ) -for atom in molecule: - ax.text( - *atom.coord, atom.element, - fontsize=ELEMENT_FONT_SIZE, color="black", - ha="center", va="center", zorder=100 + # Gradient of ray strength + ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE + ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 + ax.scatter( + *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 + ) + ax.scatter( + *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 ) + + fig.tight_layout() + -# Gradient of ray strength -ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE -ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 -ax.scatter( - *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 -) -ax.scatter( - *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 -) +molecule = info.residue(MOLECULE_NAME) +fig = plt.figure(figsize=(8.0, 8.0)) + +charges = np.stack(struc.partial_charges(molecule, num) for num in range(7)) +max_charge = np.max(np.abs(charges)) +color_map = plt.get_cmap(CMAP_NAME) + +#draw_peoe_plot(fig, molecule, color_map, charges[0], max_charge) color_bar = fig.colorbar(ScalarMappable( norm=Normalize(vmin=-max_charge, vmax=max_charge), @@ -84,5 +89,11 @@ for label in color_bar.ax.get_yticklabels(): label.set_color("white") -fig.tight_layout() +def update(iteration_num): + fig.gca().remove() + draw_peoe_plot( + fig, molecule, color_map, charges[iteration_num], max_charge + ) + +animation = FuncAnimation(fig, update, range(7), interval=1000, repeat=True) plt.show() \ No newline at end of file From e27289efe1252c6f9241d24f1b6f0acb63003d71 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Thu, 26 Nov 2020 09:55:56 +0100 Subject: [PATCH 11/18] Revert "Add animation" This reverts commit e2d0af953b2e81f943955a560e94f2defb4ff84f. --- .../scripts/structure/peoe_visualization.py | 107 ++++++++---------- 1 file changed, 48 insertions(+), 59 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 2e15fae8d..dcdc73d6f 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -18,66 +18,61 @@ CMAP_NAME = "bwr_r" -def draw_peoe_plot(fig, molecule, color_map, charges, max_charge): - # Align molecule with principal component analysis: - # The component with the least variance, i.e. the axis with the lowest - # number of atoms lying over each other, is aligned to the z-axis, - # which points into the plane of the figure - pca = PCA(n_components=3) - pca.fit(molecule.coord) - molecule = struc.align_vectors(molecule, pca.components_[-1], [0, 0, 1]) - - # Balls should be colored by partial charge - # Later this variable stores values between 0 and 1 for use in color map - charge_in_color_range = charges.copy() - # Norm charge values to highest absolute value - charge_in_color_range /= max_charge - # Transform range (-1, 1) to range (0,1) - charge_in_color_range = (charge_in_color_range + 1) / 2 - # Calculate colors - colors = color_map(charge_in_color_range) +# Get an atom array for the selected molecule +# Caffeine has the PDB reside name 'CFF' +molecule = info.residue(MOLECULE_NAME) - # Ball size should be proportional to VdW radius of the respective atom - ball_sizes = np.array( - [info.vdw_radius_single(e) for e in molecule.element] - ) * BALL_SCALE +# Align molecule with principal component analysis: +# The component with the least variance, i.e. the axis with the lowest +# number of atoms lying over each other, is aligned to the z-axis, +# which points into the plane of the figure +pca = PCA(n_components=3) +pca.fit(molecule.coord) +molecule = struc.align_vectors(molecule, pca.components_[-1], [0, 0, 1]) + +# Balls should be colored by partial charge +charges = struc.partial_charges(molecule, ITERATION_NUMBER) +# Later this variable stores values between 0 and 1 for use in color map +charge_in_color_range = charges.copy() +# Norm charge values to highest absolute value +max_charge = np.max(np.abs(charge_in_color_range)) +charge_in_color_range /= max_charge +# Transform range (-1, 1) to range (0,1) +charge_in_color_range = (charge_in_color_range + 1) / 2 +# Calculate colors +color_map = plt.get_cmap(CMAP_NAME) +colors = color_map(charge_in_color_range) +# Ball size should be proportional to VdW radius of the respective atom +ball_sizes = np.array( + [info.vdw_radius_single(e) for e in molecule.element] +) * BALL_SCALE - ax = fig.add_subplot(111, projection="3d") - graphics.plot_ball_and_stick_model( - ax, molecule, colors, ball_size=ball_sizes, line_width=3, - line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 - ) - for atom in molecule: - ax.text( - *atom.coord, atom.element, - fontsize=ELEMENT_FONT_SIZE, color="black", - ha="center", va="center", zorder=100 - ) +fig = plt.figure(figsize=(8.0, 16.0)) +ax = fig.add_subplot(111, projection="3d") +graphics.plot_ball_and_stick_model( + ax, molecule, colors, ball_size=ball_sizes, line_width=3, + line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 +) - # Gradient of ray strength - ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE - ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 - ax.scatter( - *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 - ) - ax.scatter( - *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 +for atom in molecule: + ax.text( + *atom.coord, atom.element, + fontsize=ELEMENT_FONT_SIZE, color="black", + ha="center", va="center", zorder=100 ) - - fig.tight_layout() - -molecule = info.residue(MOLECULE_NAME) -fig = plt.figure(figsize=(8.0, 8.0)) - -charges = np.stack(struc.partial_charges(molecule, num) for num in range(7)) -max_charge = np.max(np.abs(charges)) -color_map = plt.get_cmap(CMAP_NAME) - -#draw_peoe_plot(fig, molecule, color_map, charges[0], max_charge) +# Gradient of ray strength +ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE +ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 +ax.scatter( + *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 +) +ax.scatter( + *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 +) color_bar = fig.colorbar(ScalarMappable( norm=Normalize(vmin=-max_charge, vmax=max_charge), @@ -89,11 +84,5 @@ def draw_peoe_plot(fig, molecule, color_map, charges, max_charge): for label in color_bar.ax.get_yticklabels(): label.set_color("white") -def update(iteration_num): - fig.gca().remove() - draw_peoe_plot( - fig, molecule, color_map, charges[iteration_num], max_charge - ) - -animation = FuncAnimation(fig, update, range(7), interval=1000, repeat=True) +fig.tight_layout() plt.show() \ No newline at end of file From 99638cd792dd91ae6bb016a92b8d8ef092f963ed Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Thu, 26 Nov 2020 10:08:54 +0100 Subject: [PATCH 12/18] Handle NaN charges --- .../scripts/structure/peoe_visualization.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index dcdc73d6f..631d1a41b 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -9,7 +9,7 @@ import biotite.structure.graphics as graphics -MOLECULE_NAME = "CFF" +MOLECULE_NAME = "AIN" ITERATION_NUMBER = 6 ELEMENT_FONT_SIZE = 10 @@ -33,15 +33,16 @@ # Balls should be colored by partial charge charges = struc.partial_charges(molecule, ITERATION_NUMBER) # Later this variable stores values between 0 and 1 for use in color map -charge_in_color_range = charges.copy() +normalized_charges = charges.copy() +normalized_charges[np.isnan(normalized_charges)] = 0 # Norm charge values to highest absolute value -max_charge = np.max(np.abs(charge_in_color_range)) -charge_in_color_range /= max_charge +max_charge = np.max(np.abs(normalized_charges)) +normalized_charges /= max_charge # Transform range (-1, 1) to range (0,1) -charge_in_color_range = (charge_in_color_range + 1) / 2 +normalized_charges = (normalized_charges + 1) / 2 # Calculate colors color_map = plt.get_cmap(CMAP_NAME) -colors = color_map(charge_in_color_range) +colors = color_map(normalized_charges) # Ball size should be proportional to VdW radius of the respective atom ball_sizes = np.array( From 73f687a5b3c25bea52561d88003b233ad1168f35 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Thu, 26 Nov 2020 11:07:34 +0100 Subject: [PATCH 13/18] Add docstring --- .../scripts/structure/peoe_visualization.py | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 631d1a41b..379c8cb46 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -1,9 +1,25 @@ +""" +Partial charge distribution +=========================== + +This examples shows how partial charges are distributed in a +small molecule. +The charges are calculated using the PEOE method [1]_. + +.. [1] J Gasteiger and M Marsili, + "Iterative partial equalization of orbital electronegativity - a + rapid access to atomic charges" + Tetrahedron, 36, 3219 - 3288 (1980). +""" + +# Code source: Patrick Kunzmann +# License: BSD 3 clause + import numpy as np from sklearn.decomposition import PCA import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib.cm import ScalarMappable -from matplotlib.animation import FuncAnimation import biotite.structure as struc import biotite.structure.info as info import biotite.structure.graphics as graphics @@ -13,8 +29,8 @@ ITERATION_NUMBER = 6 ELEMENT_FONT_SIZE = 10 -BALL_SCALE = 200 -RAY_SCALE = 3000 +BALL_SCALE = 500 +RAY_SCALE = 8000 CMAP_NAME = "bwr_r" @@ -51,7 +67,7 @@ -fig = plt.figure(figsize=(8.0, 16.0)) +fig = plt.figure(figsize=(8.0, 6.0)) ax = fig.add_subplot(111, projection="3d") graphics.plot_ball_and_stick_model( ax, molecule, colors, ball_size=ball_sizes, line_width=3, From dc0f01addc3c31a16ac98688b65a5403730fa33b Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Fri, 4 Dec 2020 14:38:24 +0100 Subject: [PATCH 14/18] Finished example script --- .../scripts/structure/peoe_visualization.py | 48 ++++++++++++++----- 1 file changed, 35 insertions(+), 13 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 379c8cb46..04e3802b7 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -25,17 +25,27 @@ import biotite.structure.graphics as graphics +#Acetylsalicylic acid MOLECULE_NAME = "AIN" +# The number of iterations for the PEOE algorithm ITERATION_NUMBER = 6 +# The size of the element lables ELEMENT_FONT_SIZE = 10 -BALL_SCALE = 500 -RAY_SCALE = 8000 +# The scaling factor of the atom 'balls' +BALL_SCALE = 20 +# The higher this number, the more detailed are the rays +N_RAY_STEPS = 20 +# The scaling factor of the 'ray' of charged molecules +RAY_SCALE = 100 +# The transparency value for each 'ray ring' +RAY_ALPHA = 0.03 +# The color map to use to depict the charge CMAP_NAME = "bwr_r" + # Get an atom array for the selected molecule -# Caffeine has the PDB reside name 'CFF' molecule = info.residue(MOLECULE_NAME) # Align molecule with principal component analysis: @@ -50,6 +60,8 @@ charges = struc.partial_charges(molecule, ITERATION_NUMBER) # Later this variable stores values between 0 and 1 for use in color map normalized_charges = charges.copy() +# Show no partial charged for atoms +# that are not parametrized for the PEOE algorithm normalized_charges[np.isnan(normalized_charges)] = 0 # Norm charge values to highest absolute value max_charge = np.max(np.abs(normalized_charges)) @@ -65,15 +77,27 @@ [info.vdw_radius_single(e) for e in molecule.element] ) * BALL_SCALE +# Gradient of ray strength +ray_full_sizes = ball_sizes + np.abs(charges) * RAY_SCALE +ray_sizes = np.array([ + np.linspace(ray_full_sizes[i], ball_sizes[i], N_RAY_STEPS, endpoint=False) + for i in range(molecule.array_length()) +]).T +# The plotting begins here fig = plt.figure(figsize=(8.0, 6.0)) ax = fig.add_subplot(111, projection="3d") + +# Plot the atoms +# As `axes.scatter` uses sizes in points**2, +# the VdW-radii as also squared graphics.plot_ball_and_stick_model( - ax, molecule, colors, ball_size=ball_sizes, line_width=3, + ax, molecule, colors, ball_size=ball_sizes**2, line_width=3, line_color=color_map(0.5), background_color=(.05, .05, .05), zoom=1.5 ) +# Plot the element labels for atom in molecule: ax.text( *atom.coord, atom.element, @@ -81,16 +105,14 @@ ha="center", va="center", zorder=100 ) -# Gradient of ray strength -ray_full_size = ball_sizes + np.abs(charges) * RAY_SCALE -ray_half_size = ball_sizes + np.abs(charges) * RAY_SCALE/2 -ax.scatter( - *molecule.coord.T, s=ray_full_size, c=colors, linewidth=0, alpha=0.2 -) -ax.scatter( - *molecule.coord.T, s=ray_half_size, c=colors, linewidth=0, alpha=0.2 -) +# Plots the rays +for i in range(N_RAY_STEPS): + ax.scatter( + *molecule.coord.T, s=ray_sizes[i]**2, c=colors, + linewidth=0, alpha=RAY_ALPHA + ) +# Plot the colorbar color_bar = fig.colorbar(ScalarMappable( norm=Normalize(vmin=-max_charge, vmax=max_charge), cmap=color_map From 9ba4e93b36923d562c33c279899d36e421623361 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Fri, 4 Dec 2020 14:38:34 +0100 Subject: [PATCH 15/18] Fix docstring --- src/biotite/structure/charges.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biotite/structure/charges.py b/src/biotite/structure/charges.py index f8fe7e9e7..a6ff50ccb 100644 --- a/src/biotite/structure/charges.py +++ b/src/biotite/structure/charges.py @@ -210,7 +210,7 @@ def partial_charges(atom_array, iteration_step_num=6, charges=None): Otherwise, an error will be raised. Example: - .. highlight:: python + .. code-block:: python atom_array.bonds = struc.connect_via_residue_names(atom_array) From a83f713faa1df0d56acafbdfed5fcd8ba5a7e127 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Fri, 4 Dec 2020 14:48:16 +0100 Subject: [PATCH 16/18] Docstring addition --- doc/examples/scripts/structure/peoe_visualization.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 04e3802b7..769ab3a96 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -25,7 +25,7 @@ import biotite.structure.graphics as graphics -#Acetylsalicylic acid +# Acetylsalicylic acid MOLECULE_NAME = "AIN" # The number of iterations for the PEOE algorithm @@ -78,6 +78,7 @@ ) * BALL_SCALE # Gradient of ray strength +# The ray size is proportional to the absolute charge value ray_full_sizes = ball_sizes + np.abs(charges) * RAY_SCALE ray_sizes = np.array([ np.linspace(ray_full_sizes[i], ball_sizes[i], N_RAY_STEPS, endpoint=False) From ac4283ec2b3e79785290882da47b90e191a02d74 Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Fri, 4 Dec 2020 15:10:22 +0100 Subject: [PATCH 17/18] Add docstring --- src/biotite/structure/graphics/atoms.py | 67 +++++++++++++++++++++---- 1 file changed, 58 insertions(+), 9 deletions(-) diff --git a/src/biotite/structure/graphics/atoms.py b/src/biotite/structure/graphics/atoms.py index bf790ea32..4e7bcc8b1 100644 --- a/src/biotite/structure/graphics/atoms.py +++ b/src/biotite/structure/graphics/atoms.py @@ -23,22 +23,18 @@ def plot_atoms(axes, atoms, colors, line_width=1.0, background_color=None, ---------- axes : Axes3D The *Matplotlib* 3D-axes to plot the structure on. - For a correct aspect ratio of the visualized structure, the - :class:`Axes3D` must have quadratic extents on the display. - This can be approximately achieved by giving the :class:`Figure` - equal values for the width and height in its `figsize`. - atoms : AtomArray, length=n + atoms : AtomArray, shape=(n,) The structure to be plotted. The atom array must have an associated :class:`BondList`, i.e. - its `bonds` attribute must be defined. + its ``bonds`` attribute must be defined. One line for each bond is drawn. colors : ndarray, shape=(n,3) or shape=(n,4), dtype=float An array of RGB or RGBA colors for each atom in `atoms`. The values for each color channel are in the range 0 to 1. line_width : float, optional The width of the lines to be drawn. - background_color : string or iterable object - A matplotlib compatible color (color name or RGB values). + background_color : object + A *Matplotlib* compatible color (color name or RGB values). If set, the background is colored with the given value. center : tuple of (float, float, float), optional The coordinates of the structure that are in the center of the @@ -112,6 +108,59 @@ def plot_ball_and_stick_model(axes, atoms, colors, ball_size=200, line_color="black", line_width=1.0, background_color=None, center=None, size=None, zoom=1.0): + """ + Plot an :class:`AtomArray` as *ball-and-stick* model. + + The z-axis points into the screen plane. + + UNSTABLE: This function is probably subject to future changes. + + Parameters + ---------- + axes : Axes3D + The *Matplotlib* 3D-axes to plot the structure on. + atoms : AtomArray, length=n + The structure to be plotted. + The atom array must have an associated :class:`BondList`, i.e. + its ``bonds`` attribute must be defined. + One line for each bond is drawn. + colors : ndarray, shape=(n,3) or shape=(n,4), dtype=float + An array of RGB or RGBA colors for each atom in `atoms`, that + is used to color the *balls* of the model. + The values for each color channel are in the range 0 to 1. + ball_size : int or iterable of int, shape=(n,) + The size of the *balls* in the model in *pt*. + Either a single value for all atoms or an iterable object of + values for each atom. + line_color : object + A *Matplotlib* compatible color value for the *sticks* of the + model. + line_width : float, optional + The width of the *sticks* in the model in *pt*. + background_color : object + A *Matplotlib* compatible color (color name or RGB values). + If set, the background is colored with the given value. + center : tuple of (float, float, float), optional + The coordinates of the structure that are in the center of the + plot. + By default the complete molecule is centered. + size : float, optional + The size of each dimension in the plot. + The limits of the :class:`Axes3D` are set to + ``(center - size/(2*zoom)), (center + size/(2*zoom))``. + zoom : float, optional + Zoom in or zoom out. + + - ``> 1.0``: Zoom in. + - ``< 1.0``: Zoom out. + + Notes + ----- + This is a very simple visualization tools for quick visual analysis + of a structure. + For publication-ready molecular images the usage of a dedicated + molecular visualization tool is recommended. + """ if not isinstance(axes, Axes3D): raise ValueError("The given axes mut be an 'Axes3D'") if atoms.bonds is None: @@ -169,7 +218,7 @@ def _set_box(axes, coord, center, size, zoom): axes.set_xlim(center[0] - size/(2*zoom), center[0] + size/(2*zoom)) axes.set_ylim(center[1] - size/(2*zoom), center[1] + size/(2*zoom)) axes.set_zlim(center[2] - size/(2*zoom), center[2] + size/(2*zoom)) + # Make the axis lengths of the 'plot box' equal # The 'plot box' is not visible due to 'axes.axis("off")' - axes.set_box_aspect([1,1,1]) \ No newline at end of file From 7dc6e52cc03b85f503cec621691144dd3f742d6f Mon Sep 17 00:00:00 2001 From: Patrick Kunzmann Date: Wed, 9 Dec 2020 12:57:17 +0100 Subject: [PATCH 18/18] Fix comments --- doc/examples/scripts/structure/peoe_visualization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/scripts/structure/peoe_visualization.py b/doc/examples/scripts/structure/peoe_visualization.py index 769ab3a96..2ef238acb 100644 --- a/doc/examples/scripts/structure/peoe_visualization.py +++ b/doc/examples/scripts/structure/peoe_visualization.py @@ -60,13 +60,13 @@ charges = struc.partial_charges(molecule, ITERATION_NUMBER) # Later this variable stores values between 0 and 1 for use in color map normalized_charges = charges.copy() -# Show no partial charged for atoms +# Show no partial charge for atoms # that are not parametrized for the PEOE algorithm normalized_charges[np.isnan(normalized_charges)] = 0 # Norm charge values to highest absolute value max_charge = np.max(np.abs(normalized_charges)) normalized_charges /= max_charge -# Transform range (-1, 1) to range (0,1) +# Transform range (-1, 1) to range (0, 1) normalized_charges = (normalized_charges + 1) / 2 # Calculate colors color_map = plt.get_cmap(CMAP_NAME) @@ -91,7 +91,7 @@ ax = fig.add_subplot(111, projection="3d") # Plot the atoms -# As `axes.scatter` uses sizes in points**2, +# As 'axes.scatter()' uses sizes in points**2, # the VdW-radii as also squared graphics.plot_ball_and_stick_model( ax, molecule, colors, ball_size=ball_sizes**2, line_width=3,