notebook.community

Edit and run 





In [109]:



    


%matplotlib notebook
import scipy.stats as stats
from scipy.optimize import curve_fit

import numpy as np
import matplotlib.pyplot as plt
import particlesim.api
import particlesim.helpers_for_tests
import particlesim.utils.xyz
import particlesim.utils.config_parser
import particlesim.utils.conversion

from mpl_toolkits.mplot3d import Axes3D



    











In [3]:



    


def plot_nacl(traj,left,right,num_na,traj_sample = -1):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    last_na_pos = traj[traj_sample,:num_na,:]
    last_cl_pos = traj[traj_sample,num_na:,:]
    
    a =(last_na_pos <= right)*(last_na_pos >=left)
    a = a[:,0]*a[:,1]*a[:,2]
    small_box_na = last_na_pos[a]
    b = (last_cl_pos <= right)*(last_cl_pos >=left)
    b = b[:,0]*b[:,1]*b[:,2]
    
    small_box_cl = last_cl_pos[b]
    ax.scatter(small_box_na[:,0],small_box_na[:,1],small_box_na[:,2],c='r')
    ax.scatter(small_box_cl[:,0],small_box_cl[:,1],small_box_cl[:,2],c='b')
    ax.set_xlim([left,right])
    ax.set_ylim([left,right])
    ax.set_zlim([left,right])



    











In [4]:



    


def plot_systemconfig(system_config):

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    num_na = len(system_config.xyz)//2
    ax.scatter(system_config.xyz[:num_na,0],system_config.xyz[:num_na,1],system_config.xyz[:num_na,2],c='r')
    ax.scatter(system_config.xyz[num_na:,0],system_config.xyz[num_na:,1],system_config.xyz[num_na:,2],c='b')
    ax.set_xlim([0,system_config.box_size])
    ax.set_ylim([0,system_config.box_size])
    ax.set_zlim([0,system_config.box_size])



    











In [5]:



    


def beta_fct(beta_start, beta_end, iteration_number):
    x = np.linspace(0,np.pi,iteration_number)
    beta = (beta_end - beta_start)*np.exp(-x)*(0.8+0.2*np.cos(15*x)) + beta_start
    return beta[::-1]


def multiple_beta(beta_list,iteration_list, linspace=False):
    beta=np.asarray([])
    if linspace:
        for i in range(len(beta_list)):
            beta = np.append(beta, 1.0 / np.linspace(1.0 / beta_list[i][1], 1.0 / beta_list[i][0], iteration_list[i])[::-1])
    else:
        for i in range(len(beta_list)):
            beta = np.append(beta,beta_fct(beta_list[i][0], beta_list[i][1], iteration_list[i]))
    return beta

def plot_beta(beta, T=False):
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    if T:
        temp = particlesim.utils.conversion.beta_to_kelvin(beta)
        ax1.plot(temp,linestyle = 'dashed',c='r',label="T")
        ax1.set_ylabel(r"$T$ [K]", fontsize = 20)
        ax1.legend(loc=2)
        ax = ax1.twinx()
        ax.plot(beta,'+',label="beta")
        ax.set_ylabel(r"$\beta$ [mol/kcal]", fontsize = 20)
        ax.set_xlabel(r"Iteration", fontsize = 20)
        ax.legend()
    else:
        ax1.plot(beta,'+',label="beta")
        ax1.set_ylabel(r"$\beta$ [mol/kcal]", fontsize = 20)
        ax1.set_xlabel(r"Iteration", fontsize = 20)
        ax1.legend()

Config File:

[general] box-size = 12

[particle_class_1] label = Natrium type = Na charge = 1 distribution = uniform number = 4

[particle_class_2] label = Chlor type = Cl charge = -1 distribution = uniform number = 4





In [6]:



    


creator = particlesim.utils.config_parser.ProblemCreator("/home/mark/Dokumente/Studium/Master/WS1617/CompSci/compscie-mc/jupyter_notebooks/config/8_particle_nacl_rand.cfg")



    











In [7]:



    


system_config = creator.generate_problem()



    











In [8]:



    


plot_systemconfig(system_config)



    


















    




















    























In [9]:



    


sampler = particlesim.api.Sampler(system_config)

Global parameter for simulated annealing





In [10]:



    


iteration_number = 30000
step = 0.01
beta = [0.1,1.3]

Call metropolis simulated annealing function from sampler





In [11]:



    


%%time
traj,pot = sampler.metropolis_sa(iteration_number=iteration_number,step=0.01,beta=beta)



    


















    






CPU times: user 15.7 s, sys: 0 ns, total: 15.7 s
Wall time: 15.7 s

















In [12]:



    


plot_nacl(traj=traj,left=0,right=12,num_na=4,traj_sample=-1)



    


















    




















    























In [13]:



    


beta_values = 1.0 / np.linspace(1.0 / beta[1], 1.0 / beta[0], iteration_number)[::-1]



    











In [14]:



    


fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(pot)
ax.set_xlabel(r"Iteration", fontsize=15)
ax.set_ylabel(r"Potential [kcal/mol]", fontsize=15)
ax.set_ylim([pot[-1]-200,0])
ax2 = ax.twinx()
ax2.plot(beta_values,c='r')
ax2.set_ylabel(r"$\beta$ [mol/kcal]",fontsize=15)
ax.set_title('Potential curve over iterationsteps')
print(pot[-1])



    


















    




















    
















    






-1419.10580732

















In [15]:



    


plot_beta(beta_values,T=True)

create beta values for double simulated annealing of system configuration





In [16]:



    


