Introdução à dinâmica molecular clássica

só um divisor

O problema de Fermi Pasta e Ulam

A importância do problema de Fermi Pasta e Ulam (ou FPU) está além do resultado Físico (que é sumamente importante, como veremos). Antes deste problema ser atacado os computadores eram utilizados como mera ferramente de cálculo, antes do FPU você tinha equações diferenciais e um método de resolução numérico os quais eram programado no computador e obtinha um resultado. A partir do FPU passou-se a analisar os dados que o computador proporcionava de foma quase similar a como um físico experimental trabalha, a isso chamamos de uma simulação numérica. Fermi Pasta e Ulam utilizaram o computador (o MANIAC-I: Mathematical Analyzer, Numerical Integrator and Computer) para reproduzir um fenômeno físico supondo válida as equações matemáticas utilizadas para descrever o problema. Ainda o computador continuava a fazer exatamente a mesma coisa antes e depois do FPU, a forma como eram encarado os resultados foi o que realmente mudou.

A ideia do trabalho de Enrico Fermi, John Pasta, and Stanislaw Ulam (1955) foi estudar como um cristal evoluí para o equilíbrio simulando uma cadeia de partículas de massa unitária que interagem via um potencial harmônico mais uma interação não linear fraca. O hamiltoniano que descreve esse sistema está dado por \[ H = \sum_{i = 0}^{N-1}\frac{1}{2}p_i^2 + \sum_{i = 0}^{N-1}\frac{1}{2}\left( u_{i+1} - u_i \right)^2 + \frac{\alpha}{2}\sum_{i = 0}^{N-1}\frac{1}{2}\left( u_{i+1} - u_i \right)^3 \] onde $u_i$ é o deslocamento da partícula $i$ da sua posição de equilíbrio. Observe que as partículas nos extremos estão fixas. A partir deste hamiltoniano obtemos a equação de movimento do sistema dada por \[ \ddot{u} = \left( u_{i+1} - u_{i-1} + 2u_i \right) + \alpha\left[\left( u_{i+1} - u_i \right)^2 - \left( u_i - u_{i-1} \right)^2 \right] \label{EqMov} \]

A fim de realizar sua analise, eles colocaram o sistema no seu $k-esimo$ seu modos normal dado por: \[ u_i = A\sqrt{\frac{2}{N}}\sin\left( \frac{ik\pi}{N} \right) \label{UmModoNomal} \] Como o teorema de equipartição afirma que cada termo quadrático do hamiltoniano contribui com $k_BT$ para a energia do sistema então, ao analisar o problema considerando o sistema descrito em termo dos modos normais de vibração no limite em que $\alpha = 0$, o Hamiltoniano passa a ter a forma \[ H = \frac{1}{2}\sum_{k = 0}^{N-1}\left(\dot{a}_k^2 + \omega_k a_k^2\right) \label{hnormal} \] sendo que $a_k$ são as amplitudes de cada modo de vibração, dada por \[ a_k = \sqrt{\frac{2}{N}}\sum_{k=0}^{N-1}u_i\sin\left( \frac{ik\pi}{N} \right) \label{transFourier} \] e $\omega_k$ está dada pela relação de dispersão: \[ \omega_k = 2\sin\left(\frac{k\pi}{N}\right) \label{dispersao} \] de onde fica claro (eq $\ref{hnormal}$) que cada modo contribui com $k_BT$ para a energia e essa energia é própria ao modo e não sendo compartilhada com outros modos.

Sabendo isso eles modificaram o valor da interação no linear com a esperança de distribuir a energia entre os diversos modos do sistema (anos antes de fazer este experimento numérico, Fermi tinha demonstrado que utilizando uma pequena perturbação não linear o sistema deveria se tornar ergódico) mas encontraram que somente depois de ultrapassado um valor crítico de acoplamento no linear é que se obtém a esperada distribuição energética, abaixo desse valor a energia fica restrita ao modos vizinhos mas depois de um tempo volta toda para o modo original. Como esse comportamento estava em desacordo com o esperado da teoria ergódica o trabalho ficou conhecido como paradoxo de Fermi, Pasta e Ulam ou problema FPU.

O algoritmo de Verlet

O nosso problema é reproduzir o trabalho original de Fermi, Pasta e Ulam e para isso utilizaremos um algoritmo conhecido como algoritmo de Verlet. Sem demostração (futuramente retornaremos à discussão deste algoritmo), aceitaremos que o algoritmo de Verlet tem a seguente forma \[ \vec{r}_i(t + \delta t) = \vec{r}_i(t) - \vec{r}_i(t - \delta t) + \vec{a}_i(t)\delta t \label{verletR} \] e, a partir da equação anterior \[ \vec{v}_i(t) = \frac{\vec{r}_i(t + \delta t) - \vec{r}_i(t - \delta t)}{2\delta t} \label{verletV} \] este é um algoritmo onde o erro é de ordem $O(4)$ igual que o Runge-Kutta. Perceba que estamos trabalhando com vetores, por tanto teremos equações similares às de acima para cada componente.

