Embora seja teoricamente possível adquirir informações sobre o ensemble de configurações nos estados nucleares resolvendo-se a equação de onda para os núcleos (e um potencial eletrônico paramétrico), esta técnica pode ser aplicada apenas para sistemas muito simples. A solução das equações de onda nucleares rapidamente se torna impraticável à medida que o tamanho do sistema aumenta.
Felizmente, mesmo que seja muito difícil ou até mesmo impossível estudar a dinâmica dos sistemas utilizados neste trabalho aplicando-se o formalismo da mecânica quântica, ainda é possível amostrar as configurações nucleares propagando o sistema no tempo segundo equações de mecânica clássica.
Embora seja apenas uma aproximação da dinâmica nuclear real, simulações de mecânica clássica fornecem descrições bastante exatas em muitos casos, apesar de ignorar alguns fenômenos que são importantes em certas circunstâncias. Por exemplo, o tunelamento no tratamento de partículas muito leves ou a vibração de ponto zero em temperaturas muito baixas.
Por este motivo, a aplicação da mecânica clássica está restrita a sistemas onde o comportamento quântico médio do movimento nuclear possa ser confiavelmente modelado a partir da mecânica clássica. Formalmente, a aproximação clássica descreve bem a dinâmica de um sistema atômico-molecular se o comprimento de onda térmico de De Broglie (11) for muito menor que a distância média de separação entre os núcleos do sistema.
=
2 ℏ2
M kBT (11)
A mecânica clássica pode ser formalizada definindo-se um Hamiltoniano clássico do sistema, que é uma função das coordenadas de posição e de momento linear (12).
H pi, ri=
∑
i=1 N 1 2 mi pi 2 U ri (12)Os vetores pi representam as coordenadas de momento linear para cada partícula i e os vetores ri representam as coordenadas de posição. Ao todo, portanto, o Hamiltoniano clássico é uma função de 6N variáveis, onde N é o número de partículas (neste caso, núcleos atômicos).
Em um sistema dinâmico, dizemos que posições e momentos variam no tempo. Descobrir a forma como esta variação ocorre é o problema central da mecânica clássica. Na mecânica de Hamilton, o comportamento dinâmico de um sistema Hamiltoniano pode ser capturado a partir da diferencial total da função Hamiltoniana, considerando-se que as posições e momentos das partículas são função do tempo. Esse resultado é representado pelas equações de movimento de Hamilton (13,14).
˙pi=− ∂ H ∂ ri (13) ˙ri= ∂ H ∂ pi = pi mi (14)
Resolver a dinâmica do sistema envolve encontrar as velocidades e posições das partículas em um tempo posterior a partir de tempos anteriores, isto consiste em resolver estas equações diferenciais de primeira ordem para cada uma das partículas, a cada passo de tempo, e encontrar as posições e velocidades no intervalo de tempo posterior.
Para resultados de alta exatidão, é preferível resolver as duas equações diferenciais de primeira ordem de Hamilton (13,14). Entretanto, as partir da definição clássica de momento linear, é possível chegar, a partir das equações de Hamilton, a uma equação diferencial de segunda ordem
(15), que nada mais é que a equação de Newton e relaciona as forças sobre a partícula a uma
derivada temporal da posição.
Fi=mia (15)
Resolver a equação de Newton é mais eficiente computacionalmente que resolver as duas equações de Hamilton, apesar de ser menos exato, por este motivo normalmente os métodos de dinâmica molecular integram a equação de Newton.
Os algoritmos utilizados para integrar as equações de movimento da qual estamos tratando agora são historicamente conhecidos como algoritmos tipo Verlet, em homenagem ao físico francês Loup Verlet, que foi o primeiro a aplicar métodos de integração numérica às equações de movimento em simulações de dinâmica molecular (VERLET, 1967).
Considerando a posição como função do tempo (o indispensável para se pensar em um sistema dinâmico) de modo que em um tempo t a posição dos átomos no sistema seja Rt , a ideia fundamental dos algoritmos tipo Verlet está em escrever duas expansões de Taylor da posição: uma para um passo à frente (16) no tempo e outra para um passo atrás (17). Em ambas, a expansão é truncada no terceiro termo.
Rt =Rt ˙Rt 2 2 ¨Rt O 3 (16) Rt− = Rt − ˙Rt 2 2 ¨Rt −O 3 (17)
Adicionando-se estas duas expressões, conseguimos uma terceira que relaciona as coordenadas no futuro às coordenadas no presente e no passado (18).
Rt ≈2 Rt−Rt− 2M−1Ft 18
Subtraindo-se as equações (16) de (17), obtemos uma equação para as velocidades das partículas em um tempo t em função da velocidade que depende das posições um passo atrás e um passo à frente, mas não no passo atual. Perceba que à medida que o passo de tempo torna-se muito pequeno, esta definição aproxima-se da definição de derivada da posição, o que é consistente.
Vt=[ Rt −Rt− ]
2 19
Apesar de suficientes para a integração das equações de movimento, as equações (16) e (17) não são muito apropriadas. O motivo disto é que utilizando estas equações para propagar o sistema, o cálculo das velocidades está sempre defasado em relação ao das posições. Isto é, para saber
Vt é preciso antes saber Rt .
Esta característica é levemente inconveniente para as rotinas de saída de dados mas também tem uma complicação adicional de ordem prática: no estado inicial do sistema, as velocidades não podem ser calculadas porque não se tem qualquer informação sobre o estado futuro das posições e portanto é preciso usar uma forma diferente de se calcular os valores no estado inicial.
Uma variação do algoritmo de Verlet convencional, chamado criativamente de velocity-Verlet, resolve este problema integrando as velocidades no passo atual utilizando as forças do passo seguinte. Este algoritmo parte das definições de velocidade em (21') e (21).
Rt =Rt ˙Rt 2 2 ¨Rt O 3 Rt ≈R t V t 2 2 M −1 Ft 20 1 [ Rt −Rt− 2 2 M −1 Ft]≈V t Vt ≈1 [Rt −Rt − 2 2 M −1 Ft ] Vt ≈Rt − Rt − 2 M −1 Ft Vt≈[ Rt −Rt] − 2 M −1 Ft 21'
Note que a quantidade [ Rt − Rt]
é uma velocidade segundo a definição (19), essa
escrever que Vt1/ 2 =[ Rt −R t]
. Substituindo-se esta definição em (21') segue-se
que Vt ≈V t1 /2 − 2 M
−1
Ft . Que ao ser rearranjada assume a forma da equação (21):
−V t1/2 ≈−V t − 2 M −1 Ft Vt1/2 ≈V t 2 M −1 Ft Vt ≈V t1/2 2 M −1 Ft Vt t−V t≈ 2 M −1 [ Ft F t] Vt t≈V t 2 M −1[ Ft F t] 21
A equação (21) pode ser utilizada para calcular as velocidades em tempos posteriores a partir do conhecimento das velocidades no estado atual do sistema, dependendo no entanto da informação sobre as forças.
Em sistemas conservativos (onde o potencial depende apenas da posição das partículas), as forças podem ser calculadas a partir do gradiente do potencial.
Ft=−∇ U [ Rt] (22)
Como as posições posteriores são calculadas primeiro, com a equação (20), é possível calcular as forças a partir da equação (22) e então calcular as velocidades no mesmo passo temporal das posições a partir da equação (21), solucionando o problema do algoritmo básico de Verlet.
Note, no entanto, que no passo inicial t = 0, as posições são aquelas da geometria de entrada, mas as velocidades iniciais precisam ser atribuídas de forma externa ao algoritmo de integração, pois elas são necessárias para se calcular as próximas posições e portanto as próximas forças, para então encontrar as próximas velocidades.
Também é necessário saber o tamanho do passo de integração, uma vez que todas as quantidades dependem desta informação para serem apropriadamente calculadas a cada iteração do procedimento de integração.
O ciclo de integração então segue-se da seguinte forma, ilustrado na Figura 9:
2. Atribua valores iniciais de posições R0 e velocidades V0 ; 3. Defina um passo temporal de integração;
4. Determine o número nt de passos a serem integrados; 5. Calcule as forças F0 para as posições iniciais;
6. Itere ao longo de nt passos:
1. Calcule as novas posições no passo t com (20); 2. Calcule as novas forças no passo t com (22); 3. Calcule as novas velocidades em t com (21);
4. Escreva informações pertinentes sobre a integração do passo t para t ; 5. Incremente o passo de tempo;
7. Analise os resultados.
Figura 9: Ciclo de Integração do Velocity-Verlet
As posições iniciais R0 são obtidas a partir da geometria de entrada. A composição do sistema (número e tipos de átomos) também precisam ser informadas no início do procedimento. Estes parâmetros são utilizados no cálculo do potencial além de outros termos.
Uma forma conveniente de se estabelecer velocidades iniciais V0 é gerando-as aleatoriamente a partir da distribuição de Maxwell-Boltzmann (23). Este método permite que as velocidades sejam atribuídas de forma que o sistema tenha uma determinada “temperatura instantânea” no começo da simulação.
f v dv=
mi2 KBT exp
− miv2
2 KBT
dv (23)Esta forma de se atribuir velocidades é justificada por um resultado da termodinâmica estatística que relaciona a média da energia cinética K do sistema com sua temperatura macroscópica (24).
T= 1
NAKBT 〈K 〉 (24)
Analogamente, podemos pensar que uma temperatura instantânea T que está relacionada à energia cinética instantânea da trajetória (25).
T= 1
NAKBT K (25)
Esta temperatura instantânea pode ser calculada a partir das velocidades instantâneas de cada microestado por onde a dinâmica passa. E portanto, as velocidades iniciais têm influência na temperatura do microestado inicial.
O passo de integração é uma variável de entrada cuja atribuição não é muito trivial. Quanto mais curto o passo, menos tempo “real” a dinâmica será capaz de simular antes de teoricamente faltar espaço para o armazenamento de dados de saída. Entretanto, um passo muito grande compromete o algoritmo de integração porque a expansão de Taylor é uma aproximação local.
Além disso, a “resolução” da dinâmica é determinada pelo tamanho desse passo. Se o passo é muito largo, fenômenos que acontecem em intervalos de tempo inferiores ao tamanho do passo não são adequadamente amostrados do ensemble.
Por esse motivo, o passo de integração acaba sendo restrito a intervalos bastante curtos para não prejudicar a exatidão do algoritmo de integração e para amostrar apropriadamente os graus de liberdade mais rápidos do sistema.
Em geral, estes graus de liberdades mais rápidos correspondem a estiramentos de ligações envolvendo hidrogênios, particularmente a ligação carbono-hidrogênio, cuja vibração é cerca de 3000 cm-1 que ocorrem na escala de tempo da ordem de femtossegundos (10-15 s). Na prática, portanto, o intervalo de integração utilizado em simulações de dinâmica molecular usuais é de 1fs.
O potencial é mais complicado. Dissemos que, segundo a aproximação Born-Oppenheimer, é possível simular o movimento nuclear a partir do conhecimento da superfície potencial derivada da resolução paramétrica da função de onda eletrônica para cada configuração especial nuclear.
Em teoria, seria possível utilizar este potencial quântico, por assim dizer, para movimentar os núcleos. Entretanto, o custo computacional de se calcular uma função de onda eletrônica para cada passo de 1fs cresce muito rapidamente com a quantidade e os tipos de átomos no sistema.
Esta lentidão na propagação dos sistema no tempo faz com que cálculos de dinâmica molecular utilizando um potencial deste tipo sejam realizados apenas para sistemas muito pequenos, onde o tempo de convergência do cálculo do potencial não seja um fator tão limitante.