[1]:
import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from membrane_curvature.tests.datafiles import (MEMB_GRO,
MEMB_XTC)
u = mda.Universe(MEMB_GRO, MEMB_XTC)
mc_upper = MembraneCurvature(u,
select="resid 103-1023 and name PO4",
n_x_bins=12,
n_y_bins=12).run()
mc_lower = MembraneCurvature(u,
select='resid 1024-2046 and name PO4',
n_x_bins=12,
n_y_bins=12).run()
MDAnalysis : INFO MDAnalysis 2.0.0-dev0 STARTED logging to 'MDAnalysis.log'
/Users/estefania/python_playground/mdanalysis/package/MDAnalysis/topology/guessers.py:146: UserWarning: Failed to guess the mass for the following atom types: D
warnings.warn("Failed to guess the mass for the following atom types: {}".format(atom_type))
/Users/estefania/python_playground/mdanalysis/package/MDAnalysis/topology/guessers.py:146: UserWarning: Failed to guess the mass for the following atom types: G
warnings.warn("Failed to guess the mass for the following atom types: {}".format(atom_type))
/Users/estefania/python_playground/mdanalysis/package/MDAnalysis/topology/guessers.py:146: UserWarning: Failed to guess the mass for the following atom types: R
warnings.warn("Failed to guess the mass for the following atom types: {}".format(atom_type))
MDAnalysis.analysis.base: INFO Choosing frames to analyze
MDAnalysis.analysis.base: INFO Starting preparation
/Users/estefania/python_playground/membrane-curvature/membrane_curvature/surface.py:101: UserWarning: Atom outside grid boundaries. Skipping atom.
warnings.warn(msg)
MDAnalysis.analysis.base: INFO Finishing up
MDAnalysis.analysis.base: INFO Choosing frames to analyze
MDAnalysis.analysis.base: INFO Starting preparation
MDAnalysis.analysis.base: INFO Finishing up
[2]:
# surface
surf_upper = mc_upper.results.average_z_surface
surf_lower = mc_lower.results.average_z_surface
# mean curvature
mean_upper = mc_upper.results.average_mean
mean_lower = mc_lower.results.average_mean
# Gaussian curvature
gauss_upper = mc_upper.results.average_gaussian
gauss_lower = mc_lower.results.average_gaussian
[4]:
import numpy as np
print(np.max(mean_lower), np.max(mean_upper))
0.9152833670901555 0.8257232286053238
[10]:
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
results = [mean_lower*-1, mean_upper*-1]
leaflets = ["Lower", "Upper"]
fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4,3.5), dpi=200)
max_ = max([np.max(abs(h)) for h in results])
for ax, rs, lf in zip((ax1, ax2), results, leaflets):
rs = ndimage.zoom(rs, 4, mode='wrap')
levs = np.linspace(-max_, max_, 40)
im = ax.contourf(rs, cmap='bwr',
origin='lower', levels=levs,
vmin=-max_, vmax=max_)
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[-max_, 0, max_],
orientation='vertical',
ax=[ax1, ax2], shrink=0.4,
aspect=10, pad=0.05)
cbar.ax.tick_params(labelsize=4, width=0.5)
cbar.set_label('H (nm$^{-1}$)', fontsize=6, labelpad=2)
plt.savefig("/Users/estefania/python_playground/ojeda-e.github.io/assets/images/my_contours_dark.png", dpi=300)
[46]:
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
results = [mean_lower*-1, mean_upper*-1]
leaflets = ["Lower", "Upper"]
fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4,3.5), dpi=200)
max_ = np.max(results)
[print(max_)]
for ax, rs, lf in zip((ax1, ax2), results, leaflets):
rs = ndimage.zoom(rs, 5, mode='wrap')
levs = np.linspace(-max_, max_, 40)
im = ax.contourf(rs, cmap='bwr',
origin='lower', levels=levs,
vmin=-max_, vmax=max_)
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[-max_, 0, max_], orientation='vertical', ax=[ax1, ax2], shrink=0.4, aspect=10, pad=0.05)
cbar.ax.tick_params(labelsize=4, width=0.5)
cbar.set_label('H (nm$^{-1}$)', fontsize=6, labelpad=2)
plt.savefig("/Users/estefania/python_playground/ojeda-e.github.io/assets/images/my_contours_dark.png", dpi=300)
0.5751326866434266
[12]:
cms=["YlGnBu_r", "bwr", "PiYG"]
units=['$Z (\AA)$','$H$(nm$^{-1})$', '$K$(nm$^{-2})$']
results = [surf_upper, mean_upper*-1, gauss_upper]
titles = ['Surface', 'Mean Curvature', 'Gaussian Curvature']
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(7,4), dpi=200)
for ax, mc, title, cm, unit in zip((ax1, ax2, ax3), results, titles, cms, units):
mc = ndimage.zoom(mc, 5, mode='wrap')
bound = max(abs(np.min(mc)), abs(np.max(mc)))
print(bound)
if np.min(mc) < 0 < np.max(mc):
levs = np.linspace(-bound, +bound, 40)
im = ax.contourf(mc, cmap=cm, levels=levs, alpha=0.95, vmin=-bound, vmax=+bound)
tcs = [-bound, 0, bound]
else:
im = ax.contourf(mc, cmap=cm, levels=40, alpha=0.95, vmin=int(np.min(mc)), vmax=np.floor(np.max(mc)))
tcs = [int(np.min(mc)), np.floor(np.max(mc))]
ax.set_aspect('equal')
ax.set_title(title, fontsize=12)
ax.axis('off')
cbar=plt.colorbar(im, ticks=[], ax=ax, orientation='horizontal', pad=0.05)
cbar.ax.tick_params(labelsize=7, width=0.5)
cbar.set_label(unit, fontsize=9, labelpad=2)
plt.tight_layout()
plt.savefig("/Users/estefania/python_playground/ojeda-e.github.io/assets/images/my_results_dark.png", transparent=True, dpi=300)
144.2566415631751
0.8139704506587123
0.765111399161304
[9]:
surf_upper = curv_upper_leaflet.results.average_z_surface
mean_upper = curv_upper_leaflet.results.average_mean
gauss_upper = curv_upper_leaflet.results.average_gaussian
[10]:
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
mean_lower = curv_upper_leaflet.results.average_mean
curvatures = [mean_lower, mean_upper]
leaflets = ['Lower', 'Upper']
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(3,2), dpi=200)
for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
mc = ndimage.zoom(mc, 5, mode='wrap')
bound = max(abs(np.min(mc)), abs(np.max(mc)))
im = ax.contourf(mc, cmap='bwr', levels=40, alpha=0.95, vmin=-bound, vmax=+bound)
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=10)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[np.min(mc), 0, np.max(mc)], orientation='vertical',
ax=[ax1, ax2], shrink=0.6, aspect=10, pad=0.01)
cbar.ax.tick_params(labelsize=5, width=0.5)
cbar.set_label('$H$(nm$^{-1})$', fontsize=6, labelpad=2)
[6]:
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
white='white'
plt.rcParams['text.color'] = white
plt.rcParams['axes.labelcolor'] = white
plt.rcParams['xtick.color'] = white
plt.rcParams['ytick.color'] = white
cms=["YlGnBu", "bwr", "PiYG"]
units=['$Z (\AA)$','$H$(nm$^{-1})$', '$K$(nm$^{-2})$']
results = [surf_upper, mean_upper, gauss_upper]
titles = ['Surface', 'Mean Curvature', 'Gaussian Curvature']
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(7,4), dpi=200)
for ax, mc, title, cm, unit in zip((ax1, ax2, ax3), results, titles, cms, units):
mc = ndimage.zoom(mc, 5, mode='wrap')
bound = max(abs(np.min(mc)), abs(np.max(mc)))
if np.min(mc) < 0 < np.max(mc):
im = ax.contourf(mc, cmap=cm, levels=40, alpha=0.95, vmin=-bound, vmax=+bound)
else:
im = ax.contourf(mc, cmap=cm, levels=40, alpha=0.95)
ax.set_aspect('equal')
ax.set_title(title, fontsize=12)
ax.axis('off')
cbar=plt.colorbar(im, ticks=[np.min(mc), 0, np.max(mc)] if np.min(mc) < 0 < np.max(mc) else [np.min(mc), np.max(mc)], ax=ax, orientation='horizontal', pad=0.05)
cbar.ax.tick_params(labelsize=7, width=0.5)
cbar.set_label(unit, fontsize=9, labelpad=2)
plt.tight_layout()
plt.savefig("/Users/estefania/python_playground/ojeda-e.github.io/assets/images/results.png", transparent=True, dpi=300)
[9]:
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np
white='white'
plt.rcParams['text.color'] = white
plt.rcParams['axes.labelcolor'] = white
plt.rcParams['xtick.color'] = white
plt.rcParams['ytick.color'] = white
curvatures = [mean_lower, mean_upper]
leaflets = ['Lower', 'Upper']
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(3,2), dpi=200)
for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
mc = ndimage.zoom(mc, 5, mode='wrap')
bound = max(abs(np.min(mc)), abs(np.max(mc)))
im = ax.contourf(mc, cmap='bwr', levels=40, alpha=0.95, vmin=-bound, vmax=+bound)
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[np.min(mc), 0, np.max(mc)], orientation='vertical', ax=[ax1, ax2], shrink=0.5, aspect=10, pad=0.05)
cbar.ax.tick_params(labelsize=4, width=0.5)
cbar.set_label('$H$(nm$^{-1})$', fontsize=5, labelpad=2)
#plt.savefig("/Users/estefania/python_playground/ojeda-e.github.io/assets/images/mycontours.png", transparent=True, dpi=300)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-9-0b95621b571f> in <module>
10
11
---> 12 curvatures = [mean_lower, mean_upper]
13 leaflets = ['Lower', 'Upper']
14
NameError: name 'mean_lower' is not defined
print(np.min(mc))¶
[ ]:
[ ]:
[23]:
import MDAnalysis as mda
import nglview as nv
from membrane_curvature.base import MembraneCurvature
from membrane_curvature.tests.datafiles import MEMB_GRO, MEMB_XTC
from nglview.contrib.movie import MovieMaker
universe = mda.Universe(MEMB_GRO, MEMB_XTC)
view = nv.show_mdanalysis(universe)
view.add_representation('ball+stick', selection='not CHOL', radius=0.5, color="lightgrey")
view.add_representation('ball+stick', selection='.PO4', radius=2.00, color="darkcyan")
view.add_representation('ball+stick', selection='not POPC and not POPE', radius=0.5, color="lightgreen")
view.add_representation('ball+stick', selection='.ROH', radius=2.00, color="lightgreen")
view.camera='orthographic'
view.control.zoom(0.3)
view.control.rotate(
mda.lib.transformations.quaternion_from_euler(
-np.pi/2, np.pi/3, np.pi/12, 'rzyz').tolist())
view
[14]:
movie = MovieMaker(view, start=0, stop=5, fps=1, output='membrane_movie.gif')
movie.make()
[41]:
import pytraj as pt
import mdtraj as md
traj = pt.load(MEMB_XTC, top=md.load(MEMB_GRO).topology)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-41-5066ca6bc1d0> in <module>
2 import mdtraj as md
3
----> 4 traj = pt.load(MEMB_XTC, top=md.load(MEMB_GRO).topology)
~/opt/anaconda3/envs/membcurv-dev/lib/python3.7/site-packages/pytraj/io.py in load(filename, top, frame_indices, mask, stride)
122 # load to TrajectoryIterator object first
123 # do not use frame_indices_ here so we can optimize the slicing speed
--> 124 traj = load_traj(filename, top, stride=stride)
125
126 # do the slicing and other things if needed.
~/opt/anaconda3/envs/membcurv-dev/lib/python3.7/site-packages/pytraj/io.py in load_traj(filename, top, *args, **kwd)
239 if isinstance(top, string_types):
240 top = load_topology(top)
--> 241 if top is None or top.is_empty():
242 top = load_topology(filename)
243 ts = TrajectoryIterator(top=top)
AttributeError: 'Topology' object has no attribute 'is_empty'
[27]:
traj.superpose(mask='PO4', ref=0)
traj
[27]:
pytraj.Trajectory, 51 frames:
Size: 0.006323 (GB)
<Topology: 5547 atoms, 349 residues, 1 mols, PBC with box type = ortho>
[28]:
view = nv.show_pytraj(traj)
view
[39]:
topology = md.load(MEMB_GRO).topology
[ ]: