6. Noções básicas de Dinâmica Molecular Clássica

Antes de inicarmos o estudo com a Dinâmica Molecular (DM) veremos o conceito de Mecânica Molecular (MM)^6.1.

Mecânica[15] Definição: Estudo das interações entre a matéria e as forças que agem nela. Estática, está amplamente relacionada com a ação das forças quando não há nenhuma mudança do momento, enquanto que a dinâmica trata dos casos em que ocorrem um mudança do momento. Cinemática é o estudo do movimento dos corpos sem referência às forças que afetam o movimento. Ciências clássicas estão relacionadas com os corpos macroscópicos no estado sólido, enquanto que a mecânica dos fluidos é a ciência das interações entre forças e líquidos.

Baseado na definição de mecânica, no ponto de vista da física, atribuimos a MM como um tratamento estático de otimização das posições de um sistema de átomos ou moléculas, em que as forças são geradas pelos potenciais atômicos^6.2. A MM [16] tem um papel muito importante na busca da geometria molecular de sistemas com muitos átomos, por sua simplicidade comparada aos métodos quânticos. Como característica do método clássico de MM não temos a informação da parte eletrônica como no método quântico. Uma referência no uso da MM são os sistemas biológicos de proteínas [17], atualmente um dos limites na simulação atomística de sistemas orgânicos, podendo envolver centenas de milhares (ou milhões) de átomos.

A descrição mais simples do método de MM é considerar a aproximação de Bohr & Oppenheimer. Esta aproximação leva em conta que o movimento dos núcleos é mais lento que o movimento dos elétrons, então podemos separar a informação nuclear e eletrônica em duas partes e resolvê-las separadamente. Dessa aproximação (da mecânica quântica), observamos que em MM a energia total do sistema depende exclusivamente da posição dos átomos do sistema e os efeitos eletrônicos não são computados explicitamente. A energia total então é dada via potencial nuclear dependente das posições, mais conhecido pela denominação de campo de força (CF).

O campo de força é uma peça fundamental, se não a mais importante, em MM. Para se ter uma bom resultado de geometria em MM é necessário que o CF esteja adequado ao tipo de sistema no qual deseja tratar. Comunmente os CF são compostos por termos harmômicos para átomos ligados, termos de van der Walls e de Coulomb para átomos não ligados. Os termos dos átomos ligados têm a forma de $kx^{2}$, onde $x$ pode assumir valores de distância ou ângulo. Com o CF, a força do sistema é obtida via força central pela derivada da expressão espacial do potencial, que é a derivada da expressão do CF. Por fim, o sistema deve ser conservativo, uma vez que o potencial central impõe a conservação do momento angular total do sistema.

Para exemplificar um CF, na equação 6.0.1 temos um CF que é utilizado para MM também utilizado para DM, este campo de força na literatura é encontrado com o nome de Universal Force Field.

$$E_{\text{potencial}} = \begin{cases} \sum_{r} \frac{1}{2}K_{IJ}\l... ...van der Waals \par + \sum_{R} Q_{I}Q_{J}/ \epsilon R_{IJ} % coulomb \end{cases}$$ (6.0.1)

A equação 6.0.1 se refere ao CF Universal Force Field (UFF) [18]. A descrição dos termos da equação 6.0.1 é dado por:

• $\sum_{r}$, $\sum_{\theta}$, $\sum_{\phi}$, $\sum_{\omega}$, $\sum_{x}$ e $\sum_{R}$ são os termos de ligação, angular, torção, inversão, van der Waals (vdw) e coulombiano, respectivamente ;

• $K_{IJ}$, $K_{IJK}$ e $K_{IJKL}$ são constante de força de ligação, angular e torção/inversão, respectivamente;

• $r$ e $r_{IJ}$ são as distâncias calculadas e de equilíbrio entre dois átomos ligados $I$-$J$;

• $\theta$ e $\theta_{0}$ são os ângulos calculados e de equilíbrio entre três átomos $I$-$J$-$K$;

• $\phi$ e $\phi_{0}$ são os ângulos de torção calculados e de equilíbrio entre quatro átomos $I$-$J$-$K$-$L$;

• $\omega$ é o ângulo de inversão do átomo $I$ em relação ao plano formado pelos átomos $J$-$K$-$L$;

• $C_{0}^{\theta}$, $C_{1}^{\theta}$ e $C_{2}^{\theta}$ são os coeficientes de Fourier para o caso angular geral;

• $C_{0}^{\omega}$, $C_{1}^{\omega}$ e $C_{2}^{\omega}$ - são os coeficientes de Fourier do termo de inversão;

• $V_{\phi}$ e $n$ - são a altura da barreira de torção e a periodicidade do potencial de torsão;

• $x$ e $x_{IJ}$ - são as distâncias calculadas e de equilíbrio entre dois átomos não ligados $I$ e $J$;

• $D_{IJ}$ - é a profundidade do potencial de Lennard-Jones;

