Medidas dinâmicas

Em Dinâmica Molecular são chamadas de Medidas (ou Propriedades) Dinâmicas aquelas características do sistema que variam no tempo. Aqui, apresentaremos duas medidas dinâmicas de interessante importância: o deslocamento quadrático médio e a autocorrelação de velocidades.

Ambas essas medidas estão relacionadas com o coeficiente de transporte chamado difusão: a quantificação do modo em que uma substância se espalha em um meio. Particularmente, em nosso caso, as moléculas simuladas se difundem (espalham) no próprio meio, o que chamamos de autodifusão. A constante de difusão tem unidades SI de $[D]=m^{2}/s$. Em unidades reduzidas para o caso de um potencial Lennard-Jones, a constante de difusão tem unidades $[D]={\sqrt {\frac {\epsilon }{m}}}\sigma$. Talvez estas unidades fiquem mais claras após verificar o método de calculo deste coeficiente.

Outras propriedades dinâmicas podem ser citadas, como a condutividade térmica, condutividade elétrica e viscosidade.

As deduções das equações apresentadas fogem do escopo desta página. Todas elas, porém, podem ser encontradas no livro Understanding Molecular Simulation, de D. Frankel e B. Smit, páginas 87-90 (ref. [1]). As deduções partem da Lei de Fick, que estipula que a difusão de uma substância é proporcional à variação da concentração dela: ${\vec {j}}=-D\nabla c$, o sinal de menos indica que o fluxo de difusão é no sentido para onde a concentração dimunui.

Deslocamento quadrático médio

Uma introdução teórica

É possível acompanhar o movimento das partículas vermelhas nessa dinâmica. Como quantificar o quanto elas se movem em um intervalo de tempo? Um dos métodos possíveis é determinar o Deslocamento Quadrático Médio, uma função do tempo.

O Deslocamento quadrático médio, do inglês Mean square displacement (MSD), é uma medida da capacidade de dispersão (difusão) das partículas no sistema. Ela representa a média dos quadrados dos deslocamentos em relação à uma posição anterior de referência de cada partícula (Note que a média dos deslocamentos das partículas é nula, pois o centro de massa não se move, por isso a necessidade de tomarmos os quadrados). Matematicamente, podemos expressar:

${\rm {MSD}}\equiv \langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle ={\frac {1}{N}}\sum _{i=1}^{N}(x_{i}(\Delta t)-x_{0i})^{2}+(y_{i}(\Delta t)-y_{0i})^{2}$

Perceba que o MSD é uma função do tempo (na verdade, do intervalo de tempo) e portanto estamos interessados em estudar como esta função varia no tempo.

Em dinâmicas com fases líquidas e gasosas, teremos então o fenômeno da dispersão. As partículas irão "caminhar" por todo o espaço permitido em um movimento aleatório. Este movimento é conhecido como movimento Browniano e fui estudando por grandes cientistas, como Albert Einstein. Einstein conseguiu modelar o deslocamento quadrático médio como uma função linear do tempo, e determinou uma constante chamada de constante de difusão (D):

$\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle =2dD\Delta t$

Exemplos dos regimes de difusão que podem ser compreendidos com o MSD em função do tempo reduzido.

Na equação acima d é o número de dimensões simuladas ($d=2$, em nosso caso bidimensional). Assim, este é um método amplamente usado para quantificar a difusão do sistema. Qualitativamente, classifica-se os sistemas enquanto ao MSD da seguinte forma:

• Regime confinado: $\lim _{\Delta t\to \infty }\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle =\langle r_{c}^{2}\rangle$

• Regime subdifusivo: $\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle \propto \Delta t^{\alpha },0<\alpha <1$

• Regime normalmente difusivo: $\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle =2D\Delta t^{1}$

• Regime superdifusivo: $\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle \propto \Delta t^{\alpha },\alpha >1$

• Regime balístico: $\langle ({\vec {r}}-{\vec {r_{0}}})^{2}\rangle \propto \Delta t^{2}$

Em uma dinâmica com um estado sólido, esperamos que as partículas apenas vibrem, e não se deslocam grandes distâncias. Portanto, o MSD terá um limite e diremos que o regime é confinado, pois as partículas não têm espaço para se difundirem. Já o regime balístico ocorre para pequenos intervalos de tempo em simulações de fases líquidas ou gasosas, quando inicialmente, além da difusão, há também um impulso unidirecional.

