Notebook

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from scipy.integrate import quad
from scipy.special import eval_laguerre
import itertools
import copy
from image_matrix_helper import compute_master_list, imshow_list, rgb_map, color_to_rgb, list_to_matrix
import random
import time

nb_start = time.time()

Simulation for General Case of Ligand-Receptor Binding¶

In this notebook, we simulate a thermal system of particles of various colors binding onto a grid. We have $R$ different types of particles and particle of type $j$ has $n_j$ copies in the system. Particles can exist both on and off the grid and particles of type $j$ have a binding affinity of $\gamma_j$ to the grid. Each particle type also has a collection of "correct" locations on the grid. A particle of type $j$ binds to this correct location with an additional optimal binding affinity factor of $\delta_j$ (i.e., its net affinity to such a site is $\gamma_j \delta_j$ ). Here we want to use simulations of this system to affirm analytical calculations of the average number of bound particles and the average number of correctly bound particles as functions of temperature.

Numerical representations of analytical work¶

Equations of Large $N$ approximation¶

In the large $N$ limit, the order parameters for the system can be approximated as

$\begin{equation} \langle k \rangle = \sum_{j=1}^r \frac{n_j }{\bar{z} \gamma_j +1} \left(\bar{z} \gamma_j + \frac{\displaystyle L_{n_j-1} \left( \bar{\phi}_{j}\right)}{\displaystyle L_{n_j} \left( \bar{\phi}_{j} \right)}\right) \qquad \langle m \rangle = \sum_{j=1}^r \frac{n_j \delta_j}{\delta_j-1} \frac{\displaystyle L_{n_j-1} \left( \bar{\phi}_{j}\right)}{\displaystyle L_{n_j} \left( \bar{\phi}_{j} \right)}, \end{equation}$

where $\bar{z}$ and $\bar{x}$ are defined as

$\begin{equation} \bar{z} = \sum_{j=1}^{R} \frac{n_j}{\bar{z} \gamma_j +1} \left(1- \frac{\displaystyle L_{n_j-1} \left( \bar{\phi}_{j}\right)}{\displaystyle L_{n_j} \left( \bar{\phi}_{j} \right)}\right), \qquad \bar{x} = \sum_{j=1}^{R} n_j\left(1- \frac{\displaystyle L_{n_j-1} \left( \bar{\phi}_{j}\right)}{\displaystyle L_{n_j} \left( \bar{\phi}_{j} \right)}\right) \end{equation}$ with

$\begin{equation} \bar{\phi}_{j}\equiv \frac{\bar{x}}{1-\delta_j}\left(1+ \frac{1}{\bar{z} \gamma_j}\right). \label{eq:phi_def} \end{equation}$ and

$L_n(x)$ the

$n$ th Laguerre polynomial.

For these simulations we will take

$\gamma_j = (\beta E_V)^{3/2} e^{-\beta E_j}, \qquad \delta_j = e^{\beta \Delta_j}$ where

$E_V$ is a volumetric Boltzmann factor associated with free particles (e.g., for a point-particle

$E_V \equiv h^2/2\pi mV^{2/3}$ ), and

$E_j$ is the binding energy for particles to the grid. We also take where

$\Delta_j$ is the binding energy advantage for particles binding to their correct locations in the grid.

Parameter function definitions¶

In [2]:

# helper function definitions
gamma_func = lambda E0, Ev, T: 4*np.sqrt(2)*np.exp(E0/T)*(Ev/T)**(3/2)
delta_func = lambda Del, T: np.exp(Del/T)
phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta) 

Equilibrium equations¶

In [3]:

def constr_func(X, T, E0s, Dels, Evs, Ns):
    
    """
    Equations constraint equations that determine zbar and xbar
    """
    
    x = X[0]
    z = X[1]
    
    F = np.ones(2)
    
    R = len(Ns)
    gammas_ = gamma_func(E0s,Evs, T)
    deltas_ = delta_func(Dels, T)
    phis_ = phi_func(x, z, gammas_, deltas_)
    
    F[0] = z- np.sum([Ns[j]/(z*gammas_[j]+1)*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])
    F[1] = x- np.sum([Ns[j]*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) ) for j in range(R)])
                     
    return F
    
