Grupo - Lennard Jones: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com 'A evolução temporal do estado quântico <math> \Psi(\mathbf{r},t) </math> é dada pela equação de Schrödinger, a qual é postulada como: <math>i\hbar\frac{\partial}{\par...')
 
Sem resumo de edição
Linha 1: Linha 1:
A evolução temporal do estado quântico <math> \Psi(\mathbf{r},t) </math> é dada pela equação de Schrödinger, a qual é postulada como:
O potencial devido a interação entre duas partículas separadas por uma distância <math>r'</math> pode ser modelado pelo potencial de Lennar Jones:


<math>i\hbar\frac{\partial}{\partial t} \Psi(\mathbf{r},t) = \left [ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r},t)\right ] \Psi(\mathbf{r},t)</math>
<math>U(r') = 4 \epsilon \left [ \left ( \frac{\sigma}{r'} \right )^{12} + \left ( \frac{\sigma}{r'} \right )^{6} \right ].</math>


Posto em unidades atômicas (onde <math>m_e</math> e <math>\hbar</math> são unitários), o caso unidimensional de um elétron num potencial independente do tempo reduz-se a:   
Posto em unidades reduzidas (<math>r \equiv r'/ \sigma</math> e <math>U \equiv U' / \epsilon</math>), o potencial reduz-se a:   


<math>\frac{\partial}{\partial t} \Psi(x,t) = \left [ \frac{i}{2}\frac{\partial^2}{\partial x^2} - i V(x)\right ] \Psi(x,t)</math>
<math>U(r) = 4  \left [ \left ( \frac{1}{r} \right )^{12} + \left ( \frac{1}{r} \right )^{6} \right ].</math>


No caso tem-se, por conveniência, o seguintes sistema de unidades básicas:


== Método numérico ==


Buscando resolver a equação numericamente, tem-se a discretização de <math>\frac{\partial ^2 \Psi}{\partial x^2}</math> :
{| class="wikitable"
|-
! '''Grandeza'''
! '''Comprimero'''
! '''Tempo'''
! '''Massa'''
! '''Temperatura'''
! '''Energia'''
! '''Pressão'''
! '''Densidade'''
|-
| '''Unidade'''
| <math>\sigma</math>
| <math>\sigma \sqrt{m_p / \epsilon}</math>
| <math>m_p</math>
| <math>\epsilon/k_B</math>
| <math>\epsilon</math>
| <math>\epsilon / \sigma^3</math>
| <math>1 / \sigma^{3}</math>
|}
onde <math>m_p</math> é a massa da partícula e <math>k_B</math> é a constante de Boltzmann.
.
== Método Monte Carlo ==
 
Dado uma amostra com <math>N</math> partículas, a abordagem introduzida por Metropolis segue o seguinte esquema:


<math>\frac{\Psi^{n}_{j-1} - 2\Psi^{n}_{j} + \Psi^{n}_{j+1}}{\left(\Delta x \right)^2}</math>
<math>\frac{\Psi^{n}_{j-1} - 2\Psi^{n}_{j} + \Psi^{n}_{j+1}}{\left(\Delta x \right)^2}</math>

Edição das 14h04min de 10 de janeiro de 2018

O potencial devido a interação entre duas partículas separadas por uma distância pode ser modelado pelo potencial de Lennar Jones:

Posto em unidades reduzidas ( e ), o potencial reduz-se a:

No caso tem-se, por conveniência, o seguintes sistema de unidades básicas:


Grandeza Comprimero Tempo Massa Temperatura Energia Pressão Densidade
Unidade

onde é a massa da partícula e é a constante de Boltzmann. .

Método Monte Carlo

Dado uma amostra com partículas, a abordagem introduzida por Metropolis segue o seguinte esquema:

e as discretizações de (explícita e implícita, respectivamente):

(explícita)

(implícita)

Tanto no método explícito quanto no método implícito não é conservada a norma do estado (o que é estritamente necessário, já que o estado pode ser interpretado como uma onda de probabilidade). Por esse motivo, utiliza-se o método de Crank-Nicolson, o qual tem essa propriedade.

O método de Crank-Nicolson consiste em uma média aritmética dos métodos explícito e implícito. Excetuando manipulações algébricas triviais, verifica-se que a relação de recorrência do método é dada por:

onde

e .

A integração numérica depende, portanto, do potencial em que o elétron está sujeito, bem como da sua condição inicial e suas das condições de contorno.

Que condições podemos impor para a fronteira? Quando se trata do problema analiticamente, costuma-se considerar que a função de onda tende a zero no infinito. Numericamente, pode-se fazer uma transposição disso, criando uma condição para bordas em pontos suficientemente distantes do centro da distribuição da função de onda, igualando-as a zero. Outra forma de tratar o problema numericamente é criando condições de contorno periódicas, em que para as bordas vale para todo (ou, para as bordas e há a relação para todo ).

Condições de contorno iguais a zero

Para as condições de contorno , a iteração reduz-se à equação matricial:


Condições de contorno periódicas

De maneira semelhante, a iteração do caso das condições de contorno periódicas - - reduz-se à equação matricial:


Condição inicial

Já a condição inicial é arbitrária, pois define o estado inicial do sistema que queremos tratar. Fazendo uma referência ao tratamento de sistemas clássicos, seria como definir posição e momento iniciais. É claro que, para ter o sentido físico de uma função de onda, deve-se ter o cuidado de criar uma condição inicial normalizada, satisfazendo

bastando, então, inseri-la no programa.

Implementação em C


#include <stdio.h>
#include <math.h>
#include <complex.h> 

#define N_steps 100000
#define L 1000
#define dt 1
#define dx 4.0
#define w 0.002

double V(int x_rel)
{
	return - pow(w,2) * pow(x_rel - 500,2) / 2.0; //SHO

	//return 0; 
}

double complex b(int i)
{
	return 0.5 * dt * I * (1.0 / pow(dx,2.0) + V(i)) + 1.0;
}

double u_0(int x_rel)
{
	//return sqrt(2 / L) * sin(5*x_rel * M_PI / L); //eigenstate of infinite square well

	//return (1.0/ sqrt(5 * M_PI)) * exp(I * 0.5 * x_rel) * exp(-pow(x_rel - 500, 2) / (2 * 25)); //gaussian packet

	return  (w / M_PI) * exp(-w * pow(x_rel - 500, 2) / 2.0);
}

void CN(double complex *u, double complex *u_aux, double complex *u_next, double complex a, int bi)
{
	//bi = 1 (bounded); bi = 0 (periodic)

	int i;

	for(i = bi; i <= L - bi; i++) u_aux[i] = conj(a) * (u[(i-1+L)%L] + u[(i+1+L)%L]) + conj(b(i)) * u[i];

	//thomas algorithm

	double complex c_new[L+1], d_new[L+1];

	c_new[bi] = a / b(bi);
	for(i = 1 + bi; i <= L - bi; i++) c_new[i] = a / (b(i) - c_new[i-1] * a);

	d_new[bi] = u_aux[bi] / b(bi);
	for(i = 1 + bi; i <= L - bi; i++) d_new[i] = (u_aux[i] - d_new[i-1] * a) / (b(i) - c_new[i-1] * a);

	if (bi == 1) u_next[0] = u_next[L] = 0;
	u_next[L-bi] = d_new[L-bi];
	for(i = L-1-bi; i >= bi; i --) u_next[i] = d_new[i] - c_new[i] * u_next[i+1];
	
	//u = u_next

	for(i = 0; i <= L; i++) u[i] = u_next[i];
} 

int main(void)
{
	int i, j, n = 0;

	double complex u[L+1], u_next[L+1], u_aux[L+1], a = - 0.25 * I * dt / pow(dx,2.0);

	//Initial Contition
	for(i = 0; i <= L; i ++) u[i] = u_0(i);

	while(n < N_steps)
	{
		CN(u, u_aux, u_next, a, 0);
		printf("set title 'Time = %d'\nplot \'-' w lp pt 7 ps 0.1", n);
		for(i = 0; i <= L; i ++) printf("%d\t%.10lf\n",i,pow(creal(u[i]),2) + pow(cimag(u[i]),2));
		//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,cimag(u[i]));
		//for(i = 0; i <= L; i ++) printf("%d\t%lf\n",i,creal(u[i]));
		printf("e\npause 0.1\n"); 
		n ++;
	}

	return 0;
	
}



Referências