A imagem ao lado mostra o MSD em função do tempo para cada um destes casos qualitativos.

Implementação computacional

Primeiro é importante ressaltar que precisamos comprar a posição da partícula com a posição da mesma partícula em um tempo anterior de referência, sem levar em consideração o condições de contorno periódicas. Nos trechos de códigos apresentados a seguir, consideraremos que as posições não são corrigidas quanto à este aspecto. Caso forem corrigidas, é necessário manter para cada partícula uma variável que represente o número de vezes que ela atravessou a parede em um sentido.

Definiremos primeiro um período de passos em que analisaremos o MSD.

#define numMSD 1000 // Período do MSD

Alocaremos dinamicamente na memória um vetor double com o tamanho de passos máximo. Também será necessário vetores para armazenarmos as posições de referência de cada partícula (NP = Número de partículas).

double *x_zero  = (double *) malloc( NP * sizeof(double) ); //Valores de x de referência
double *y_zero  = (double *) malloc( NP * sizeof(double) ); //Valores de y de referência
double *msd = (double *) malloc( numMSD * sizeof(double) ); //Valores de deslocamento MSD(t)

Declararemos também uma assinatura da função MSDf. Para ela, será necessário enviar as posições atuais, as posições de referência, o vetor onde armazenamos os valores de deslocamento e também a quantidade de passos a frente da posição de referência que estamos.

void MSDf (double *x, double *y, double *x_zero, double *y_zero, double *msd, int ref);

No corpo principal do programa, chamaremos a função MSDf. Tomaremos o cuidado de garantir que calcularemos o MSD após atingirmos um ponto de equilíbrio do sistema (assim, ignoraremos os primeiros numMSD passos com um if ).

Manteremos também um counter, para sabermos quantas vezes passamos por esta função e acumulamos dados no vetor msd. Consideramos que o counter é igual para todos os valores de $\Delta t$. Caso isso não seja verdade, é necessário um vetor que contabiliza quantos vezes calculou-se o MSD para cada intervalo.

if( t >= numMSD*dt ){
 		MSDf(x, y, x_zero, y_zero, msd, passo % numMSD);
		counter++;
	}

Note que o valor de passo%numMSD varia de 0, quanto tomaremos nossos valores de referência, até numMSD-1 que será o nosso maior intervalo de passos. No caso deste programa calculamos o MSD para cada $\Delta t$. Geralmente programas que não apresentam essa necessidade, o que pode ser facilmente ajustado.

Então, teremos o corpo da função MSDf, que, caso o parâmetro ref é zero, toma as posições atuais como referência. Se não, calcula a diferença entre a posição atual e a posição de referência, fazendo a média sobre essa informação das partículas e salvando no vetor msd.

void MSDf (double *x, double *y, double *x_zero, double *y_zero, double *msd, int ref){

int i;

if(ref == 0){ //Se estamos no tempo de referência, salvamos a posição atual
  for(i=0; i<NP; i++){
    x_zero[i]=x[i];
    y_zero[i]=y[i];
	}
}else{
double aux=0; // Essa variável irá acumular a soma dos deslocamentos quadráticos

for(i=0; i<NP; i++)
	aux += (x[i]-x_zero[i])*(x[i]-x_zero[i]) + (y[i]-y_zero[i])*(y[i]-y_zero[i]);

msd[ ref ] += aux/NP; // Então acumularemos o deslocamento quadrático médio
                      // por partícula para um intervalo ref de passos no vetor msd.
}

}

Por fim, depois do laço temporal, dividiremos o valor de counter por numMSD, dando assim o valor de vezes em que se acumulou o valor do MSD por intervalo de passos. Para imprimirmos, usaremos o intervalo de passo em intervalo de tempo reduzido, multiplicando o número de passos pelo tamanho de cada passo.

counter = counter/numMSD; //Counter marcará quantas vezes se passou pela função para cada 
                          //valor de intervalo de passos ref.