def m_avg(T, E0s, Dels, Evs, Ns):
    
    """
    Function that computes m_avg
    """    
    
    x, z =  fsolve(constr_func, x0 = (50,500), args = (T, E0s, Dels, Evs, Ns))
    
    R = len(Ns)
    gammas_ = gamma_func(E0s,Evs, T)
    deltas_ = delta_func(Dels, T)
    phis_ = phi_func(x, z, gammas_, deltas_)
    
    return np.sum([Ns[j]*deltas_[j]/(deltas_[j]-1)*eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) for j in range(R)] )
    
def k_avg( T, E0s, Dels, Evs, Ns):
    
    """
    Function that computes k_avg
    """       
    
    x, z = fsolve(constr_func, x0 = (50,500), args = (T, E0s, Dels, Evs, Ns))

    R = len(Ns)
    gammas_ = gamma_func(E0s,Evs, T)
    deltas_ = delta_func(Dels, T)
    phis_ = phi_func(x, z, gammas_, deltas_)    
    
    return np.sum([Ns[j]/(z*gammas_[j]+1)*(z*gammas_[j] + eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])

Important temperatures¶

The temperature $k_BT_{\text{crit}}$ at which the system settles into its completely correct configuration is defined by

$\begin{equation} 1 = \sum_{j=1}^{R}n_j e^{-\beta_{\text{crit}} \Delta_j}\left( 1+ (\beta_{\text{crit}} E_V)^{-3/2} e^{-\beta_{\text{crit}} E_j}\right) + O(e^{-2\beta_{\text{crit}} \Delta_j}). \label{eq:master_therm_subs} \end{equation}$

In [4]:

# general thermal constraint
def master_them_constr(T, E0s, Dels, Evs, Ns):
    
    F =  1-np.sum(Ns*delta_func(Dels, T)**(-1)*(1+gamma_func(E0s, Evs, T)**(-1)))
    
    return F

# critical temperature
kBTcrit_master = lambda E0s, Dels, Evs, Ns: fsolve(master_them_constr, x0 = 100.5, args = (E0s, Dels, Evs, Ns))[0]

Example analytical plot¶

In [5]:

# temperature vals
Tvals = np.linspace(0.1, 3.0, 50)

# parameters for the integral
np.random.seed(42)

R = 50
E0_bar, sigma_E = 10.0, 2.0
Del_bar, sigma_D = 3.0, 1.0
E0s = np.random.randn(R)*sigma_E+E0_bar 
Dels = np.random.randn(R)*sigma_D+Del_bar
Nelems = np.random.randint(1,10,R)
Evs = np.ones(R)*0.001

# computing analytical values of k and m
avg_k_approx_vals = [k_avg(T, E0s, Dels, Evs, Nelems) for T in Tvals]
avg_m_approx_vals = [m_avg(T, E0s, Dels, Evs, Nelems) for T in Tvals]

<ipython-input-2-a582446d2d66>:4: RuntimeWarning: divide by zero encountered in true_divide
  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)
<ipython-input-3-08aba1dbaf99>:17: RuntimeWarning: invalid value encountered in double_scalars
  F[0] = z- np.sum([Ns[j]/(z*gammas_[j]+1)*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j])) for j in range(R)])
<ipython-input-3-08aba1dbaf99>:18: RuntimeWarning: invalid value encountered in double_scalars
  F[1] = x- np.sum([Ns[j]*(1-eval_laguerre(Ns[j]-1, phis_[j])/eval_laguerre(Ns[j], phis_[j]) ) for j in range(R)])
<ipython-input-2-a582446d2d66>:4: RuntimeWarning: invalid value encountered in multiply
  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)
/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)

In [6]:

## plotting order parameters
ax = plt.subplot(111)
ax.plot(Tvals, avg_k_approx_vals/np.sum(Nelems), label = r'$\langle k \rangle$ (large $N$)')
ax.plot(Tvals, avg_m_approx_vals/np.sum(Nelems), label = r'$\langle m \rangle$ (large $N$)')
ax.set_xlabel(r'$k_BT$', fontsize = 18, labelpad = 10.5)
ax.grid(alpha = 0.5)
ax.axvline(x = kBTcrit_master(E0s, Dels, Evs, Nelems), color = 'g', linestyle = '--')

# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

# Only show ticks on the left and bottom spines
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

Metropolis Hastings simulation code¶

Microstate transitions¶

In [7]:

## dissociation operator
def trans_dissoc(free_objs, bound_objs):
        
        # indices of non-empty
        indxs = [i for i, x in enumerate(bound_objs) if x != "-"]
        
        # random choice for bound object
        random_indx = random.choice(indxs)

        ## new state vector 
        free_objs_new = copy.deepcopy(free_objs)
        bound_objs_new = copy.deepcopy(bound_objs)

        # putting empty slot
        bound_objs_new[random_indx] = '-'
        # appending previously bound object to free objects
        free_objs_new.append(bound_objs[random_indx])

        return free_objs_new, bound_objs_new


## association operator
def trans_assoc(free_objs, bound_objs):
    
        # random element to associate
        elem = random.choice(free_objs)

        # indices of empty spaces
        indxs = [i for i, x in enumerate(bound_objs) if x == "-"]  

        # random choice for empty space
        random_indx = random.choice(indxs)

        ## new state vector 
        free_objs_new = copy.deepcopy(free_objs)
        bound_objs_new = copy.deepcopy(bound_objs)    

        ## state
        free_objs_new.remove(elem)
        bound_objs_new[random_indx] = elem      

        return free_objs_new, bound_objs_new

## permutation operator
def trans_perm(free_objs, bound_objs):
    
    Ncomp = len(bound_objs)
    i1 = int(random.choice(range(Ncomp)))
    i2 = int(random.choice(range(Ncomp)))

    ## new omega vector 
    bound_objs_new = copy.deepcopy(bound_objs)
    bound_objs_new[i2] = bound_objs[i1]
    bound_objs_new[i1] = bound_objs[i2]
    
    return free_objs, bound_objs_new

Logarithm of Botlzmann factor¶

The logarithm of the Botlzmann factor for a microstate (i.e., the temperature normalized negative energy of the microstate) is defined as

$\begin{equation} \beta E(\boldsymbol{k}, \boldsymbol{m}) = \sum_{i=1}^R(m_i \ln \delta_i + k_i \ln \gamma_i). \label{eq:sim_en} \end{equation}$

In [8]:

def log_boltz(free_objs, bound_objs, mstr_vec, deltas, gammas, name_key):
    
    elem_set = list(set(mstr_vec))
    count_dict = dict()
    for elem in elem_set:
        count_dict[elem] = bound_objs.count(elem)

    bind_log_factor = 0
    for elem in elem_set:
        key = name_key[elem]
        bind_log_factor += count_dict[elem]*np.log(gammas[key])
        
        
    corr_log_factor = 0
    for j in range(len(bound_objs)):
        if bound_objs[j] == mstr_vec[j]:
            elem = bound_objs[j]
            key = name_key[elem]
            corr_log_factor+=np.log(deltas[key])        
    
    return bind_log_factor+corr_log_factor    

Function to count the number of correctly bound particles¶

In [9]:

def m_calc(bound_objs, mstr_vec):
    
    num = 0
    for k in range(len(mstr_vec)):
        if mstr_vec[k] == bound_objs[k]:
            num += 1
    return num

Checking logarithm of Boltzmann factor definition¶

In [10]:

# defining name key
name_key0 = dict()
key_list = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', ]
for j in range(len(key_list)):
    name_key0[key_list[j]] = j
    
# random energies    
np.random.seed(2)
q1 = np.random.rand(10)
q2 = np.random.rand(10)