• $Q_{I}$ e $\epsilon$ - são as cargas parciais do átomo $I$ e a constante dielétrica.

Em um sistema molecular dizemos que o sistema está em equilíbrio quando o somatório de todas as forças que atuam em cada átomo tendem para zero. Para completar o processo de MM é necessário um método de comparação entre os passos de otimização. Cada passo de otimização está relacionado com um pequeno movimento do átomo, em geral na direção do ponto de equilíbrio, devido a ação da força proveniente do potencial. Em um sistema multi-atômico modificar as posições dos átomos e calcular a energia total até encontrar um mínimo entre as possíveis configurações pode ser um tanto quanto trabalhoso. Para isso utilizamos outros recursos que nos auxiliam na busca da geometria como o método do gradiente conjugado [16]. Além do método do gradiente conjugado, outros métodos matemáticos de minização utilizam a expressão de energia potencial (energia total) como critério de convergência, onde, um sistema com uma boa geometria, apresenta em geral, diferença de energia entre um ciclo e outro menor que $10^{-1} kcal/mol$ Å [16] após os passos de otimização.

Do processo da MM utilizamos o conceito de CF para resolver as equações de movimento no tempo, e a resolução dessas equações é então chamada de DM. A DM [19,20,21] apresenta conceitualmente algumas características somadas ao processo dinâmico do sistema molecular e inclui o termo de energia cinética na energia total, o que não acontece em MM que pois esta representa a otimização de geometria de um sistema molecular à temperatura de zero Kelvin. As principais características da DM são: resolução das equações de movimento (segunda lei de Newton) com dependência temporal e ajuste da temperatura do sistema (no caso do ensamble seja canônico (NVT) a temperatura é constante).

Como já sabemos se o sistema for conservativo podemos obter a força do sistema a partir da energia potencial do sistema como descrito pela equação 6.0.2.

$$\vec{F}=-\vec{\nabla}U$$ (6.0.2)

$$\vec{F}=m\vec{a}=m\frac{d\vec{v}}{dt}=m\frac{d^{2}\vec{r}}{dt^{2}}$$ (6.0.3)

$$\vec{a}=-\frac{\vec{\nabla}U}{m}$$ (6.0.4)

$$\frac{d\vec{v}}{dt}=\vec{a}$$ (6.0.5)

$$\frac{d\vec{r}}{dt}=\vec{v}$$ (6.0.6)

Pela segunda lei de Newton (equação 6.0.3) podemos encontrar a aceleração (equação 6.0.4) utilizando e derivada da expressão do potencial, ou como chamamos acima de CF, a equação 6.0.2. Para se obter as velocidades e posições é necessário integrar as equações 6.0.5 e 6.0.6. A integração das equações de movimento no processo de dinâmica molecular é realizada numericamente ou por expressões analíticas. Dentre as técnicas de integração uma das mais utilizadas é o algorítmo de Verlet, para o desenvolvimento das nossas atividades utilizaremos o algortimo Velocity Verlet (eq. 5.3.14).

Com esse breve histórico de MM e DM, podemos iniciar a construção do nosso programa de dinâmica molecular. No programa a ser construído levaremos em conta os seguintes aspectos:

1. O sistema que iremos simular é um gás de Argônio (gás inerte - Ar) dentro de uma caixa. Assim teremos um número fixo de átomos de Ar dentro de uma caixa de dimensões a ser definida. A inclusão da caixa indica que teremos que utilizar condições periódicas de contorno;
2. Utilizaremos todas as dimensões baseado no sistema internacional de unidades (SI);
3. Utilizaremos somente interação de van der Walls (vdW) como o potencial entre os átomos, a implementação de outros tipos de potenciais demanda muito tempo por isso não faremos aqui. Este tipo de potencial se encaixa na interação entre átomos não ligados;
4. Uma vez definida uma caixa de dimensões $x \times y \times z$ atribuiremos as posições aleatoriamente dos átomos de Ar (x, y, z) dentro das dimensões da caixa.
5. Com as posições definidas, atribuiremos as velocidades aleatoriamente baseado em uma distribuição gaussiana;
6. Com as velocidades teremos uma energia cinética, assim as velocidades deverão ser normalizadas para uma temperatura em Kelvin baseada no Princípío da Equipartição de Energia, onde o pico da distribuição gaussiana das velocidades deve estar relacionado com a temperatura desejada;
7. O banho térmico será simples dado pelo reescalamento das velocidades;
8. Durante a resolução das equações no tempo iremos guardar (escreve em um arquivo) a energia cinética, energia potencial, temperatura, passo, (em um outro arquivo) as posições atômicas.
9. É desejável que o programa da geração das coordendas seja independente do programa de DM.
10. Também é desejável que o programa de DM utilize dois arquivos de entrada um com as coordenadas e outro com os dados de controle da DM, como por exemplo, temperatura, $\Delta t$, números de passos, tamanho da caixa, números de passos para escrever os dados, número de passos para escrever a trajetória, etc.