for(i=1; i<numMSD; i++){
	fprintf(arq, "%lf\t%lf\n", i*dt, (double)msd[i]/counter); //Impressão

No arquivo apontado por arq, constará ao final do programa:

0.010000	0.000145
0.020000	0.000580
0.030000	0.001303
0.040000	0.002314
...             ...

Onde a primeira coluna é a diferença entre o tempo da medição e o tempo de referência. A segunda coluna é o valor de $\langle r^{2}\rangle$ para este intervalo de tempo.

Análise de resultados e exemplos

A seguir alguns exemplos de gráficos do Mean Square Displacement podem ser observados para diferentes densidades do sistema simulado com potencial do tipo Lennard-Jones. Para produzir os dados de todas as figuras, utilizou-se um passo de tamanho 0.01, esperando-se 5000 passos para começar as medições, e então fez-se a média sobre 9 medições por partícula (correspondendo à um tempo final de 500 e um numMSD de 5000).

• 

  MSD para uma densidade de 5%. É possível perceber o comportamento balístico inicial e então o comportamento linear. O ajuste linear foi realizado no intervalo [10, 50].

• 

  MSD para uma densidade de 20%. É possível perceber um breve comportamento balístico inicial e então o comportamento linear. O ajuste linear foi realizado no intervalo [10, 50].

• 

  MSD para uma densidade de 56%. Para altas densidades o comportamento balístico inicial já não é tão notável, mas o comportamento linear continua ocorrendo. O ajuste linear foi realizado no intervalo [10, 50].

• 

  MSD para uma densidade de 81%. Para altas densidades nota-se o comportamento linear altamente ruídoso, devido a vibração das partículas. O tempo para atingir um comportamento linear começa a ser maior. O ajuste linear foi realizado no intervalo [15, 50]. Na escala do eixo ordenado é possível observar que os deslocamentos são pequenos em unidades de comprimento reduzido.

• 

  MSD para uma densidade de 100%. Para altas densidades nota-se o comportamento linear altamente ruídoso, devido a vibração das partículas. O ajuste linear foi realizado no intervalo [10, 50]. Na escala do eixo ordenado é possível observar que os deslocamentos são pequenos em unidades de comprimento reduzido.

Análises do expoente da relação do MSD com o tempo nos diferentes intervalos pode ser obtida através de um gráfico log x log.

Note que, pela fórmula apresentada anteriormente, podemos obter o coeficiente de difusão D através das inclinações das retas ajustadas. Note que quanto maior a densidade, menor este coeficiente, pois menos difuso é o sistema.

Autocorrelação de velocidades

Pressupostos teóricos

A Autocorrelação de velocidades é uma medida da variação das velocidades das partículas em um sistema em função do intervalo do tempo. Ela nos diz o quanto a velocidade está correlacionada com a velocidade a $\tau$ passos atrás. Se definirmos $\tau =t''-t'$ , ou seja, uma grandeza de intervalo de tempo, então a função correlação das velocidades é dada por:

$Z_{x}(\tau )={\frac {1}{N}}\sum _{i=1}^{N}v_{x,i}(\tau )v_{x,i}(0)$

O mesmo vale para a componente cartesiana y da velocidade. $\tau =0$ é um instante arbitrariamente escolhido.

Sua interpretação pode ser dividida em dois momentos:

• Para intervalo de tempo curto: em média, a velocidade das partículas tem a mesma direção que ela tinha alguns passos atrás. Assim, o produto destas velocidades é positivo e as contribuições da soma sobre as correlações sobre todas as partículas se somam.

• Para intervalo de tempo longo: a velocidade das partículas pode ou não estar na mesma direção que ela tinha a vários passos atrás. Assim, o produto hora é positivo e hora é negativo, e então a soma sobre todas as partículas se anulam. A autocorrelação é aproximadamente zero.

Note que, quanto mais tempo as partículas andam para a mesma direção, mais elas se deslocam do ponto de origem (se difundem!). Assim, podemos relacionar a correlação das velocidades das partículas de nosso sistema de tal forma que quanto maior a correlação entre as velocidades, maior a difusão.

$D=\int _{0}^{\infty }Z_{x}(\tau )d\tau$

Perceba que esta equação relaciona uma grandeza macroscópica (a difusão) com a integral sobre o tempo de uma autocorrelação (das velocidades de cada partícula). Tais relações são chamadas relações de Green-Kubo e podem ser aplicadas para outros coeficientes (por exemplo: viscosidade, condutividade térmica, condutividade elétrica. Mais detalhes verifique ref. [1], página 90).

Dada a integral imprópria, seria impossível obter dados infinitamente por razões óbvias. Felizmente a correlação decaí a zero, e as contribuições à integral tornam-se desprezíveis. Assim, é possível observar o sistema em tempo finito e obter, a partir da integral, um bom valor para a constante de difusão, desde que garanta-se que a correlação de velocidade já tenha chegado à zero (vide exemplos).

Num potencial como Lennard-Jones, uma maior densidade de partículas significa mais forças atuando e portanto a variação da velocidade (aceleração) é maior. Assim, as velocidades em diferentes instantes logo se descorrelacionam. Para sistemas menos densos, a correlação de velocidades se mantém por mais tempo.

Interessante notar que a correlação de velocidades em 0 é a energia cinética da componente x em unidades reduzidas: $Z_{x}(0)=v_{x}(0)v_{x}(0)=K_{x}(0)$. Também é interessante o fato de que a função de correlação de velocidades pode assumir valores negativos, e deste caso, as velocidades estão anti-correlacionadas, ou seja, é mais provável que depois de um intervalo de tempo a componente da velocidade tenha invertido seu sentido.

O limite $Z(t\to \infty )=0$, ou seja, nada podemos prever sobre a velocidade da partícula depois de um tempo infinitamente grande baseado na velocidade dela neste instante.

Implementação computacional

Vamos apenas medir a correlação entre a componente x da velocidade. Nosso programa apenas fornecerá os pontos do gráfico de $Z(\tau ):\tau$ , e, para integrá-lo, é possível usar um algoritmo numérico, ou um programa como o XM Grace.

Declararemos um vetor Zvel, que tenha numZ posições(cada posição corresponderá há um intervalo de passos diferente). Também será necessário um vetor para armazenar a componente x da velocidade no instante definido como zero.

#define numZ 1000 //Intervalo máximo de tempo medido
double *Zvel = (double *)malloc(numZ*sizeof(double)); // Vetor com os valores de Z(t);
double *vx_zero = (double *)malloc(NP*sizeof(double)); // Vetor com os valores de vx(0);

Criaremos uma função vel_correl que mede a correlação entre as velocidades. Para isso, a declaremos da seguinte forma, enviando as velocidades vx atuais, as salvas como $v_{x}(0)$, quantos passos estamos desde a definição de 0 e o vetor que armazenamos os valores calculados:

void vel_correl (double *vx, double *vx_zero, int ref, double *Zvel);

No corpo principal do programa, chamaremos a função vel_correl. Tomaremos o cuidado de garantir que calcularemos a correlação de velocidades após atingirmos um ponto de equilíbrio do sistema (assim, ignoraremos os primeiros numZ passos com um if).

Manteremos também um counter, para sabermos quantas vezes passamos por esta função e acumulamos dados no vetor Zvel. Consideramos que o counter é igual para todos os valores de $\tau$. Caso isso não seja verdade, é necessário um vetor que contabiliza quantos vezes calculou-se o $Z(\tau )$ para cada intervalo.

if(passo >= numZ){
 		vel_correl(vx, vx_zero, passo % numZ, Zvel);
		counter++;
	}

Note que o valor de passo%numZ varia de 0, quanto tomaremos nossos valores de referência, até numZ-1 que será o nosso maior intervalo de passos. No caso deste programa calculamos o Z a cada passo. Geralmente programas não apresentam essa necessidade, o que pode ser facilmente ajustado.

Então, teremos o corpo da função vel_correl, que, caso o parâmetro ref é zero, toma as velocidades atuais como referência. Se não, calcula o produto entre a velocidade atual e a velocidade de referência, fazendo a média sobre essa informação das partículas e salvando no vetor Zvel.

void vel_correl (double *vx, double *vx_zero, int ref, double *Zvel){

int i;

if(ref==0){ //Se a ref é 0, então...

  for(i=0; i<NP; i++)
    vx_zero[i]=vx[i]; //... é o momento de tomar as velocidades atuais como referência

}

double Zaux=0; // Essa variável vai acumular a soma das correlações das partículas

for(i=0; i<NP; i++)
	Zaux = Zaux + vx[i]*vx_zero[i];

Zvel[ ref ]+= Zaux/NP; //Aqui acumulamos a correlação média por partícula das velocidades

}

Por fim, depois do laço temporal, dividiremos o valor de counter por numZ, dando assim o valor de vezes em que se acumulou o valor de Z(t) por intervalo de passos. Para imprimirmos, usaremos o intervalo de passo em intervalo de tempo reduzido, multiplicando o número de passos pelo tamanho de cada passo.

counter = counter/numZ;    //Counter marcará quantas vezes se passou pela função para cada 
                          //valor de intervalo de passos ref.

for(i=0; i<numZ; i++){
	fprintf(arq, "%lf\t%lf\n", i*dt, (double)Zvel[i]/counter); //Impressão
}

No arquivo então constará:

0.000000	0.484037
0.010000	0.465207
0.020000	0.463353
0.030000	0.460682
0.040000	0.457317
...             ...

Onde a primeira coluna é a diferença de tempo reduzido entre o valor de referência e o valor de medição, e a segundo coluna é o valor de $Z(\tau )$ médio do sistema.

Para integrar esta função no XM Grace, por exemplo, para obter o valor do coeficiente de difusão, primeiro é necessário garantir que tomou-se numZ grande o suficiente para a correlação ter caído para próximo de 0. Então é possível abrir o programa no terminal:

> xmgrace

E ir no menu Data>Import, selecione o arquivo recém gerado, dê Ok e Cancel. Depois, vá em Data>Transformation>Integration, selecione em Load a opção Sum only e o valor registrado na caixa abaixo da opção é o valor da integral sobre toda a curva importada.

Exemplos e análises dos resultados

A seguir alguns exemplos de gráficos da Correlação de Velocidades em função do intervalo de tempo $\tau$ podem ser observados para diferentes densidades do sistema simulado com um potencial Lennard-Jones. Para produzir os dados de todas as figuras, utilizou-se um passo de tamanho 0.01, esperando-se 5000 passos para começar as medições, e então fez-se a média sobre 9 medições por partícula (correspondendo à um tempo final de 500 e um numZ de 5000).

• 

  Z(t) para uma densidade de 5%. É possível perceber o comportamento correlacionado da velocidade para pequenos intervalos de tempo. A integração foi realizada no intervalo [0, 50].

• 

  Z(t) para uma densidade de 20%. É possível perceber o comportamento correlacionado da velocidade para pequenos intervalos de tempo e que esta correlação decaí rapidamente à zero, onde aparece ruidoso. A integração foi realizada no intervalo [0, 50].

• 

  Z(t) para uma densidade de 56%. É possível perceber um brevíssimo comportamento correlacionado da velocidade e em seguida um rápido decaimento à zero. A integração foi realizada no intervalo [0, 50].

• 

  Z(t) para uma densidade de 81%. Não é possível perceber um comportamento correlacionado da velocidade, senão para um tempo brevíssimo em que as velocidades aparecem anti-correlacionadas, ou seja, as partículas vibram com sentidos opostos. A integração foi realizada no intervalo [0, 50].

• 

  Z(t) para uma densidade de 100%. Para altas densidades nota-se a falta de correlação e o comportamento altamente ruidoso, pois devido a vibração das partículas, as velocidades trocam de orientação rapidamente. A integração foi realizada no intervalo [0, 50].

As integrações foram realizadas com o programa XM Grace.

Note que, pela fórmula apresentada anteriormente, podemos obter o coeficiente de difusão D através do valor da integral. Note que quanto maior a densidade, menor este coeficiente, pois menos difuso é o sistema.

Constante de difusão

Scatterplot para a constante de difusão quando calculadas via os dois métodos sugeridos nesta página.

É possível assim comparar as constantes de difusão obtidas para as diferentes densidades de partículas num sistema com potencial do tipo Lennard-Jones. Lembrando que, neste caso, sua unidade é $[D]={\sqrt {\frac {\epsilon }{m}}}\sigma$.

Para calculá-la via MSD, basta dividir por 4 o coeficiente da reta ajustada. Já pelo método da autocorrelação de velocidades, o próprio valor da integral é a constante de difusão.

É possível verificar que os valores em geral se aproximam, especialmente para densidade maiores (neste caso deve haver o cuidado na escolha da região para o ajuste da reta, dado que o coeficiente tende a se anular e tonar constante), já para poucas partículas, onde a difusão é maior, a diferença se tornou um pouco mais relevante e acredita-se que pode ser reduzida com médias sobre mais dados.

Referências

[1] Frenkel, Daan and Smit, Berend (2001). Understanding Molecular Simulation. Academic Press.