# sample master list
sample_master = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', ]
bound_init_0 = ['-', '-', '-', '-', '-', '-', 'G', 'H', 'I', 'J', ]
free_init_0 = ['A', 'B', 'C', 'D', 'E', 'F' ]

print('Energy for original ordering:', log_boltz( free_init_0,bound_init_0, sample_master, deltas = np.exp(-q1), gammas = np.exp(-q2), name_key =  name_key0))
e1 = -np.sum([q1[k] for k in range(len(sample_master)) if sample_master[k] == bound_init_0[k]])
e2 = -np.sum([q2[k] for k in range(len(sample_master)) if sample_master[k] in bound_init_0])
print('Checking  energy value:', e1+e2)
print('Number of correctly placed elements:', m_calc(bound_init_0, sample_master))
print('Number of bound elements:', np.sum([1 for elem in bound_init_0 if elem!='-']))
print('-----')
random.seed(1)
free_init_0, perm_bound = trans_perm(free_init_0, bound_init_0)
free_init_new, bound_init_new = trans_assoc(free_init_0, perm_bound)
print('Energy after permutation and associaation values:', log_boltz(free_init_new, bound_init_new, sample_master, deltas = np.exp(-q1), gammas = np.exp(-q2), name_key =  name_key0))
e1 = -np.sum([q1[k] for k in range(len(sample_master)) if sample_master[k] == bound_init_new[k]])
e2 = -np.sum([q2[k] for k in range(len(sample_master)) if sample_master[k] in bound_init_new])
print('Checking  energy value:', e1+e2)
print('Number of correctly placed elements:', m_calc(bound_init_new, sample_master))
print('Number of bound elements:', np.sum([1 for elem in bound_init_new if elem!='-']))

Energy for original ordering: -3.664820641553842
Checking  energy value: -3.664820641553841
Number of correctly placed elements: 4
Number of bound elements: 4
-----
Energy after permutation and associaation values: -4.01912719922027
Checking  energy value: -4.01912719922027
Number of correctly placed elements: 3
Number of bound elements: 5

Metropolis Hastings algorithm¶

In [11]:

### Metropolis Monte Carlo Algorithm 

## loads uniform random sampling 
runif = np.random.rand

def met_assembly_grid(Niter, free_objs, bound_objs, mstr_vec, deltas, gammas, name_key):
    '''
    #################################################################
    # function to sample using Metropolis 
    #  
    # n_iter:  number of iterations
    # initial_state: initial state for the start position for our chain
    # gamma: energy cost for incorrect component
    # temp: temperature 
    ##################################################################
    '''
    
    # Initialize state values
    free_objs_vals = [0]*(Niter+1)
    bound_objs_vals = [0]*(Niter+1)
    
    # Set initial values
    free_objs_vals[0] = free_objs[:]
    bound_objs_vals[0] = bound_objs[:]
        
    # Initialize acceptance counts
    # We can use this to tune our number of steps
    accepted = 0
    
    # debugging code
    debug_assoc, debug_dissoc, debug_perm = 0, 0, 0
    
    for i in range(Niter):
    
        # get current monomer and dimer states
        current_free_objs = copy.deepcopy(free_objs_vals[i])
        current_bound_objs = copy.deepcopy(bound_objs_vals[i])
        
        N_free = len(current_free_objs)
        N_bound = len(current_bound_objs)-len(current_free_objs)
        
        u_trans = runif()
        
        if u_trans < 1/3: #first type of transition; monomer association 
            
            if N_free < 1:
                
                log_alpha = np.log(1e-15) 
            
            else:
                # proposed new monomer and dimer states
                new_free_objs, new_bound_objs = trans_assoc(current_free_objs, current_bound_objs)

                # transition elements
                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)
                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)

                # weight
                num = N_free*N_free
                den = N_bound+1

                # Log-acceptance rate
                log_alpha = log_final-log_init+np.log(num/den) 
            
        elif 1/3 <= u_trans < 2/3: #second type of transition; bound monomer dissociation
            
            if N_bound <1:
                
                log_alpha = np.log(1e-15) 
                
            else: 
                
                # proposed new monomer and dimer states
                new_free_objs, new_bound_objs = trans_dissoc(current_free_objs, current_bound_objs)

                # transition elements
                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)
                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)

                # weight
                num = N_bound 
                den = (N_free+1)*(N_free+1)   

                # Log-acceptance rate
                log_alpha = log_final-log_init+np.log(num/den) 
            
        elif 2/3 <= u_trans: #third type of transition; switching bounded elements
            
            if N_bound <2:
                
                log_alpha = np.log(1e-15) 
            
            else:
                
                # proposed new monomer and dimer states
                new_free_objs, new_bound_objs = trans_perm(current_free_objs, current_bound_objs)

                # transition elements
                log_init = log_boltz(current_free_objs, current_bound_objs, mstr_vec, deltas, gammas, name_key)
                log_final = log_boltz(new_free_objs, new_bound_objs, mstr_vec, deltas, gammas, name_key)

                # Log-acceptance rate
                log_alpha = log_final-log_init
            
        # Sample a uniform random variate
        u = runif()

        # Test proposed value
        if np.log(u) < log_alpha:
            # Accept
            free_objs_vals[i+1] = new_free_objs
            bound_objs_vals[i+1] = new_bound_objs
            
            #log_current_prob = log_proposed_prob
            accepted += 1
        else:
            # Stay put
            free_objs_vals[i+1] = free_objs_vals[i]
            bound_objs_vals[i+1] = bound_objs_vals[i]

    # return our samples and the number of accepted steps
    return free_objs_vals, bound_objs_vals, accepted