beta_list=[[0.1,0.9]]+[[0.3,1.]]
iter_list=[15000,15000]
#beta_list=[beta]
#iter_list= [iteration_number]
beta_many = multiple_beta(beta_list,iter_list, linspace=True)
plot_beta(beta_many,T=True)



    


















    




















    























In [17]:



    


traj_many, pot_many = sampler.metropolis_sa(iteration_number=iteration_number,step=0.01,beta=beta_many)



    











In [18]:



    


fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(pot_many)
ax.set_xlabel(r"Iteration", fontsize=15)
ax.set_ylabel(r"Potential", fontsize=15)
ax.set_ylim([pot_many[-1]-200,0])
ax2 = ax.twinx()
ax2.plot(beta_many,c='r')
ax2.set_ylabel(r"$\beta$",fontsize=15)
print(pot_many[-1])



    


















    




















    
















    






-1328.47233756

















In [19]:



    


plot_nacl(traj=traj_many,left=0,right=12,num_na=4,traj_sample=-1)



    


















    




















    























In [20]:



    


system_config_mc = creator.generate_problem()



    











In [21]:



    


system_config_mc.xyz = traj[-1]



    











In [22]:



    


sampler_mc = particlesim.api.Sampler(system_config_mc)

Do plain metropolis monte carlo with last configuration from simulated annealing with three different temperatures

  • T_1 = 200 Kelvin
  • T_2 = 800 Kelvin
  • T_3 = 1200 Kelvin




In [23]:



    


beta_mc = particlesim.utils.conversion.kelvin_to_beta(200)
iteration_mc = 10000
step_mc = 0.001
traj_mc_200, pot_mc_200 = sampler_mc.metropolis(iteration_number=iteration_mc,step = step_mc, beta= beta_mc) 

beta_mc = particlesim.utils.conversion.kelvin_to_beta(800)
traj_mc_800, pot_mc_800 = sampler_mc.metropolis(iteration_number=iteration_mc,step = step_mc, beta= beta_mc) 

beta_mc = particlesim.utils.conversion.kelvin_to_beta(1200)
traj_mc_1200, pot_mc_1200 = sampler_mc.metropolis(iteration_number=iteration_mc,step = step_mc, beta= beta_mc)



    











In [122]:



    


fig, axes = plt.subplots(1, 2, figsize=(10, 4))
#fig = plt.figure()
#ax = fig.add_subplot(111)
axes[0].plot(pot_mc_200, c='r', label = r"200 $K$")
axes[0].plot(pot_mc_800, c='b', label = r"800 $K$")
axes[0].plot(pot_mc_1200,c='g', label = r"1200 $K$")
axes[0].set_xlabel(r"Iteration", fontsize=15)
axes[0].set_ylabel(r"Potential [$kcal$/$mol$]", fontsize=15)
axes[0].legend()

#ax.set_ylim([pot_mc[-1]-200,0])

#ax2 = fig.add_subplot(222)
def create_hist(pot_mc):
    hist, edges = np.histogram(pot_mc, bins=50)
    rho = hist.astype(np.float64) / np.sum(hist)
    return edges,rho
edges200, rho200 = create_hist(pot_mc_200)
edges800, rho800 = create_hist(pot_mc_800)
edges1200, rho1200 = create_hist(pot_mc_1200)

axes[1].plot(edges200[1:], rho200, c='r', label = r"200 $K$")
axes[1].plot(edges800[1:], rho800, c='b', label = r"800 $K$")
axes[1].plot(edges1200[1:], rho1200, c='g', label = r"1200 $K$")

axes[1].set_xlabel(r"$E$ [$kcal$/$mol$]", fontsize=15)
axes[1].set_ylabel(r"$\rho(E)$", fontsize=15)
axes[1].legend()
fig.tight_layout()

Comparison to Boltzmann Deviation





In [185]:



    


def boltzmann(x, a, b, x_0):
    return b * np.sqrt(2/np.pi) * (x-x_0)**2 * np.exp(-(x-x_0)**2/(2 * a**2))/a**3



    











In [177]:



    


# Fitting of parameters
popt200, pcov = curve_fit(boltzmann, edges200[1:], rho200, p0=(1, 1, edges200[0]))
popt800, pcov = curve_fit(boltzmann, edges800[1:], rho800, p0=(1, 1, edges800[0]))
popt1200, pcov = curve_fit(boltzmann, edges1200[1:], rho1200, p0=(1, 1, edges1200[0]))



    











In [186]:



    


fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(edges200[1:], rho200, c='r', lw="2", alpha=0.5, label = r"200 $K$")
ax.plot(edges800[1:], rho800, c='b', lw="2",alpha=0.5,label = r"800 $K$")
ax.plot(edges1200[1:], rho1200, c='g', lw="2", alpha=0.5,label = r"1200 $K$")

# Plotting fits
x = np.linspace(edges200[1], edges200[-1], 200)
ax.plot(x, boltzmann(x, *popt200),'k--', lw='2')

x = np.linspace(edges800[1], edges800[-1], 200)
ax.plot(x, boltzmann(x, *popt800), 'k--', lw='2')

x = np.linspace(edges1200[1], edges1200[-1], 200)
ax.plot(x, boltzmann(x, *popt1200), 'k--', lw='2')

ax.set_xlabel(r"$E$ [$kcal$/$mol$]", fontsize=15)
ax.set_ylabel(r"$\rho(E)$", fontsize=15)
ax.set_title("Comparison with Boltzmann")
ax.legend()
fig.tight_layout()



    


















    




















    























In [ ]:

Content source: maqnius/compscie-mc

Similar notebooks:

notebook.community | gallery | about