Na figura acima se mostra esquematicamente o funcionamento do algoritmo:

  1. Ter armazenado $\vec{r}(t-\delta t)$ e $\vec{r}(t)$, utilizando este último calcular $\vec{a}(t)$. No primeiro passo de integração é necessário calcular $\vec{r}(t-\delta t)$ utilizando uma expansão em serie de Taylor até II ordem com o $\delta t$ negativo: \[ \vec{r}_i(t - \delta t) = \vec{r}_i(t) - \vec{v}_i(t)\, \delta t \nonumber \]
  2. Utiliza-se a equação $\ref{verletR}$ para calcular $\vec{r}_i(t + \delta t)$
  3. Finalmente fazemos \[ \left\{ \begin{align*} \vec{r}(t-\delta t) =& \vec{r}(t)\\ \vec{r}(t) =& \vec{r}(t+\delta t) \end{align*} \right. \] estado assim prontos para repetir o passo inicial.

Continuando com FPU

Para o caso do problema de FPU nossa aceleração está dada pela equação $\ref{EqMov}$. Como pode ser lido no trabalho original de FPU devemos monitorar a energia de alguns modos de vibração, principalmente aqueles adjacentes ao modo inicialmente exitado, para isso relembramos que utilizando uma serie de Fourier, no caso discreta, podemos expressar qualquer curva, em especial estas que são harmônicas: \[ u_i = \sqrt{\frac{2}{N}}\sum_{k=0}^{N-1}a_k\sin\left( \frac{ik\pi}{N} \right) \] Note está serie não é mais do que uma superposição dos diferentes modos normais que o sistema pode adotar. Do ponto de vista mecânico temos uma forma de expressar as coordenadas cartesianas a partir das coordenadas dos modos normais. Se utilizamos as coordenadas normais, o hamiltoniano fica escrito na forma indicada pela equação $\ref{hnormal}$ onde os $a_k$ são as amplitudes dadas por $\ref{transFourier}$. Sabemos que o sistema descrito pelo hamiltoniano $\ref{hnormal}$ é conservativo ($\alpha = 0$) então ele é a energia do sistema, portanto podemos, a partir dele, calcular a energia de cada modo: \[ E_k = \frac{1}{2}\left[ \left( \frac{da_k}{dt} \right)^2 + \omega_k a_k^2 \right] \]

Caso consideremos o acoplamento no linear, o hamiltoniano $\ref{hnormal}$ passa a ser escrito como \[ H = \frac{1}{2}\sum_{k = 0}^{N-1}\left(\dot{a}_k^2 + \omega_k a_k^2\right) + \frac{\alpha}{3}\sum_{k, l, m = 0}^{N-1} c_{klm}a_ka_la_m \omega_k \omega_l \omega_m \] onde os coeficiente $c_{klm}$ são constantes e podem ser encontradas no trabalho de D. Sholl

Do FPU aos Solitons e o Caos

Porque Fermi não pode demostrar sua hipótese? Alguns anos depois do trabalho de FPU, Kolmogorov propôs um teorema que foi demostrado por Moser e Arnol (teorema KAM) o qual dizia que para hamiltonianos integráveis (com solução) algumas orbitas que são ligeiramente perturbadas se mantem quase periódicas dessa forma isto explica qualitativamente o porque o sistema de FPU apresenta as recorrências observadas mas não da uma resposta definitiva ao problemas.

Uma melhor explicação explicação foi dada por Korteweg e de Vries, eles partiram da equação de movimento onde consideraram explicitamente a consta de mola $k$ e a massa: \[ m\ddot{u} = k\left( u_{i+1} + u_{i-1} - 2u_i\right)\left[ 1 + \alpha \left( u_{i+1} - u_{i-1} \right) \right] \nonumber \] mudando de variáveis: \[ k \rightarrow \frac{\kappa}{\Delta u}\;\;\;\;m \rightarrow \rho\Delta u\;\;\;\; \frac{\rho}{\kappa} \rightarrow c \equiv 1\;\;\;\; u_i(t) \rightarrow x_u(u,t) \nonumber \] chegamos à equação \[ \ddot{x}(u,t) = \left( 1 + \alpha x_u \right)x_{uu} \nonumber \] Como essa equação apresentou problemas de multiplas solução antes do tempo de recorrência, Kruskal e Zabusky introduziram um termo dispersivo para melhorar a solução \[ \ddot{x}(u,t) = x_{uu} + \alpha x_ux_{uu} + \beta x_{uuuu} \nonumber \] que trabalhando com ondas propagantes somente em uma direção e considerando as transformações \[ u \rightarrow u - t\;\;\;\;\; t \rightarrow t/\alpha\;\;\;\;\; x_u \rightarrow x \nonumber \] obtém-se as equações \[ x_{t} + xx_{u} + \delta x_{uuu} = 0\;\;\;\;\; x_t - 6xx_u - x_{uuu} = 0 \nonumber \] com $\delta = \beta / \alpha$. Está equação, conhecida como equação kdV, possui soluções de ondas solitárias que se propagam sem modificação de suas formas (solitons), solução está primeira vez observada por Diederik Korteweg e Gustav de Vries estudando ondas num canal raso com líquidos de baixa viscosidade. A solução formal para a equação anterior é \[ x(u,t) = \frac{c}{2}\text{sech} ^2 \left[ \frac{\sqrt{c}}{2} \left( x - ct + \delta \right) \right] \] onde $c$ e $\delta$ são constantes. Kruskal e Zabusky observaram que os sólitons mantinham sua forma após interagirem. A partir destes resultados uma interpretação possível ao problema de FPU é que os sólitons são equivalentes ao "fonôns clássicos" da rede não linear dessa forma a energia fica "distribuída" nos diversos sólitons da rede e como ele não se perturbam após interagirem então é fácil entender o porque essa energia não é compartilhada com outro modo solitônico.

Foi demostrado depois que esse tipo de equação possui infinitas constantes de movimento e consequentemente os sistemas se lhes denomina de integráveis.

Tarefa

Construa um programa que simule o FPU utilizando $N=33$, exitando inicialmente o modo $k=1$ e colocando $\alpha = 0$. Utilize $\delta t = 0.05$ e $a_1 = 10$. Faça o gráfico da energia dos modos 1, 3, 5. Modifique o programa para a situação em que $\alpha = 0.3,\; 1\; 3$ e além de plotar a energia dos modos anteriores plote $\dot{a}_1$ vs $a_1$ para cada caso.