Computing microstate averages from simiulations¶

In [12]:

def avg_k(bound_objs_vals, Nmc):
    
    """
    Microstate average of number of bound objects
    We only consider microstates near the end of theh chain to ensure
    that the system has equilibrated
    """
    
    length = int(Nmc/50)
    
    ls = [0]*length
    ls = np.array(ls)
    for k in range(length):
        ls[k] = len(bound_objs_vals[Nmc-length+k]) - bound_objs_vals[Nmc-length+k].count('-')
    
    return(np.mean(ls))

# average number of correctly bound objects
def avg_m(bound_objs_vals, mstr_vec, Nmc):
    
    """
    Microstate average of number of correctly bound objects
    We only consider microstates near the end of theh chain to ensure
    that the system has equilibrated
    """ 
    
    length = int(Nmc/50)
    
    ls = [0]*length
    ls = np.array(ls)
    for k in range(length):
        ls[k] =  np.sum([1 for j in range(len(mstr_vec)) if bound_objs_vals[Nmc-length+k][j]==mstr_vec[j]])
    
    return(np.mean(ls))

Image grid for completely correct configuration¶

In [13]:

# defining master_list
master_list =compute_master_list()
# testing plot
imshow_list(master_list, title = 'Completely Correct Configuration');
# defining Nelems
Nelems = np.zeros(8)
key_list = list(rgb_map.keys())[:-1]
name_key_ = dict()
for j in range(len(key_list)):
    name_key_[key_list[j]] = j
    Nelems[j] = master_list.count(key_list[j])

In [14]:

# displaying copy-number counts of the various elements
Nelems

Out[14]:

array([ 9.,  9., 10.,  5.,  7.,  6.,  3., 51.])

Simulating system¶

In [15]:

## Generate lf for each temperature from .03 to 2.0 in npoints steps
t0 = time.time()

# number of steps for MC algortihm
Nmc = 30000
    
# parameter definitions    
R = 8
np.random.seed(24)
Del_bar, sigma_D = 3.0, 1.0
Dels = np.random.randn(R)*sigma_D+Del_bar
E0_bar, sigma_E = 14.0, 2.0
E0s = np.random.randn(R)*sigma_E+E0_bar
Evs = np.ones(R)*0.001    
    
