Potencial de Lennard-Jones

Nome: André Bracht Burmeister

Fluidos estão em todos os lugares e, portanto, estudá-los é muito importante. Como são objetos de estudo extremamente complicados, se fazem necessários modelos e simplificações.

Uma possível simplificação é a análise do sistema em duas dimensões apenas, excluindo a terceira dimensão espacial. A segunda simplificação da simulação

O Modelo

O Potencial de Lennard Jones

O potencial $U'$ que uma partícula faz na outra, pode ser aproximado a uma distribuição em termos da distancia $r'$ entre elas:

$U'(r')=4\epsilon \left[\left({\frac {\sigma }{r'}}\right)^{12}-\left({\frac {\sigma }{r'}}\right)^{6}\right]$

Código

"""lennard_jones.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1jSPjBz0RXXWLEAzf3zcHgdlQpzjyO_ni

# Potencial de Lennard-Jones
# Monte Carlo no contínuo

#### Delete variables and import packages
"""

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation 
import random
from scipy.optimize import bisect as find_root

"""## Simulação:

### Definições Numéricas:
"""

densidade = 0.1 #(2 ** (2/3)) * (3 ** (-1/2))  # partícula(s) por sigma^2 (partículas por area)
densidade_round = round(densidade, 2)
n_particles = 64
L = np.sqrt(n_particles / densidade)  # unidade de sigma
print(L)
dLmax = 0.1
MCS = 10000
eps = 1
T = .1
k = 1
beta = 1/(k*T)
raio = 2**(-5/6)

"""### Funcões:"""

def lennard_jones(r1, r2, L):  # r1 e r2 em unidades de sigma (na vdd seria r1/sig e r2/sig)
    dr = np.absolute(np.subtract(r2, r1))
    for i in np.arange(2):
        if dr[i] > 0.5:
            dr[i] = L - dr[i]

    r2 = (dr[0] ** 2) + (dr[1] ** 2)
    sr2 = 1/r2
    return 4 * (((sr2) ** 6) - ((sr2) ** 3))


def initial_energy(pos_particles, n_particles, L):
    U0 = 0
    for xxx in range(n_particles - 1):
        for yyy in range(xxx+1, n_particles):
            U0 += lennard_jones(pos_particles[:, xxx], pos_particles[:, yyy], L=L)
    return U0


def energy_var(p, pos_particles, n_particles, L, new_pos=np.zeros(2)):
    U1 = 0
    U2 = 0
    for i in range(n_particles - 1):
        j = (i + p + 1) % n_particles  # j vai ser todas as partículas 
        # exceto a partícula p (porque p interage com todas, menos com ela mesma)
        U1 += lennard_jones(pos_particles[:, j], pos_particles[:, p], L=L)
        U2 += lennard_jones(pos_particles[:, j], new_pos, L=L)
    return U2 - U1


def main_loop(dLmax, beta, n_particles, MCS, L, acceptance_array=False, 
              last_acceptance=False):
    
    pos_particles = L * np.random.random([2, n_particles])
    # Random positions:
    pos_MCS = np.zeros([MCS+1, 2, n_particles])
    pos_particles = pos_particles
    pos_MCS[0] = pos_particles

    # Energy:
    energy_arr = np.zeros(MCS+1)
    energy = initial_energy(pos_particles, n_particles, L=L)
    energy_arr[0] = energy

    # Acceptance:
    accepted = 0
    if acceptance_array:
        acceptance_arr = np.zeros(MCS)

    # Main Loop:
    for m in range(MCS):
        for _ in range(n_particles):
            p = random.randint(0, n_particles - 1)
            # change position:
            displacement = dLmax * ((2 * np.random.random([2])) - 1)
            new_pos = pos_particles[:, p] + displacement
            for i in range(2):
                if new_pos[i] < 0:
                    new_pos[i] += L
                elif new_pos[i] > L:
                    new_pos[i] -= L 
            dU = energy_var(p, pos_particles, n_particles, L, new_pos)
            if dU <= 0:
                pos_particles[0][p] = new_pos[0]
                pos_particles[1][p] = new_pos[1]
                energy += dU
                accepted += 1
            else:
                if random.random() < np.exp(-beta * dU):
                    pos_particles[0][p] = new_pos[0]
                    pos_particles[1][p] = new_pos[1]
                    energy += dU
                    accepted += 1
        if not (acceptance_array):
            pos_MCS[m+1] = pos_particles
            energy_arr[m+1] = energy

        if acceptance_array:
            acceptance_arr[m] = accepted / (n_particles*(m + 1))
    if acceptance_array:
        return acceptance_arr
    elif last_acceptance:
        return accepted / (n_particles*MCS)
    else:
        return pos_MCS, energy_arr, (accepted / (n_particles*(m + 1)))


def save_data(dLmax=dLmax, n_particles=n_particles, MCS=MCS, T=T, densidade=densidade):
    beta = 1 / (k*T)
    L = np.sqrt(n_particles / densidade)  # unidade de sigma
    print(L)
    pos, energy, acceptance = main_loop(dLmax=dLmax, beta=beta, n_particles=n_particles, 
                                        MCS=MCS, L=L)
    
    densidade_round = round(densidade, 3)
    # save positions:
    pos_reshaped = pos.reshape(pos.shape[0], -1)
    np.savetxt(f"positions_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}.txt", 
               pos_reshaped)

    # save energy:
    np.savetxt(f"energy_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}.txt", 
               energy)
    return pos, energy, acceptance


def open_pos(file, n_particles):
    loaded_pos = np.loadtxt(file)
    load_original_pos = loaded_pos.reshape(loaded_pos.shape[0], loaded_pos.shape[1] // n_particles, n_particles)
    return load_original_pos


"""### Rodar e salvar dados:

#### Rodar uma vez:
"""

dLmax = round(0.5*raio, 3)
pos, energy, acceptance = save_data(dLmax=dLmax)
print(f'dLmax={dLmax} \n beta={beta} \n n_particles={n_particles} \n MCS={MCS} \n\n pos[0]: \n {pos[0]}')
print(acceptance)

"""#### Rodar para vários T e rô:"""

dLmax = round(0.5*raio, 3)
n_particles = 32 
T_arr = np.array([0.1, 1.5, 3])
dens_arr = np.array([0.1, 0.5, 0.9])
for T in T_arr:
    print()
    print(f'T = {T}')
    for densidade in dens_arr:
        print(f'd = {densidade}')
        pos, energy, acceptance = save_data(dLmax=dLmax, n_particles=n_particles, MCS=MCS, T=T, densidade=densidade)
        #print(f'dLmax={dLmax} \n beta={beta} \n n_particles={n_particles} \n MCS={MCS} \n\n pos[0]: \n {pos[0]}')
        print(f'acceptance = {acceptance}')


"""
## Testes:

### Encontrar deslocamento máximo:
Testando valores:
"""

def testing_dL(possible_dL, times):
    for l in possible_dL:
        accept = np.zeros(times)
        for d in range(times):
            accept[d] = main_loop(dLmax=l*L, last_acceptance=True, acceptance_array=False)
        print(f'dLmax = {l}*L, acceptance = {np.mean(accept)}')

# Roda mais geral:
possible_dL = np.round(0.1*np.arange(1, 10), 1)
testing_dL(possible_dL, 5)

# Ajuste fino:
possible_dL = np.round(0.01*np.arange(10, 14), 3)
testing_dL(possible_dL, 50)

# Novo dLmax:
dLmax = 0.12*L
print(dLmax)

"""### Comportamento da acceptance:"""

def acceptance_behaviour(possible_dL, times):
    for d in range(times):
        fig, axes = plt.subplots(1, len(possible_dL), figsize=(20, 5))
        fig.suptitle(f'T = {T}   N = {n_particles}   d = {densidade_round}    dLmax = {round(dLmax, 2)}', fontsize=14)
        for l in range(len(possible_dL)):
            print(f'{l+1}/{len(possible_dL)}, {d+1}/{times}')
            acceptance = main_loop(dLmax=possible_dL[l]*raio, acceptance_array=True)
            axes[l].plot(np.arange(len(acceptance)), acceptance)
            axes[l].set_ylim([0, 1])
            axes[l].set_title(f'dLmax = {possible_dL[l]}*2^(-5/6)')
        plt.show()
        plt.close()

# Plotando acceptance
possible_dL = np.round(0.2*np.arange(1, 5), 1)
acceptance_behaviour(possible_dL, 1)

# Mais fino:
possible_dL = np.round(np.concatenate([0.0001*np.arange(1, 6)]), 4)
acceptance_behaviour(possible_dL, 1)


"""
## Figuras

### Intial and final positions:
"""

def scatter_pos(pos, indices, T, L, circles=True):
    graphs = len(indices)
    fig, axes = plt.subplots(1,graphs)
    fig.suptitle(f'T = {T}   N = {n_particles}   d = {densidade_round}    dLmax = {round(dLmax, 2)}', fontsize=14)

    for j in range(graphs):
        k = indices[j]
        if circles:
            # draw the circles where they are
            for i in range(n_particles):
                c = plt.Circle([pos[k][0][i], pos[k][1][i]], raio, fill=False, edgecolor='black', alpha = 0.5)
                axes[j].add_patch(c)
            # draw the circle's projections:
            for a in range(-1, 2):
                for b in range(-1, 2):
                    for i in range(n_particles):
                        c = plt.Circle([pos[k][0][i] + a*L, pos[k][1][i] + b*L], raio, fill=False, edgecolor='black', alpha=0.25)
                        axes[j].add_patch(c)
        
        points = axes[j].scatter(pos[k][0], pos[k][1])
        margem = raio
        axes[j].plot([0, L, L, 0, 0], [0, 0, L, L, 0], color='black')
        axes[j].axis([-margem, L + margem, -margem, L + margem], 'equal')
        axes[j].grid(visible=False)
        axes[j].set_aspect("equal")
        axes[j].set_title(f'passo {k}')
    return fig, axes, points


for T in T_arr:
    print()
    print(f'T = {T}')
    for densidade in dens_arr:
        print(f'd = {densidade}')
        densidade_round = round(densidade, 3)
        file_end = f"_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}"
        file_name = f"positions{file_end}.txt"
        pos = open_pos(file_name, n_particles=32)

        L = np.sqrt(n_particles / densidade)  # unidade de sigma
        plot = scatter_pos(pos=pos, indices=[0, pos.shape[0] - 1], T=T, L=L, circles=True)
        plt.show()
        plt.savefig(f"first_last_positions{file_end}.png")
        plt.close()

"""### Energy:"""

def y_limits(energy):
    energy_hist_range = energy[int(round(0.1*len(energy), 0)):]
    return [min(energy_hist_range)-1 , max([0, max(energy_hist_range) + 1])]

"""#### Energy graph:"""

for T in T_arr:
    print()
    print(f'T = {T}')
    for densidade in dens_arr:
        print(f'd = {densidade}')
        file_end = f"_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}"
        energy = np.loadtxt(f"energy{file_end}.txt")

        # Plot:
        fig, ax = plt.subplots(1, 1)
        titulo = f'T = {T}   N = {n_particles}   d = {densidade_round}    L = {round(L, 2)} \n dLmax = {round(dLmax, 2)}    acceptance = {round(acceptance, 2)}'
        fig.suptitle(titulo,
                    fontsize=14)
        energy_range = energy
        ax.plot(np.arange(len(energy_range)), energy_range)
        limites = y_limits(energy)
        print(limites)
        ax.set_ylim(limites)
        plt.show()
        plt.savefig(f"energy_graph{file_end}.png")
        plt.close()


        # Histogram:
        fig, ax = plt.subplots(1, 1)
        fig.suptitle(titulo,
                    fontsize=14)
        ax.hist(energy[100: ], density=True)
        ax.set_xlabel('energy')
        plt.show()
        plt.savefig(f"energy_histogram{file_end}.png")
        plt.close()

"""#### Energy histogram:"""



"""## Animação"""

file_name = f"positions_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}.txt"
print(file_name)

start = 0
end = len(pos)
pos = open_pos(file_name)[start:end]
fig, ax = plt.subplots(1,1, figsize=(5, 6))
fig.suptitle(f'T = {T}   N = {n_particles}   d = {densidade_round}    L = {round(L, 2)} \n dLmax = {round(dLmax, 2)}    acceptance = {acceptance}',
             fontsize=14)
margem = raio
ax.axis([-margem, L + margem, -margem, L + margem], 'equal')
ax.grid()
ax.set_aspect("equal")
ax.set_title(f'passo {k}')

c = []
for i in range(n_particles):
    c.append(plt.Circle([pos[0][0][i], pos[0][1][i]], raio, fill=False, color='lightgrey'))
    ax.add_patch(c[i])
scatter_points = ax.scatter(pos[0][0], pos[0][1], s=1000)


def animate(i):
    for k in range(n_particles):
        c[k].center = (pos[i][0][k], pos[i][1][k])
    scatter_points.set_offsets(pos[i].T)
    return scatter_points


anim = FuncAnimation(fig, animate, interval = 10, frames = end-start-1)

anim_file = f'animacao_T{T}_N{n_particles}_d{densidade_round}_MCS{MCS}_dL_{dLmax}.mp4'
print(anim_file)
anim.save(anim_file, writer = 'ffmpeg', fps = 20)

"""## Imagens:

### Raio:
"""

x0 = 1
y0 = 1
l = 2 * (2 ** 1/6)
triX = [x0, x0+l, x0+(l/2)]
h = np.sqrt(3)*l/2
triY = [y0, y0, y0 + h]

triX.append(triX[0])
triY.append(triY[0])

fig, ax = plt.subplots()

for i in range(3):
    c = plt.Circle([triX[i], triY[i]], l, alpha=0.5)
    c = plt.Circle([triX[i], triY[i]], l, alpha=0.5)
    c = plt.Circle([triX[i], triY[i]], l, alpha=0.5)
    ax.add_patch(c)

ax.plot(triX, triY, color = 'b')
ax.axis([x0-(l/2), x0 + (3*l/2), y0 - (l/2), y0 + h + l/2], 'equal')
ax.grid()
ax.set_aspect("equal")
plt.show()
plt.close()

"""### Potencial Lennard Jones:"""

def lennard_jones1d(r):
    sr6 = 1/(r ** 6)
    return 4 * (((sr6) ** 2) - (sr6))


xx = np.linspace(0, 4, 100)
yy = lennard_jones1d(xx)
fig, ax = plt.subplots(1, 1, figsize=(10, 4))
ax.plot([-10, 10], [0, 0], color='black')
ax.plot(xx, yy)
ax.set_ylabel('U(r)')
ax.set_xlabel('r')
ax.set_ylim([-1.5, 5])
ax.set_xlim([0, 3])
plt.show()
fig.savefig('lennard_jones_potential.png')
plt.close()