<jupyter_start><jupyter_text>Modules and Constants<jupyter_code># Mathematical operations
import numpy as np
import scipy as sp
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
%matplotlib inline

# My own functions
from diffusion import *
# Natural constants
kB    = 1.3806485*10**-23  # Boltzmann constant [J/K]

# Simulation parameters
T     = 310<jupyter_output><empty_output><jupyter_text>Box properties<jupyter_code>old_box_data = np.loadtxt("simulations-ANT1-constant-density/results-translation/dc_Protein.dat").T
new_box_data = np.loadtxt("simulations-ANT1-dilute/results-translation/dc_Protein.dat").T
print('number of old simulations:', old_box_data.shape[1])
print('number of new simulations:', new_box_data.shape[1])
old_box_edge   = old_box_data[0]
old_box_height = (old_box_data[1]-4.5)/2.0
new_box_edge   = new_box_data[0]
new_box_height = (new_box_data[1]-4.5)/2.0<jupyter_output><empty_output><jupyter_text>MSD and diffusion coefficients<jupyter_code>numblocks = 10
fitstart  = 30 # 50 # 200
fitend    = 60 # 125 # 400
old_msd = []
old_dc  = []
old_err = []

mult   = ['01','01a','01b','02','02a','02b','03','03a','03b','04','05','07','10','15','20','25','30']

# Iterate through the old data
for n_, m_ in enumerate(mult):
    
    data = np.loadtxt('simulations-ANT1-constant-density/results-rotation/Rotation_Protein_PC60PE40-PC50PE40CL10-1OKC-piece12-mult'+m_+'.xvg',comments=['#','@']).T
    
    # The first column is time
    time = data[0]
    
    # overwrite the increments with their cumulative sum ( = total rotation in radians)
    rota = np.cumsum(data[1:],axis=1)
    
    # Calculate MSD and diffusion coefficient with uncertainty
    totdc, diffcoeff, err_diffc, t, msd_i = compute_dc(time,rota,numblocks,fitstart,fitend)
    
    old_msd.append( msd_i     )
    old_dc.append(  totdc     )
    old_err.append( err_diffc )

# Convert to arrays
old_msd = np.array(old_msd)
old_dc  = np.array(old_dc )
old_err = np.array(old_err)
new_msd = []
new_dc  = []
new_err = []

piece  = ['07','08','09','10','12','14','16','18','20','24','28']
letter = ['A','B','C','D','E','F']

# Iterate through the new data
for np_, p_ in enumerate(piece):
    for nl_, l_ in enumerate(letter):
        
        data = np.loadtxt('simulations-ANT1-dilute/results-rotation/Rotation_Protein_PC60PE40-PC50PE40CL10-1OKC-piece'+p_+'-mult01-'+l_+'.xvg',comments=['#','@']).T

        # The first column is time
        time = data[0]

        # overwrite the increments with their cumulative sum ( = total rotation in radians)
        rota = np.cumsum(data[1:],axis=1)

        # Calculate MSD and diffusion coefficient with uncertainty
        totdc, diffcoeff, err_diffc, t, msd_i = compute_dc(time,rota,numblocks,fitstart,fitend)

        new_msd.append( msd_i   )
        new_dc.append(  totdc     )
        new_err.append( err_diffc )

# Convert to arrays
new_msd = np.array(new_msd)
new_dc  = np.array(new_dc )
new_err = np.array(new_err)<jupyter_output><empty_output><jupyter_text>Average the values from the different trajectories at same size<jupyter_code>def old_average_init(arr):
    arr_av = []
    arr_av.append( np.mean(arr[ :3], axis=0 ) )
    arr_av.append( np.mean(arr[3:6], axis=0 ) )
    arr_av.append( np.mean(arr[6:9], axis=0 ) )
    for i in range(9,17):
        arr_av.append( arr[i] )
    return np.array(arr_av)
        
old_edge_av = old_average_init(old_box_edge)
old_msd_av  = old_average_init(old_msd)
old_dc_av   = old_average_init(old_dc)
old_err_av  = old_average_init(old_err)

# Correct the error estimate for sizes with multiple simulations
old_err_av[:3] = old_err_av[:3]/np.sqrt(3)
def new_average_init(arr):
    arr_av = []
    for j in range(0,66,6):
        arr_av.append( np.mean(arr[j:j+6], axis=0) )
    arr_av = np.array(arr_av)
    return np.array(arr_av)

new_edge_av = new_average_init(new_box_edge)
new_msd_av  = new_average_init(new_msd)
new_dc_av   = new_average_init(new_dc)
new_err_av  = new_average_init(new_err)

