Giulio Tesei, Erik Hellstrand, Sara Linse, Emma Sparr, Robert Vacha, and Mikael Lund
This Jupyter Notebook was originally run in MacOS 10.11 with GCC 4.8, cmake, Python3, matplotlib, pandas within the Anaconda environment. Contemporary Linux distributions such as Ubuntu 14.04 will also work.
In [1]:
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np, pandas as pd
import os.path, os, sys, json, shutil
from pathlib import Path
try:
workdir
except NameError:
workdir=%pwd
else:
%cd $workdir
To simulate the interaction of rigid fibrils with a charged, planar surface we use a custom MC program build within the Faunus framework. The sections below will fetch the Faunus project and compile the appropriate version of the program.
In [2]:
%%bash -s "$workdir"
cd $1
if [ ! -d "faunus" ]; then
git clone https://github.com/mlund/faunus.git
cd faunus
git checkout dfdddf3
else
cd faunus
fi
cmake . -DCMAKE_BUILD_TYPE=Release -DENABLE_APPROXMATH=on -DMYPLAYGROUND=`pwd`/../mc/ &>/dev/null
make r2s -j4
cd $1
In [2]:
%%writefile mc/submit.sh
#!/bin/bash
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 01:00:00
../r2s > out
This will generate all input for MC simulations as well as loop over the full parameter space. In the inner loop you may specify if the jobs should be run locally (not recommended) or be submitted to a cluster, assuming you are running this Notebook from a cluster environment.
In [3]:
%cd $workdir
js = {
"atomlist" : {
"Na" : { "q": 1, "r":1.9, "mw":22.99 },
"Cl" : { "q":-1, "r":1.7, "mw":35.45 },
"I" : { "q":-1, "r":2.0, "mw":1 },
"SCN" : { "q":-1, "r":2.0, "mw":1 },
"ASP" : { "q":-1, "r":3.6, "mw":110 },
"HASP" : { "q":0, "r":3.6, "mw":110 },
"LASP" : { "q":2, "r":3.6, "mw":110 },
"CTR" : { "q":-1, "r":2.0, "mw":16 },
"HCTR" : { "q":0, "r":2.0, "mw":16 },
"GLU" : { "q":-1, "r":3.8, "mw":122 },
"HGLU" : { "q":0, "r":3.8, "mw":122 },
"LGLU" : { "q":2, "r":3.8, "mw":122 },
"HIS" : { "q":0, "r":3.9, "mw":130 },
"HHIS" : { "q":1, "r":3.9, "mw":130 },
"NTR" : { "q":0, "r":2.0, "mw":14 },
"HNTR" : { "q":1, "r":2.0, "mw":14 },
"TYR" : { "q":-1, "r":4.1, "mw":154 },
"HTYR" : { "q":0, "r":4.1, "mw":154 },
"LYS" : { "q":0, "r":3.7, "mw":116 },
"HLYS" : { "q":1, "r":3.7, "mw":116 },
"CYS" : { "q":-1, "r":3.6, "mw":103 },
"HCYS" : { "q":0, "r":3.6, "mw":103 },
"ARG" : { "q":0, "r":4.0, "mw":144 },
"HARG" : { "q":1, "r":4.0, "mw":144 },
"ALA" : { "q":0, "r":3.1, "mw":66 },
"ILE" : { "q":0, "r":3.6, "mw":102 },
"LEU" : { "q":0, "r":3.6, "mw":102 },
"MET" : { "q":0, "r":3.8, "mw":122 },
"PHE" : { "q":0, "r":3.9, "mw":138 },
"PRO" : { "q":0, "r":3.4, "mw":90 },
"TRP" : { "q":0, "r":4.3, "mw":176 },
"VAL" : { "q":0, "r":3.4, "mw":90 },
"SER" : { "q":0, "r":3.3, "mw":82 },
"THR" : { "q":0, "r":3.5, "mw":94 },
"ASN" : { "q":0, "r":3.6, "mw":108 },
"GLN" : { "q":0, "r":3.8, "mw":120 },
"GLY" : { "q":0, "r":2.9, "mw":54 }
}
}
with open('mc/slitpolymer.json', 'w+') as f:
f.write(json.dumps(js, indent=4))
def mkinput():
d = {
'atomlist': '../slitpolymer.json',
'loop_macrosteps': 10,
'loop_microsteps': micro,
'temperature': 298, # Kelvin
'epsilon_r': 78.9, # Water dielectric const
'dh_ionicstrength': salt, # mol/l
'gouychapman_qarea': rho, # unitless i.e. e*phi/kT
'gouychapman_linearize': 'yes',
'gouychapman_offset': 3,
'cuboid_len': -1, # Box side length Angstrom
'cuboid_xlen': 2000,
'cuboid_ylen': 2000,
'cuboid_zlen': 500,
'transrot_transdp': 100, # Molecular translation parameter
'transrot_rotdp': 6, # Molecular rotation parameter
'oligomer': '../orderedrod.aam',
'numpep': numpep
}
f = open('r2s.input', 'w')
for key,val in d.iteritems():
f.write('%-30s %s\n' %(key, val))
f.close()
%cd -q $workdir'/mc'
# parameter space
rho_range = [266] # surface charge density (angstrom^2/charge)
salt_range = [0.06, 0.08, 0.1, 0.15, 0.2, 0.25, 0.27, 0.3, 1e10] # ionic strength (mol/l)
numpep_range = [2, 4, 8, 16, 32, 64, 128] # fibril length
# run through all parameter space
micro=20000000
if not Path('monte-carlo-data.p').is_file(): # run only if this file is missing (see next section)
for rho in rho_range:
for salt in salt_range:
for numpep in numpep_range:
pfx='__rho'+str(rho)+'-N'+str(numpep)+'-I'+str(salt)
if not os.path.isdir(pfx):
%mkdir -p $pfx
%cd $pfx
mkinput()
if shutil.which('sbatch') is not None:
!sbatch ../submit.sh # run on slurm cluster
else:
!../r2s > out # run locally
%cd -q ..
We now load all distribution functions into a Pandas object for further analysis. We also integrate the potential of mean force, $w(r)$ to a surface excess,
$$ \Gamma / \rho = \int_{contact}^{\infty} (e^{-w(r)/k_BT} - 1) dr $$Note: All loaded MC data will be saved to disk as a single pickle file, monte-carlo-data.p. The next time this cell is run, this file will be loaded and no raw MC data is needed.
In [4]:
%cd -q $workdir'/mc'
if Path('monte-carlo-data.p').is_file():
data = pd.read_pickle('monte-carlo-data.p')
else:
d = {}
rho = rho_range[0]
for salt in salt_range:
d[salt] = {}
for numpep in numpep_range:
pfx='__rho'+str(rho)+'-N'+str(numpep)+'-I'+str(salt)
distfile=pfx+'/cmdist.dat'
if os.path.isfile(distfile):
r, P = np.loadtxt(distfile, unpack=True) # load data set
w = -np.log(P) # potential of mean force (PMF)
# normalize such that w(r)->0 for large r
m = (r>200) & (r<220)
if len(w[m])>0:
w = w - w[m].mean()
# remove tail
w = w[r<200]
r = r[r<200]
# integrate PMF to surface excess
gamma = np.trapz( np.expm1(-w), r )
d[salt][numpep] = { 'r': r, 'pmf': w, 'gamma': gamma }
if len( d[salt] ) > 0:
d[salt] = pd.Series(d[salt])
data = pd.DataFrame(d)
data = data.rename(columns={1e10: r'$\infty$'})
data.to_pickle('monte-carlo-data.p')
data
Out[4]:
In [5]:
plt.rcParams.update({'font.size': 16, 'figure.figsize': [6.0, 5.0]})
for N, d in data[0.25].iteritems(): # at 0.25 mol/l salt
if N in [2,8,16,32,64,128]:
if type(d)!=type(float()):
plt.plot( d['r']*0.1, d['pmf'], label='$N$='+str(N), lw=2)
plt.legend(loc=0, frameon=False, fontsize='medium')
plt.xlabel(r'Surface-Fibril Separation, $r$ (nm)')
plt.ylabel(r'Free energy, $\beta w(r)$')
plt.ylim((-4,3))
plt.xlim((0,18))
plt.text(15,1.7,'b',fontsize=40)
plt.savefig('mc_pmf.pdf', bbox_inches='tight')
In [6]:
for salt, d in data.items():
if str(salt) in ['0.1', '0.15', '0.2', '0.25', '0.3', r'$\infty$']:
N = d.keys()
gamma = np.array( [ i['gamma'] for i in d.get_values() if type(i)==dict] )
plt.plot(N, gamma*0.1, label=salt, linestyle='-', marker='o', lw=2, ms=6)
plt.plot( [0,128], [0,0], 'k--')
plt.legend(loc='lower right',frameon=False, fontsize='medium', title=r'$c_{s}$ (M)')
plt.xlabel(r'Number of monomers, $N$')
plt.ylabel(r'Surface excess, $\Gamma/\rho$ (nm)')
plt.xscale('log')
plt.ylim((-8,10))
plt.xlim(1, 1000)
plt.text(300,7,'a',fontsize=40)
print(os.getcwd())
plt.savefig('mc_excess.pdf', bbox_inches='tight')
Here we calculate the surface excess using a simplified model, considering Gouy-Chapman electrostatics and the excluded volume of a rod approaching a planar surface. As detailed in the main text, the free energy can be approximated with $$ \beta w(r)^{\scriptsize{\mbox{rod}}} \approx - 2 z \ln \left[ \frac{1+\tanh (\beta \phi e /4)e^{-\frac{r+s}{\lambda_D}}}{1-\tanh (\beta \phi e /4)e^{-\frac{r+s}{\lambda_D}}} \right] - \ln(2r/L) \vee 0
In [7]:
%cd -q $workdir
from math import asinh, sqrt, pi
s=10 # closest contact (Å)
lm=5 # monomer length (Å)
zm=2 # monomer charge (e)
lB=7.1 # Bjerrum length (Å)
NA=6.022e23 # Avogadros number (1/mol)
def TdS(r, L):
""" rod-2-surface excluded volume contribution to free energy (kT) """
sel = r<0.5*L
w = -np.log(2*r[sel]/L)
return np.concatenate((w, 0*r[~sel]))
for rho in [266]:
for N in [2,8,16,32,64,128]: # number of monomers
for salt in [0.1]: # ionic strength (mol/l)
L = N*lm # rod length (Å)
z = zm*N # rod charge (e)
lD = 3.04/sqrt(salt) # Debye length (Å)
phi = 2*asinh( (1/rho)*sqrt(pi*lB/(2*salt*NA*1e-27)) ) # surf. potential (kT/e)
r = np.arange(s, 400, 0.5)
#w = -z * phi * np.exp(-(r+s)/lD) + TdS(r,L)
w = - 2 * z * np.log( ( 1 + np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) )
/ ( 1 - np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) ) ) + TdS(r,L)+ (s/r)**9
gamma = np.trapz( np.expm1(-w), r )
plt.plot(r*0.1, w, label=r'$N=$'+str(N), lw=2)
plt.legend(loc=0,frameon=False, fontsize='medium')
plt.xlabel('Surface-Fibril Separation, $r$ (nm)')
plt.ylabel(r'Free energy, $\beta w(r)$')
plt.xlim((0,35))
plt.ylim((-4,2))
plt.show()
In [8]:
from pylab import cm, imshow, contour, clabel, colorbar
x = np.linspace(0.05,0.45, 221) # salt concentration
y = np.linspace(1, 100, 160) # number of monomers
Z = np.ndarray(shape=(y.size,x.size))
rho = 266
for N in y: # number of monomers
for salt in x: # ionic strength (mol/l)
L = N*lm # rod length (Å)
z = zm*N # rod charge (e)
lD = 3.04/sqrt(salt) # Debye length (Å)
phi = 2*asinh( (1/rho)*sqrt(pi*lB/(2*salt*NA*1e-27)) ) # surf. potential (kT/e)
r = np.arange(0, 1000, 0.2)
# eq 62. DOI: http://dx.doi.org/10.1016/S0006-3495(02)75375-7
w = - 2 * z * np.log( ( 1 + np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) )
/ ( 1 - np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) ) ) + TdS(r,L)+ (s/r)**9
#w = -z * phi * np.exp(-(r+s)/lD) + TdS(r,L)+ (s/r)**9
exc = np.trapz( np.expm1(-w), r )*0.1
Z[y.tolist().index(N)][x.tolist().index(salt)]=exc
In [238]:
vmin=-10
vmax=10
print(Z.min())
im = imshow(Z,cmap=cm.RdBu,extent=[x.min(),x.max(),y.min(),y.max()],origin='lower',aspect='auto', vmin=vmin, vmax=vmax)
# contour lines
cset = contour(x, y, Z, (-1,0), linewidths=4, cmap=cm.RdBu)
manual_locations = [(0.3, 10),(0.18, 10)]
clabel(cset, inline=True, fmt='%1.f', fontsize=15, manual=manual_locations)
plt.figtext(0.3,0.8,'Attraction',rotation='66', color='white',fontsize='x-large', style='italic')
plt.figtext(0.63,0.8,'Repulsion',rotation='82',color='white',fontsize='x-large', style='italic')
# labels
colorbar(im,ticks=np.arange(-10,12,2),label=r'Surface Excess, $\Gamma/\rho$ (nm)') # adding the colorbar on the right
plt.ylabel('Number of Monomers, $N$')
plt.xlabel('Salt Concentration, $c_s$ (M)')
plt.xticks(np.arange(.05,.50,.1))
plt.text(0.077,83,'b',fontsize=40,color='w')
plt.savefig('map.png')
plt.show()
In [9]:
y = np.linspace(1, 100, 160) # number of monomers
salt = [0.1, 0.15, 0.2, 0.25, 0.3, 1]
for cs in salt:
salt = 0.25
Y = []
X = []
for N in y: # number of monomers
L = N*lm # rod length (Å)
z = zm*N # rod charge (e)
lD = 3.04/sqrt(cs) # Debye length (Å)
phi = 2*asinh( (1/rho)*sqrt(pi*lB/(2*salt*NA*1e-27)) ) # surf. potential (kT/e)
r = np.arange(0, 1000, 0.2)
#w = - 2 * z * np.log( ( 1 + np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) )
# / ( 1 - np.tanh( phi / 4 ) * np.exp(-(r+s)/lD) ) ) + TdS(r,L) + (s/r)**9
w = -z * phi * np.exp(-(r+s)/lD) + TdS(r,L) + (s/r)**9
exc = np.trapz( np.expm1(-w), r )
Y.append(N)
X.append(exc*0.1)
plt.plot(Y,X, label=cs, linestyle='-', marker='o', lw=2, ms=6, markevery = 8)
plt.legend(loc='lower right',frameon=False, fontsize='medium', title=r'$c_{s}$ (M)')
plt.ylim(-8,10)
plt.xlim(1,1000)
plt.xlabel(r'Number of monomers, $N$')
plt.ylabel(r'Surface excess, $\Gamma/\rho$ (nm)')
plt.xscale('log')
plt.text(1.5,7,'a',fontsize=40)
plt.savefig('linesegment_excess.pdf')
plt.show()
In [14]:
F50fib = np.loadtxt('exp/50fib_1.dat',comments={'T'},usecols=(0,1,6,9))
D50fib = np.loadtxt('exp/50fib_1.dat',comments={'T'},usecols=(0,3,7,11))
Fit50fib = np.loadtxt('exp/50fib_1.dat',comments={'T'},usecols=(0,2,4,6,8,10,12))
shear50fib = np.loadtxt('exp/par50fib_1.dat',comments={'t'},usecols=(0,1,3,5))
mass50fib = np.loadtxt('exp/par50fib_1.dat',comments={'t'},usecols=(0,2,4,6))
F50fib2 = np.loadtxt('exp/50fib_2.dat',comments={'T'},usecols=(0,1,6,9))
D50fib2 = np.loadtxt('exp/50fib_2.dat',comments={'T'},usecols=(0,3,7,11))
Fit50fib2 = np.loadtxt('exp/50fib_2.dat',comments={'T'},usecols=(0,2,4,6,8,10,12))
shear50fib2 = np.loadtxt('exp/par50fib_2.dat',comments={'t'},usecols=(0,1,3,5))
mass50fib2 = np.loadtxt('exp/par50fib_2.dat',comments={'t'},usecols=(0,2,4,6))
F50mon = np.loadtxt('exp/50mon.dat',comments={'T'},usecols=(0,1,3,5))
D50mon = np.loadtxt('exp/50mon.dat',comments={'T'},usecols=(0,2,4,6))
F200fib = np.loadtxt('exp/200fib.dat',comments={'T'},usecols=(0,3,5,7))
D200fib = np.loadtxt('exp/200fib.dat',comments={'T'},usecols=(0,4,6,8))
F200mon = np.loadtxt('exp/200mon.dat',comments={'T'},usecols=(0,1,3,5))
D200mon = np.loadtxt('exp/200mon.dat',comments={'T'},usecols=(0,2,4,6))
In [17]:
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
f = plt.figure()
xticks = np.arange(0,19,3)
# 50 mM fibrils
ax = f.add_subplot(221)
#ax.patch.set_alpha(0.05)
#ax.patch.set_color('b')
plt.ylim(-6.4,5)
ax.text(-28,0.5,'Fibrils',rotation=90,fontsize=20, style='italic')
ax.text(23,6,'$c_s = 0$ M',fontsize=20)
plt.xlim(0,71)
plt.ylabel("$\Delta$F/n / Hz",color='blue')
plt.plot(F50fib[1:,0]/60,F50fib[1:,:],color='blue', lw=2)
plt.plot(D50fib[1:,0]/60,D50fib[1:,:],color='red', lw=2)
plt.plot(Fit50fib[1:,0]/60,Fit50fib[1:,:],color='black', lw=1)
# 200 mM fibrils
ax = f.add_subplot(222)
plt.xticks(xticks)
#ax.patch.set_alpha(0.05)
#ax.patch.set_color('r')
ax.text(4.5,6,'$c_s = 0.15$ M',fontsize=20)
ax.yaxis.tick_right()
plt.ylim(-6.4,5)
plt.xlim(0,18)
plt.ylabel("$\Delta$D / 10$^{-6}$",color='red')
ax.yaxis.set_label_position("right")
plt.plot(F200fib[1:,0]/60,F200fib[1:,:],color='blue', lw=2)
plt.plot(D200fib[1:,0]/60,D200fib[1:,:],color='red', lw=2)
# 50 mM monomers
ax = f.add_subplot(223)
plt.xticks(xticks)
#ax.patch.set_alpha(0.05)
#ax.patch.set_color('r')
ax.text(-7,2,'Monomers',rotation=90,fontsize=20, style='italic')
plt.ylim(-6,5)
plt.xlim(0,18)
plt.ylabel("$\Delta$F/n / Hz",color='blue')
plt.xlabel("$t$/min")
plt.plot(F50mon[1:,0]/60,F50mon[1:,:],color='blue', lw=2)
plt.plot(D50mon[1:,0]/60,D50mon[1:,:],color='red', lw=2)
# 200 mM monomers
ax = f.add_subplot(224)
plt.xticks(xticks)
#ax.patch.set_alpha(0.05)
#ax.patch.set_color('r')
ax.yaxis.tick_right()
plt.ylim(-6,5)
plt.xlim(0,18)
plt.ylabel("$\Delta$D / 10$^{-6}$",color='red')
ax.yaxis.set_label_position("right")
plt.xlabel("$t$/min")
plt.plot(F200mon[2:,0]/60,F200mon[2:,:],color='blue', lw=2)
plt.plot(D200mon[2:,0]/60,D200mon[2:,:],color='red', lw=2)
f.tight_layout()
f.savefig('qcm-d.pdf', bbox_inches='tight', pad_inches=0.3)
plt.show()
In [16]:
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
f = plt.figure()
xticks = np.arange(0,19,3)
# 50 mM fibrils, batch 1
ax = f.add_subplot(221)
plt.ylim(-11,5)
ax.text(25,6,'batch 1',fontsize=16)
plt.xlim(0,72)
plt.ylabel("$\Delta$F/n / Hz",color='blue')
plt.plot(F50fib[1:,0]/60,F50fib[1:,:],color='blue', lw=2)
plt.plot(D50fib[1:,0]/60,D50fib[1:,:],color='red', lw=2)
plt.plot(Fit50fib[1:,0]/60,Fit50fib[1:,:],color='black', lw=1)
# 50 mM fibrils, batch 2
ax = f.add_subplot(222)
ax.text(25,6,'batch 2',fontsize=16)
ax.yaxis.tick_right()
plt.ylim(-11,5)
plt.xlim(0,72)
plt.ylabel("$\Delta$D / 10$^{-6}$",color='red')
ax.yaxis.set_label_position("right")
plt.plot(F50fib2[1:,0]/60,F50fib2[1:,:],color='blue', lw=2)
plt.plot(D50fib2[1:,0]/60,D50fib2[1:,:],color='red', lw=2)
plt.plot(Fit50fib2[1:,0]/60,Fit50fib2[1:,:],color='black', lw=1)
# mass per area and shear modulus, batch 1
ax = f.add_subplot(223)
plt.yticks(np.arange(0,1.8,0.4))
plt.ylim(0,1.7)
plt.xlim(0,72)
plt.ylabel("mass per\n area / ng mm$^{-2}$",color='blue')
plt.xlabel("$t$/min")
plt.plot(mass50fib[20:,0]/60-mass50fib[20,0]/60,
(mass50fib[20:,1:]-mass50fib[20,1]).mean(axis=1),color='blue',lw=2)
plt.fill_between(mass50fib[20:,0]/60-mass50fib[20,0]/60,
(mass50fib[20:,1:]-mass50fib[20,1]).mean(axis=1)
+(mass50fib[20:,1:]-mass50fib[20,1]).std(axis=1),
(mass50fib[20:,1:]-mass50fib[20,1]).mean(axis=1)
-(mass50fib[20:,1:]-mass50fib[20,1]).std(axis=1),color='blue',lw=2,alpha=.3)
plt.plot(shear50fib[20:,0]/60-shear50fib[20,0]/60,
(shear50fib[20:,1:]-shear50fib[20,1]).mean(axis=1),color='red', lw=2)
plt.fill_between(shear50fib[20:,0]/60-shear50fib[20,0]/60,
(shear50fib[20:,1:]-shear50fib[20,1]).mean(axis=1)
+(shear50fib[20:,1:]-shear50fib[20,1]).std(axis=1),
(shear50fib[20:,1:]-shear50fib[20,1]).mean(axis=1)
-(shear50fib[20:,1:]-shear50fib[20,1]).std(axis=1),color='red',lw=2,alpha=.3)
# mass per area and shear modulus, batch 2
ax = f.add_subplot(224)
plt.yticks(np.arange(0,1.8,0.4))
ax.yaxis.tick_right()
plt.ylim(0,1.7)
plt.xlim(0,72)
plt.ylabel("bilayer shear\n modulus / 10$^{5}$ N m$^{-2}$",color='red',size=16)
ax.yaxis.set_label_position("right")
plt.xlabel("$t$/min")
plt.plot(mass50fib2[40:,0]/60-mass50fib2[40,0]/60,
(mass50fib2[40:,1:]-mass50fib2[40,1]).mean(axis=1),color='blue', lw=2)
plt.fill_between(mass50fib2[40:,0]/60-mass50fib2[40,0]/60,
(mass50fib2[40:,1:]-mass50fib2[40,1]).mean(axis=1)
+(mass50fib2[40:,1:]-mass50fib2[40,1]).std(axis=1),
(mass50fib2[40:,1:]-mass50fib2[40,1]).mean(axis=1)
-(mass50fib2[40:,1:]-mass50fib2[40,1]).std(axis=1),color='blue',lw=2,alpha=.3)
plt.plot(shear50fib2[40:,0]/60-shear50fib2[40,0]/60,
(shear50fib2[40:,1:]-shear50fib2[40,1]).mean(axis=1),color='red', lw=2)
plt.fill_between(shear50fib2[40:,0]/60-shear50fib2[40,0]/60,
(shear50fib2[40:,1:]-shear50fib2[40,1]).mean(axis=1)
+(shear50fib2[40:,1:]-shear50fib2[40,1]).std(axis=1),
(shear50fib2[40:,1:]-shear50fib2[40,1]).mean(axis=1)
-(shear50fib2[40:,1:]-shear50fib2[40,1]).std(axis=1),color='red',lw=2,alpha=.3)
f.tight_layout()
f.savefig('qcm-d_SI.pdf', bbox_inches='tight', pad_inches=0.3)
plt.show()
In [ ]: