Equação de Langevin

Artur Uhlik Fröhlich e Leonardo Dasso Migotto

O objetivo deste trabalho é resolver computacionalmente a equação de Langevin utilizando o método BAOAB[LEIMKUHLER.] Serão explorados os casos de partículas individuais livres ou sujeitas a um campo potencial, estudando os efeitos da variação do coeficiente de atrito no desvio quadrático médio e na transisão de fases.

Equação de Langevin

Esta equação diferencial estocástica descreve a evolução de um sistema quando sujeito a forças do tipo determinísticas e estocásticas simultâneamente. A sua aplicação mais popular é relativa ao movimento Browniano, o movimento de uma partícula macroscópica imersa em um fluído, sujeita à força de atrito excercida pelas partículas microscópicas do fluído. Neste caso, a equação pode ser escrita como:

\frac{d \vec{v}}{d t} = - γ \vec{v} / m + ξ (t) .

Na equação acima, $γ$ é o coeficiente de atrito e $ξ (t)$ é um ruído estocástico branco, que segue o Teorema Central do Limite com média 0 e desvio padrão relacionado à temperatura, a Constante de Boltzmann, $γ$ e a massa da partícula. A partir desta expressão, é possível descobrir a relação do coeficiente de difusão do fluído e os valores envolvidos na equação:

D = \frac{2 k_{B} T}{γ m} (I)

Onde $D$ é o coeficiente de difusão do meio, $k_{B}$ é a constante de Boltzmann, $T$ é a temperatura e $m$ é a massa da partícula macroscópica. Outra relação presente no livro do Frenkel ^[1], desenvolvida teoricamente, é a do coeficiente de difusão e o deslocamento quadrático médio de uma partícula no meio:

⟨ {(r (t) - r_{0} (t))}^{2} ⟩ = 2 d D t (I I)

Método BAOAB

O método numérico escolhido para realizar a integração da equação é conhecido como BAOAB, desenvolvido por Leimkuhler e Mattews ^[2] utilizado para resolver equações diferenciais estocásticas.

Ele é baseado na solução exata para o momentum,

d \vec{p} = - γ \vec{p} d t + \sqrt{2 γ m k_{B} T} d \vec{W},

e faz o uso de um método de separação das equações entre as denominadas A, B e O, respectivamente representadas:

\vec{r} (t + \frac{Δ t}{2}) = \vec{r} (t) + \frac{Δ t}{2} \vec{p} (t + \frac{Δ t}{2}) \frac{1}{m}

\vec{p} (t + \frac{Δ t}{2}) = \vec{p} (t) + \frac{Δ t}{2} \vec{f} (t)

\vec{p^{'}} (t + \frac{Δ t}{2}) = e x p (- γ Δ t) \vec{p} (t + \frac{Δ t}{2}) + \sqrt{1 - e x p (- 2 γ Δ t)} \sqrt{m k_{B} T} \vec{G}

O $\vec{G}$ aqui representa um número aleatório Gaussiano que faz o papel do ruído estocástico.

A equação "A" realiza meio passo no tempo da distância, a "B" realiza um meio passo para o momentum e o "O" contabiliza a contribuição estocástica equação.

Essas equações podem formar vários algoritmos de integração mas o utilizado nesse trabalho será o BAOAB:

\vec{p} (t + \frac{Δ t}{2}) = \vec{p} (t) + \frac{Δ t}{2} \vec{f} (t) (1)

\vec{r} (t + \frac{Δ t}{2}) = \vec{r} (t) + \frac{Δ t}{2} \vec{p} (t + \frac{Δ t}{2}) \frac{1}{m} (2)

\vec{p^{'}} (t + \frac{Δ t}{2}) = e x p (- γ Δ t) \vec{p} (t + \frac{Δ t}{2}) + \sqrt{1 - e x p (- 2 γ Δ t)} \sqrt{m k_{B} T} \vec{G} (3)

\vec{r} (t + Δ t) = \vec{r} (t + \frac{Δ t}{2}) + \frac{Δ t}{2} \vec{p^{'}} (t + \frac{Δ t}{2}) \frac{1}{m} (4)

\vec{p} (t + Δ t) = \vec{p^{'}} (t + \frac{Δ t}{2}) + \frac{Δ t}{2} \vec{f} (t) (5)

É importante lembrar que entre os dois últimos passos é necessário atualizar o termo $\vec{f}$ , já que ele pode depender de termos já atualizados como $\vec{p}$ ou $\vec{r}$ .

Implementação

Um exemplo de implementação desse método feito foi para a partícula livre. Nota-se que nesse caso a partícula não é afetada por um campo potencial então os passos do método que envolvem a força $\vec{f}$ serão desconsiderados.

A função a seguir se encontra em um código (orientado a objeto) de autoria de do integrante Artur e que pode ser acessado aqui.

    def baoab_livre(self, dt, exp, sqexp, sqt, G):
        '''
        dt: discretização do tempo
        exp: termo referente a primeira exponencial da eq.3
        sqexp: termo da raiz quadrada com exponencial da eq.3
        sqt: termo da raiz quadrada com a Temperatura da eq.3
        G: o vetor G da eq.3
        '''
        # 1/2 passo da distância
        self.pos[0] += self.vel[0]*(dt/2)
        self.pos[1] += self.vel[1]*(dt/2)
        # Passo estocástico
        self.vel[0] = exp*self.vel[0] + sqexp*sqt*G[0]
        self.vel[1] = exp*self.vel[1] + sqexp*sqt*G[1]
        # Atualização final da posição
        self.pos[0] += self.vel[0]*(dt/2)
        self.pos[1] += self.vel[1]*(dt/2)

Essa função representa um passo do laço temporal utilizando o método BAOAB. Os atributos em forma de lista ".pos" e ".vel" representam as componentes x e y respectivamente da posição e velocidade do objeto partícula.

Um exemplo de funcionamento do código é a animação a seguir:

. O código foi executado usando a semente de números aleatórios 420.

Análise do MSD

Os resultados a serem analisados são originados do código (procedural) do integrante Leonardo.

Foram executadas diferentes vezes a simulação para compararmos a influência do $γ$ e da temperatura no deslocamento quadrático médio.

Temperaturas de 1-3 e $γ = 0.1$ :

Temperaturas de 1-3 e $γ = 1$ :

Temperaturas de 1-3 e $γ = 10$ :

De início pode-se notar que o $Δ r^{2}$ escala de maneira proporcional a $t^{2}$ (balístico) no início do processo, e depois de um determinado tempo obtemos uma relação linear que é a parte difusiva do processo. Essas inclinações das retas nos regimes balístico e difusivo são constantes em relação ao $γ$ e à temperatura. Outro aspecto constante é a temperatura, que não gera aparente diferença nos gráficos de MSD

Uma referência da literatura a respeito do MSD (mas no contexto de dinâmica molecular usando o potencial de Lennard-Jones) se encontra no livro do Frenkel, expressa por essa imagem:

A partir da figura pode-se notar que a inclinação das retas representantes dos regimes balístico e difusivo são 2 e 1, valores esses encontrados pelas simulações realizadas.

Outro aspecto importante a se ressaltar é a concordância do cálculo do coeficiente de difusão através da equação (I) e da relação de Einstein em (II), obtida através dos pontos presentes na parte difusiva.

Uma relação interessante encontrada foi a maneira como escala a mudança de regimes, indicada pela reta perpendicular ao eixo do tempo. Podemos encontrar uma aproximação para o tempo em que o regime difusivo inicia na simulação a partir da relação dependente de $γ$ :

Ω = \frac{10}{γ}

Referências

↑ Daan Frenkel and Berend Smit. 2001. Understanding Molecular Simulation (2nd. ed.). Academic Press, Inc., USA.
↑ Leimkuhler, B., & Matthews, C. (2015). Molecular Dynamics: With Deterministic and Stochastic Numerical Methods. (Interdisciplinary Applied Mathematics; Vol. 39). Springer. https://doi.org/10.1007/978-3-319-16375-8

[1] Daan Frenkel and Berend Smit. 2001. Understanding Molecular Simulation (2nd. ed.). Academic Press, Inc., USA.

[2] Leimkuhler, B., & Matthews, C. (2015). Molecular Dynamics: With Deterministic and Stochastic Numerical Methods. (Interdisciplinary Applied Mathematics; Vol. 39). Springer. https://doi.org/10.1007/978-3-319-16375-8

[1]

[2]

Equação de Langevin

Índice