# Correct the error estimate for sizes with multiple simulations
new_err_av[-11:] = np.array([np.std(new_dc[i:i+6])/np.sqrt(6) for i in range(11)])
print( 'Total simulation time:', time[-1] )
print( 'Total number of frames:', len(time) ) 
ts = time[-1]/(len(time)-1) 
print( 'Time step between two frames:', ts, 'ps' )
blocksize = (len(time)-1)/numblocks
blocktime = blocksize*ts
print( 'Size of one block:', blocksize, 'frames = ', blocktime, 'ps' )<jupyter_output>Total simulation time: 2000000.0
Total number of frames: 20001
Time step between two frames: 100.0 ps
Size of one block: 2000.0 frames =  200000.0 ps<jupyter_text>Save the data<jupyter_code>np.savetxt('rotation_dc_ant1_old.dat',np.array([old_edge_av,old_dc_av,old_err_av]).T,header="# box width [nm], diff. coeff. [rad^2/ps], err(diff. coeff.) [rad^2/ps]" ) 
np.savetxt('rotation_dc_ant1_new.dat',np.array([new_edge_av,new_dc_av,new_err_av]).T,header="# box width [nm], diff. coeff. [rad^2/ps], err(diff. coeff.) [rad^2/ps]" )<jupyter_output><empty_output><jupyter_text>Fit to Theory<jupyter_code>def dc_rot_theory(width,radius,d0):
    return d0*(1 - np.pi*radius**2/width**2)
old_edge_av,old_dc_av,old_err_av = np.loadtxt('rotation_dc_ant1_old.dat',comments='#').T
old_res_opt = sp.optimize.curve_fit(dc_rot_theory,old_edge_av,old_dc_av,sigma=old_err_av,absolute_sigma=True)
old_rad_opt = old_res_opt[0][0]
old_dc0_opt = old_res_opt[0][1]
old_rad_err = old_res_opt[1][0,0]
old_dc0_err = old_res_opt[1][1,1]
print(r"R_H: %2.2f +/- %2.2f nm"%(old_rad_opt,old_rad_err))
print(r"D_0: %2.11f +/- %2.11f rad^2/us"%(old_dc0_opt*1e6,old_dc0_err*1e6))
old_etah = 1e11*(kB*T)/(4*np.pi*old_dc0_opt*1e-6*old_rad_opt**2)
print(r"eta*h: %2.11f * 10^-11 Pa s m" %old_etah )
new_edge_av,new_dc_av,new_err_av = np.loadtxt('rotation_dc_ant1_new.dat',comments='#').T
new_res_opt = sp.optimize.curve_fit(dc_rot_theory,new_edge_av,new_dc_av,sigma=new_err_av,absolute_sigma=True)
#res_opt = sp.optimize.curve_fit(dc_rot_theory,edge_av[11:],dc_av[11:],sigma=err_av[11:],absolute_sigma=True)
new_rad_opt = new_res_opt[0][0]
new_dc0_opt = new_res_opt[0][1]
new_rad_err = new_res_opt[1][0,0]
new_dc0_err = new_res_opt[1][1,1]
print(r"R_H: %2.2f +/- %2.2f nm"%(new_rad_opt,new_rad_err))
print(r"D_0: %2.11f +/- %2.11f rad^2/us"%(new_dc0_opt*1e6,new_dc0_err*1e6))
new_etah = 1e11*(kB*T)/(4*np.pi*new_dc0_opt*1e-6*new_rad_opt**2)
print(r"eta*h: %2.11f * 10^-11 Pa s m" %new_etah )<jupyter_output>R_H: 2.53 +/- 0.00 nm
D_0: 1.62739014216 +/- 0.00000000049 rad^2/us
eta*h: 3.27883025978 * 10^-11 Pa s m<jupyter_text>Figure<jupyter_code>old_cmap = plt.cm.viridis
new_cmap = plt.cm.viridis
old_labels = [ 'L = %3.1f nm'%e for e in old_edge_av ]
new_labels = [ 'L = %3.1f nm'%e for e in new_edge_av ]
fig, ax = plt.subplots(2,2, figsize=[8,7.8], dpi=300)

# Diffusion coefficients
for i in range(len(old_edge_av)):
    ax[0,1].errorbar(old_edge_av[i],old_dc_av[i]*1e6,yerr=old_err_av[i]*1e6,ls='none',marker='o',capsize=5,c=old_cmap(0.1+0.9*float(i)/len(old_edge_av)))