# initial monomer and dimer states; 
# system in microstate of all correct dimers
random.seed(0)
free_objs_0 = []
bound_objs_0 = random.sample(master_list, len(master_list))
mstr_vec = copy.deepcopy(master_list)

# temperature limits
Tmin = .05
Tmax = 3.0

npoints = 15 #number of temperature values
navg = 5 # number of times we run simulation at each temperature; 50 in paper
temp_vals = np.linspace(Tmin, Tmax, npoints).tolist()

# list of dimer values 
sim_k_vals = [0]*npoints
# list of correct dimer values
sim_m_vals = [0]*npoints
# accepted list 
accepted_list = [0]*npoints
# saved list for plotting
saved_list = dict()

for k in range(npoints):
    
    fin_k_vals = [0]*navg
    fin_m_vals = [0]*navg
    fin_accepted = [0]*navg
    
    for j in range(navg): 
        
        # make copy of initial monomer and dimer states 
        free_objs_copy = copy.deepcopy(free_objs_0)
        bound_objs_copy = copy.deepcopy(bound_objs_0)
        
        # defining helper functions
        gammas_ = gamma_func(E0s, Evs, temp_vals[k])
        deltas_ = delta_func(Dels, temp_vals[k])  

        # metroplois generator
        free_list, bound_list, accepted = met_assembly_grid(Nmc, 
                                                            free_objs_copy,
                                                            bound_objs_copy,
                                                            mstr_vec,
                                                            deltas_,
                                                            gammas_,
                                                            name_key_) 
        
        

        
        # averaging final states to compute observables
        fin_k_vals[j] = avg_k(bound_list, Nmc)
        fin_m_vals[j] = avg_m(bound_list, mstr_vec, Nmc)
        fin_accepted[j] = accepted
    
    # saving every 5 temperatures
    if (k+1)%5 ==0 or k ==0:
        saved_list[k] = ['white' if x=='-' else x for x in bound_list[-1]]   
    
    # averaging over computed equilibrium averages
    sim_k_vals[k] = np.mean(np.array(fin_k_vals))
    sim_m_vals[k] = np.mean(np.array(fin_m_vals))
    accepted_list[k] = np.mean(np.array(fin_accepted))
    
    t_prelim = time.time()
    print("Temperature Run:",str(k+1),"; Current Time:", round(t_prelim-t0,2),"secs")
    

t1 = time.time()
print("Total Simulation Run Time:",t1-t0,"secs")

Temperature Run: 1 ; Current Time: 62.85 secs
Temperature Run: 2 ; Current Time: 125.33 secs
Temperature Run: 3 ; Current Time: 184.85 secs
Temperature Run: 4 ; Current Time: 236.61 secs
Temperature Run: 5 ; Current Time: 283.27 secs
Temperature Run: 6 ; Current Time: 341.44 secs
Temperature Run: 7 ; Current Time: 406.55 secs
Temperature Run: 8 ; Current Time: 471.81 secs
Temperature Run: 9 ; Current Time: 529.68 secs
Temperature Run: 10 ; Current Time: 583.34 secs
Temperature Run: 11 ; Current Time: 633.34 secs
Temperature Run: 12 ; Current Time: 682.38 secs
Temperature Run: 13 ; Current Time: 730.92 secs
Temperature Run: 14 ; Current Time: 780.55 secs
Temperature Run: 15 ; Current Time: 831.61 secs
Total Simulation Run Time: 831.6078369617462 secs

Simulated image grid at various temperatures¶

In [16]:

# figure parameters
rows, cols, idx = 2, 2, 0
fig, axes = plt.subplots(nrows=rows, ncols=cols, figsize=(9,9))

# list of keys for saved snapshots of image
img_key_list = list(saved_list.keys())

