Medidas dinâmicas

De Física Computacional
Ir para navegação Ir para pesquisar

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 correlação entre velocidades.

Ambas essas medidas estão relacionadas com a difusão, o modo em que uma substâncias 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.

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. [1]


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 do 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:

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 consegiu 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):

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. 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:
  • Regime subdifusivo:
  • Regime normalmente difusivo:
  • Regime superdifusivo:
  • Regime balístico:

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.

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 anteior 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.

Difiniremos 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 referencia 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 . Caso isso não seja verdade, é necessário um vetor que contabiliza quantos vezes calculou-se o MSD para cada intervalo.


if(passo >= numMSD){
 		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 . 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 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.


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.

Correlação de velocidades

Pressupostos teóricos

A Correlação de velocidades é uma medida da variação das velocidades das partículas em um sistema em função do intervalo do tempo. Se definirmos , ou seja, uma grandeza de intervalo de tempo, então a função correlação das velocidades é dada por:

O mesmo vale para a componente cartesiana y da velocidade. é 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 horas é positivo e horas negativo, e então a soma sobre todas as partículas se anulam. A correlaçã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, podes 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.

Note que esta equação relaciona uma grandeza macroscópica (a difusão) com uma grandeza microscópica (a velocidade de cada partícula). Tais relações são chamadas relações de Green-Kubo.

Dada a integral imprópria, seria impossível obter dados infinitamente. 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.

Num potencial como Lennard-Jones, uma maior densidade de partículas significa mais forças atuando e portanto a variação da velocidade é 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.

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 , 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 , 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 . Caso isso não seja verdade, é necessário um vetor que contabiliza quantos vezes calculou-se o 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 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 podem ser observados para diferentes densidades do sistema simulado com um potencial Lennard-Jones.

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.

Referências

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