for j in range(len(new_edge_av)):
    ax[0,0].errorbar(new_edge_av[j],new_dc_av[j]*1e6,yerr=new_err_av[j]*1e6,ls='none',marker='o',capsize=5,c=new_cmap(0.1+0.9*float(j)/len(new_edge_av)))
    
# Fit to theory and limiting value
xax = np.arange(501)
ax[0,1].plot(xax[1:],dc_rot_theory(xax[1:],old_rad_opt,old_dc0_opt)*1e6,c='black',linestyle='--',lw=2,label='fit to theory')
ax[0,0].plot(xax[1:],dc_rot_theory(xax[1:],new_rad_opt,new_dc0_opt)*1e6,c='black',linestyle='--',lw=2,label='fit to theory')
ax[0,1].plot(xax,np.ones(len(xax))*old_dc0_opt*1e6,c='gray',linestyle=':',lw=2,label=r'$D_0$ from fit')
ax[0,0].plot(xax,np.ones(len(xax))*new_dc0_opt*1e6,c='gray',linestyle=':',lw=2,label=r'$D_0$ from fit')

# Titles
title0 = ax[0,0].text(0.03,1.15,'One ANT1 protein in the box',   transform=ax[0,0].transAxes, fontsize=14, va='top',color='black')
title1 = ax[0,1].text(0.01,1.15,'Constant ANT1 protein density', transform=ax[0,1].transAxes, fontsize=14, va='top',color='black')
ax[0,1].legend(loc='upper left')

# Axis options
for axi in ax[0,:]:
    #axi.set_facecolor('lightgray')
    axi.set_alpha(0.5)
    axi.set_xscale('log')
    axi.set_xlim(5,500)
    #ax[0,0].set_xlim(0,400)
    axi.set_ylim(0.65,1.85)
    axi.set_xlabel(r'box width [nm]',size=12)
    axi.set_ylabel(r'$D$ [rad$^2$/$\mathrm{\mu}$s]',size=12)
    #axi.grid('on')

# MSD curves old
for i,msdi in enumerate( old_msd_av ):
    ax[1,1].plot(time*1e-3,msdi*1e3,c=old_cmap(0.1+0.9*float(i)/len(old_msd_av)),lw=2,label=old_labels[i])
ax[1,1].axvspan(fitstart*ts*1e-3,fitend*ts*1e-3, facecolor='lightgray', alpha=1)
#ax[1,0].set_yticks(np.arange(0,0.0,0.02))
ax[1,1].set_xlabel('lag time [ns]',size=12)
ax[1,1].set_ylabel('MSD [$10^{-3}$ rad$^2$]',size=12)

# MSD curves new
for i,msdi in enumerate( new_msd_av ):
    ax[1,0].plot(time*1e-3,msdi*1e3,c=new_cmap(0.1+0.9*float(i)/len(new_msd_av)),lw=2,label=new_labels[i])
ax[1,0].axvspan(fitstart*ts*1e-3,fitend*ts*1e-3, facecolor='lightgray', alpha=1)
ax[1,0].set_xlabel('lag time [ns]',size=12)
ax[1,0].set_ylabel('MSD [$10^{-3}$ rad$^2$]',size=12)


for axi in ax[1,:]:
    axi.set_xlim(0.0,7.0)
    axi.set_ylim(0.0,32.5)
    axi.legend(loc='upper left',fontsize=8)

# Panel labels
ax[0,0].text(0.9, 0.95, 'A', transform=ax[0,0].transAxes, fontsize=16, fontweight='bold', va='top')
ax[0,1].text(0.9, 0.95, 'B', transform=ax[0,1].transAxes, fontsize=16, fontweight='bold', va='top')
ax[1,0].text(0.9, 0.95, 'C', transform=ax[1,0].transAxes, fontsize=16, fontweight='bold', va='top')
ax[1,1].text(0.9, 0.95, 'D', transform=ax[1,1].transAxes, fontsize=16, fontweight='bold', va='top')

fig.tight_layout()

# TO DO: Include the fit to MSD curves

fig.savefig( 'rotational_diffusion_ANT1.pdf', bbox_inches='tight', bbox_extra_artists=(title1,title0), dpi=300 )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18154003059663437, 'fraction_numerical': 0.06160122386537481, 'mean_word_length': 5.762758620689655, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '388643', 'license': 'MIT', 'n_tokens_mistral': 4599, 'n_tokens_neox': 4310, 'n_words': 641, 'path': 'Rotation_ANT1.ipynb', 'repo_name': 'bio-phys/rotmemdiff', 'size': '9805'}