for i in range(rows):
    for j in range(cols):
        if idx < 4:
            axes[i, j].imshow(color_to_rgb(list_to_matrix(saved_list[img_key_list[idx]])))
            ax = plt.gca()
            
            # making ticks invisible
            axes[i, j].set_xticks([])
            axes[i, j].set_yticks([])

            # Minor ticks
            axes[i, j].set_xticks(np.arange(-0.5, 11, 1), minor=True)
            axes[i, j].set_yticks(np.arange(-0.5, 10, 1), minor=True)
            axes[i, j].tick_params(axis='y', colors='red')
            
            # labeling images
            itimes = 'i'*(1+idx) if idx<3 else 'iv'
            
            # Gridlines based on minor ticks
            axes[i, j].grid(which='minor', color='w', linestyle='-', linewidth=3)
            axes[i, j].set_title(fr'({itimes})  $k_BT = {round(temp_vals[img_key_list[idx]],2)}$', fontsize = 18, y = -.2)

            # making spines invisible
            axes[i, j].spines['right'].set_visible(False)
            axes[i, j].spines['top'].set_visible(False)
            axes[i, j].spines['left'].set_visible(False)
            axes[i, j].spines['bottom'].set_visible(False)  
            
            idx +=1
    
# plt.savefig('general_grid_assembly_grid_plots.png', bbox_inches='tight', format = 'png')      
plt.show()

Comparing analytical and simulation results¶

In [17]:

plt.figure(figsize = (7,5))
ax = plt.subplot(111)

# simulation results
plt.plot(temp_vals,np.array(sim_k_vals)/np.sum(Nelems), 
        label = r'Sim. $\langle k \rangle$/N',
        markersize = 7.5,
       marker = 'D',
       linestyle = '')
plt.plot(temp_vals,np.array(sim_m_vals)/np.sum(Nelems), 
        label = r'Sim. $\langle m \rangle$/N',
        markersize = 7.5,
       marker = 's',
       linestyle = '')

# large N analytical results
k_avg_approx_vals = [k_avg(T, E0s, Dels, Evs, Nelems)/np.sum(Nelems) for T in Tvals]
m_avg_approx_vals = [m_avg(T, E0s, Dels, Evs, Nelems)/np.sum(Nelems) for T in Tvals]
plt.plot(Tvals, k_avg_approx_vals, label = r'Large $N$ $\langle k \rangle$/N', linestyle= '--', linewidth = 3.0)
plt.plot(Tvals, m_avg_approx_vals, label = r'Large $N$ $\langle m \rangle$/N', linewidth = 2.0 )
ax.axvline(x = kBTcrit_master(E0s, Dels, Evs, Nelems), color = 'k', linestyle = '-.', linewidth = 2)


plt.legend(loc = 'best', fontsize = 12)
# plot formatting
ax.set_xlabel(r'$k_B T$', fontsize = 18)
plt.xlim([-0.01,3.2])
plt.ylim([0,1.1])
plt.grid(alpha = 0.45)


# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

# Only show ticks on the left and bottom spines
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

# increase label size
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=12)
ax.text(kBTcrit_master(E0s, Dels, Evs, Nelems)-.2, 0.25, r'$k_BT_{crit}$', color='black', fontsize = 14.5,
        bbox=dict(facecolor='white', edgecolor='none', pad=5.0))


for i in range(4):
    ax.text(temp_vals[img_key_list[i]], sim_k_vals[img_key_list[i]]/np.sum(Nelems)+.05,'('+'i'*(1+i)+')' if i<3 else '(iv)', fontsize = 14 ) 
    
# plt.savefig(f'general_grid_assembly.png', bbox_inches='tight', format = 'png')         

<ipython-input-2-a582446d2d66>:4: RuntimeWarning: divide by zero encountered in true_divide
  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)
<ipython-input-2-a582446d2d66>:4: RuntimeWarning: invalid value encountered in multiply
  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)
/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
<ipython-input-2-a582446d2d66>:4: RuntimeWarning: overflow encountered in true_divide
  phi_func = lambda x, z, gamma, delta: x*(1+ 1/(z*gamma))/(1-delta)
/Users/mobolajiwilliams/opt/miniconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

In [18]:

print('Total Notebook Runtime: %.3f mins' % ((time.time()-nb_start)/60))

Total Notebook Runtime: 13.918 mins