Método de Temperaturas Paralelas (Troca de Réplicas) em Dinâmica Molecular e Monte Carlo

Prévia do material em texto

Introdução ao Método de Replica Exchange: Uma Abordagem Usando os
Métodos de Monte Carlo e Dinâmica Molecular
Diego Guedes-Sobrinho∗
Doutorando em Físico-Química no Instituto de Química de São Carlos, Universidade de
São Paulo.
d_guedes@iqsc.usp.br; diego_jem555@hotmail.com
A proposta desse artigo consiste em apresentar alguns desenvolvimentos elementares do método
Parallel Tempering aplicado aos métodos de Monte Carlo e Dinâmica Molecular. A partir de um conciso
apanhado histórico, as limitações do mapeamento do espaço de fase de sistemas complexos − somado à
impossibilidade de uma amostragem completa dos estados − podem ser tomadas como motivação inicial,
além da contextualização de alguns conceitos elementares para estudantes de química e física interessados
no tema.
1 Uma Breve História de Algumas Coisas
O desenvolvimento da termodinâmica a partir do século XVII se constitui um dos
marcos históricos na ciência moderna, uma vez que logrou a junção de duas categorias de
difícil distinção substancial, a saber, os corpos estão em movimento no espaço (dinâmica),
ao passo que possuem propriedades térmicas. O formalismo termodinâmico envolvendo a
dinâmica microscópica da matéria, como legado do determinismo matemático de Galileo
Galilei e Isaac Newton, seria adequadamente proposto somente no final do século XIX
com os trabalhos de Ludwig Boltzmann (1844 − 1906) no ensejo da popularização e do
reconhecimento da física teórica como campo de investigação útil e palatável ao método
científico. Todavia, de Newton a Boltzmann, muitas tensões ocorreram não apenas no
desenvolvimento de ferramentas racionais para investigação da natureza, mas também, na
modelagem que essas transformações proporcionaram à concepção humana no que tange
à atividade científica.
Em meio à revolução industrial, personalidades como Conde de Rumford (1753−1814),
Sadi Carnot (1796−1832), James Prescott Joule (1808−1889) e Rudolf Clausius (1822−
1888), alicerçaram as fundações da termodinâmica com estudos voltados para o desenvol-
vimento de processos de transferência de calor em macroescala − tanto no sentido das
dimensões da matéria quanto no uso desses processos na otimização das máquinas a va-
por. A assertiva de Lawrence Joseph Henderson (1878− 1942)† resume bem o sentimento
científico da época, em que Science owes much more to the steam engine than the steam
engine owes to science (A ciência deve muito mais ao motor a vapor do que o motor a
vapor deve à ciência). Nesse ponto, cabe realçar que a verdade da proposição de Hen-
derson é unilateral, uma vez que um segundo momento na história do desenvolvimento
da termodinâmica, sustentado pelo fortalecimento da concepção atômica, estabeleceria os
sustentáculos para revolução quântica do século XX.
O momento ulterior à termodinâmica clássica foi marcado pela proposta de uma teoria
que estabelecesse uma relação direta entre dois mundos ininteligíveis entre si: de um lado,
a observação experimental; do outro, a previsão a partir do nível molecular. Os trabalhos
∗Para sugestões e melhoramentos: d_guedes@iqsc.usp.br; diego_jem555@hotmail.com
†Da famosa equação de Henderson-Hasselbalch para ajuste entre pH e concentração de um dado ácido
e sua base conjugada na síntese de uma solução tampão.
1
de Boltzmann e, mais tarde, de Josiah Willard Gibbs (1839−1903) utilizaram a base for-
malística da mecânica clássica sob um viés estatístico, resolvendo parcialmente o problema
da complexidade devido à grande quantidade de graus de liberdade de um sistema real.
Desde então, a partir da estatística de Boltzmann e do conceito de ensemble‡ de Gibbs,
fundou-se uma nova disciplina responsável por uma interpretação molecular de sistemas
macroscópicos em equilíbrio à luz dos microssistemas: a termodinâmica estatística. [1–4]
No contexto do que ficou conhecido população de partículas, propriedades como energia
e entropia foram tratadas de modo matematicamente refinado nessa nova abordagem. A
primeira delas, sob um aspecto determinístico; a segunda, admitida a partir do conceito
da aleatoriedade, uma vez que tem por fundamentação etimológica a palavra Entropie, do
alemão, como "ato de virar-se para". A concepção da entropia, antes de ser inteligível sob o
viés configuracional ligado à quantidade de microestados possíveis num sistema,§ pode ser
admitida como rebento direto da evolução do conceito de calor, se aproximando das proto-
definições de temperatura já em Hipócrates (460−370 a.C.) e Claudius Galeno (133−200
d.C.), os quais admitiam uma relação direta da influência do clima na mistura dos fluidos
corpóreos, determinando, assim, questões subjetivas relacionadas ao julgamento humano
e seu intelecto. Desse modo, o próprio termo latino aludia à temperatura − temperare
− a concepção de "mistura", trazendo a forma mnemônica dos gregos e romanos no
entendimento do temperamento humano em suas manifestações melancólicas, coléricas,
fleumáticas e saguíneas. [6]
O entendimento da energia e da entropia como características de todo sistema físico
macroscópico resultou em dois postulados basilares, conhecidos como a primeira e a se-
gunda lei da termodinâmica, nas quais envolve não apenas as relações entre trabalho e
calor − entendidas como formas de transferência de energia e diretamente relacionadas à
ideia de estrutura configuracional (entropia) −, como também, alocam a termodinâmica
como ciência de caráter fenomenológico. Como resultado da contínua e ininterrupta ob-
servação dos dois postulados na natureza, tateamos a concepção de "verdade"no sentido
mais restrito do termo. Em seu Molecular Thermodynamics, Richard E. Dickerson de-
clara: This is the way the world is. (...) Why these two laws are true ir irrelevant (o
mundo é assim, o por quê dessas leis serem verdade é irrelevante). [7] Todavia, há quem
entenda a problemática em torno da primeira e segunda leis sob o aspecto das meras
projeções da mente humana aplicadas à realidade, ou seja, não haveria fenomenologia na
natureza sem a concepção humana. Nas palavras do professor Ingo Müller, da Technische
Universität Berlin, em uma clara expressão da influência do nominalismo, seriam apenas
palavras soltas, pensamento subjetivo e cultura científica. Para sustentar sua assertiva,
Müller utiliza das formulações de Boltzmann e Gibbs como expressões dessa cultura cien-
tífica, possibilitando, assim, a obtenção de resultados tangíveis. [6] Todavia, apesar do uso
da abstração matemática como conhecimento seguro na descrição dos eventos naturais,
o problema é demasiadamente profundo, sendo mais sensato tomar a posição de Dicker-
son em que algumas coisas devem se restringir à especulação filosófica, enquanto outras,
competem ao aperfeiçoamento teórico.
Entre as questões que competem ao método científico, está o uso da estatística boltz-
manniana aplicada a uma determinada população de partículas. No cerne desse for-
‡Coleção mental de um grande número de sistemas, cada um deles construído para ser uma réplica de
um sistema em um nível termodinâmico macroscópico.
§Ou ainda, na interpretação de Ilya Progogine como flecha temporal dos eventos na natureza, em
oposição à proposição de tempo simétrico da mecânica clássica. [5]
2
malismo, está o conceito de distribuição de probabilidades, no qual se define um dado
sistema físico possuindo ni partículas com energia �i dispostas de uma determinada forma
(estado), podendo ser contadas em um ou mais estados de energia definidos pelo nú-
mero de degenerecência gi. De modo conciso, a conexão entre o estado macroscópico de
um dado sistema, também chamado de macroestado − sendo definido pelos parâmetros
relativos às condições do sistema (temperatura, volume, pressão e concentração) −, e o
estado microscópico quanto à disposição molecular, foi uma das contribuições basilares da
termodinâmica estatística. O resultado dessa abordagem culminou na proposta de Gibbs
baseada na coleção de sistemas noensemble, no qual cada microestado é representativo do
sistema macroscópico em questão. Portanto, a partir de determinadas operações (previa-
mente definidas no formalismo da termodinâmica estatística) aplicadas sobre uma função
que descreve o ensemble, chamada função de partição Z, vários parâmetros termodinâ-
micos podem ser conhecidos como médias sobre esse ensemble. Nesse ponto, se postula o
que é conhecido como hipótese ergódica, onde a média de um determinado parâmetro no
tempo corresponde à média do ensemble, uma vez que estamos trabalhando com sistemas
no equilíbrio termodinâmico. A hipótese ergódica corresponde a um exemplo do caráter
fenomenológico da termodinâmica, onde a validade do método de Gibbs se afirma diante
da sistemática confirmação empírica. [?, 2]
A distribuição de probabilidades de Boltzmann, necessária para as operações sobre o
ensemble, requer o conhecimento da amostragem completa do sistema em cada microes-
tado possível. No entanto, o conhecimento analítico da função de partição ficou restrito
para sistemas muito simples, como, por exemplo, os chamados sistemas diluídos de par-
tículas distinguíveis (boltzons) definidos pela estatística de Maxwell-Boltzmann. Nesses
sistemas, o grau de degenerecência gi para cada estado é muito grande com relação ao
número de partículas ni (distinguíveis), de modo que, gi >> ni, possibilitando a deter-
minação exata do ensemble mais simples conhecido como microcanônico (NV E)¶. Desse
modo, para sistemas mais complexos com um número impraticável de estados quânticos
acessíveis e compatíveis com o estado mascroscópico em questão, além da grande pos-
sibilidade de disposição configuracional das partículas no espaço ocupado pelo sistema,
somente no século XX a publicação de uma série de trabalhos colocou a termodinâmica
estatística no rol das ciências teóricas beneficiárias às aspirações experimentalistas.
Essa nova etapa foi marcada pela publicação dos trabalhos de Nicholas Metrópolis e co-
laboradores (entre eles o físico Edward Teller),‖ com o método de Monte Carlo (1953), [9]
e de Alder e Wainwringht com o método de Dinâmica Molecular (1957 e 1959) [10]. Esses
métodos possibilitaram o estudo de sistemas mais complexos, onde a simetria espacial
pode ser ignorada ao passo que aspectos estruturais de cada molécula e as interações de
diferentes naturezas (coulombianas ou de curto alcance) podem ser levadas em conside-
ração. Desde então, se tornou possível (dentro de algumas limitações) o mapeamento do
espaço de fase do sistema por meio de dois tipos de cinéticas: uma estocástica, para o
método Monte Carlo, e outra temporal, para o método de Dinâmica Molecular. Portanto,
a partir da amostragem restrita aos pontos de estabilidade do espaço de fase do sistema
por meio dos dois métodos, as médias termodinâmicas podem ser calculadas de acordo
com a amostragem realizada no ensemble. [11–15]
¶Nesse ensemble, o número de partículas, N , o volume do sistema, V , e a energia, E, são constantes.
‖Teller foi contado entre os cientistas convocados para o Projeto Manhattan. J.M. Hammersley e D.
C. Handscomb, na obra Monte Carlo Methods, de 1964, afirmam que o nome do método foi dado durante
o Projeto Manhattan, homenageando um dos distritos de Mônaco famoso por seus cassinos. [8]
3
No entanto, existe um inconveniente no uso dos métodos de Monte Carlo e Dinâmica
Molecular no que se refere ao aprisionamento do sistema em mínimos locais do espaço
de fase do sistema. Isso significa que em muitas temperaturas de interesse, estados de
mais baixa energia se tornam inacessíveis pelos métodos devido às altas barreiras de
energias, as quais o sistema não pode transpor. Portanto, quanto maior a complexidade do
sistema, mais rugosa se torna sua superfície de energia potencial, dificultando o acesso às
configurações representativas no cálculo de processos termodinâmicos. Nesse ponto, uma
estratégia envolvendo trocas de réplicas simuladas em diferentes temperaturas (parallel
tempering) tem uma interessante utilidade, uma vez que, em princípio, diferentes regiões
do espaço de fase podem ser acessadas a partir dos critérios de troca num tipo genérico
de ensemble considerado.
Na próxima seção, serão apresentadas as bases do método Monte Carlo, discutindo
apenas seus aspectos funcionais no sentido de tornar a discussão envolvendo o algoritmo
parallel tempering mais adequada.
2 Fundamentos do Método Monte Carlo
Considerando o ensemble canônico (NV T ), em que o número de partículas (N), o
volume (V ) e a temperatura (T ) são constantes, o qual também é conhecido como ensemble
NV T , a energia do sistema pode ser calculada com base no Hamiltoniano H(qi, pi). O
Hamiltoniano é definido pelos momentos e posições das partículas do sistema, em que
H(Ri,pi) = K(pi) + V (Ri) =
N∑
i=1
p2i
2m
+ V (Ri), ⇒

Ri ≡ Ri(R1,R2, . . . ,RN)
pi ≡ pi(p1,p2, . . . ,pN)
m ≡ mk for k = 1, 2, . . . , N
(1)
onde, as partículas possuem posições Ri(R1,R2, . . . ,RN), momentos pi(p1,p2, . . . ,pN)
e massas mk para k = 1, 2, . . . , N , os quais definem os termos de energia cinética (K(pi))
e potencial de interação das partículas (V (Ri)). Desse modo, é possível calcular a média
de uma dada propriedade do sistema a partir de uma relação entre o fator de Boltzmann
e a função de partição canônica, ou seja,
〈A(Ri,pi)〉T =
∫
1
∫
2
. . .
∫
N
A(Ri,pi) e
−βH(Ri,pi)dRi dpi∫
1
∫
2
. . .
∫
N
e−βH(Ri,pi)dRi dpi
, onde β =
1
kBT
. (2)
Utilizando o Hamiltoniano definido na equação 1, em que K(pi) e V (Ri) são respecti-
vamente dependentes apenas dos momentos e das posições, a energia média do sistema
〈E(Ri,pi)〉 = 〈A(Ri,pi)〉T pode ser calculada a partir das contribuições individuais de
K(pi) e V (Ri) utilizando a equação :
〈H(Ri,pi)〉T = 〈K(pi)〉T + 〈V (Ri)〉T
=
∫
1
∫
2
. . .
∫
N
K(pi) e
−βK(pi)dpi∫
1
∫
2
. . .
∫
N
e−βK(pi) dpi
+
∫
1
∫
2
. . .
∫
N
V (Ri) e
−βV (Ri)dRi∫
1
∫
2
. . .
∫
N
e−βV (Ri) dRi
.(3)
4
Uma vez que o primeiro termo da equação 3, relativo à média da energia cinética,
depende apenas dos momentos, sua solução pode ser obtida utilizando um modelo de gás
ideal (onde não há interação entre as partículas do sistema), a qual tem sua resolução
discutida em vários livros texto. A partir das integrais envolvidas no segundo, o cálculo
de várias propriedades termodinâmicas dependem do acesso às configurações espaciais
possíveis para resolução dessas integrais multidimensionais. Consequentemente, uma vez
determinado o conjunto de todas as configurações espaciais possíveis, a função de partição
canônica (Z) pode ser obtida de forma exata e, assim, levando ao cálculo de qualquer
propriedade termodinâmica. Logo, uma vez que
〈H(Ri)〉T =
∫
1
∫
2
. . .
∫
N
H(Ri) e
−βV (Ri)dRi∫
1
∫
2
. . .
∫
N
e−βV (Ri) dRi
. ⇒ Z =
∫
1
∫
2
. . .
∫
N
e−βH(Ri) dRi, (4)
temos que Z é justamente a função de partição canônica configuracional. No entanto,
é impossível determinar todas as configurações espaciais possíveis (consequentemente,
uma distribuição completa de Z) em um sistema de N partículas mais complexas (como
moléculas), e é nesse ponto que o método Monte Carlo se mostra extremamente útil a
partir do trabalho de Metropolis et al. [9].
A ideia por trás do algoritmo Metrópolis consiste, basicamente, no uso de um meca-
nismo de evolução probabilística na qual os valores das propriedades médias podem ser
obtidas sem a necessidade de conhecimento de todo espaço configuracional. Isso signi-
fica que as propriedades são calculadas sobre um espaço amostral, gerado exatamente a
partir do uso do algoritmo Metrópolis seguindo algumas regras para escolha das configu-
rações. Desse modo, a convergência na formação de um conjunto satisfatório é acelerada
a partir do que é chamado tradicionalmente de importance sampling(amostragem por
importância). Como consequência disso, a primeira aproximação ocorre na função de
partição canônica, a qual pode ser determinada pela soma de ν configurações acessadas
no processo, ou seja:
Z =
∫
1
∫
2
. . .
∫
N
e−βV (Ri) dRi ≈
ν∑
ν=1
e−βVν(Ri). (5)
A partir dessa aproximação, é possível determinar a probabilidade de observação para
cada configuração por
P (Rνi ) =
e−βVν(Ri)
ν∑
ν=1
e−βVν(Ri)
. (6)
Portanto, cada configuração ν é selecionada pelo algoritmo Metrópolis utilizando a amos-
tragem por importância a fim de concentrar as configurações nas regiões mais densas do
espaço de fase, ou seja, termodinamicamente mais prováveis. Uma vez que o processo de
escolha sucessiva de cada configuração se constitui uma evolução probabilística, não há
nenhuma relação dinâmica temporal no processo. Isso significa que as configurações são
discretas, mas não sucessivas no tempo, justamente porque não há evolução temporal no
método Monte Carlo. Essa característica intrínseca do método se constitui um caso de
processo estocástico, no entanto, o algoritmo Metrópolis avança na formação do elencado
5
de probabilidades conforme a distribuição de Boltzmann por um processo de memória
markoviana. Essa estratégia é definida como o uso do método Monte Carlo de Cadeia de
Markov, desenvolvido pelo matemático Andrei Andreyevich Markov, em que os estados
anteriores ao atual não são considerados na predição dos estados seguintes.
Vejamos a seguir como esse processo é feito. Inicialmente, temos a transição do estado
Ain inicial para o estado Bse seguinte:
Ain ≡ {Rνi } Prob. de Transição−−−−−−−−−−−−−−→
pi(Ain→Bse)
Bse ≡ {Rν′i }, (7)
onde, a mudança de estado ocorre a partir da mudança da configuração espacial das par-
tículas de {Rνi } para {Rν′i } com probabilidade dada por pi(Ain → Bse). A probabilidade
de transição pode ser decomposta em
pi(Ain → Bse) = ξ(Ain → Bse) × acc(Ain → Bse)
{
ξ(Ain → Bse) ≡ Prob. do processo ocorrer
acc(Ain → Bse) ≡ Prob. de aceitação
Nesse ponto, a condição de balanço detalhado precisa ser obedecida, no qual garante que
o sistema parmaneça em equilíbrio, onde:
Ase ≡ {Rνi } Prob. de Transição←−−−−−−−−−−−−−−−−
pi(Bin→Ase)
Bin ≡ {Rν′}, (8)
e assim,
pi(Bin → Ase) = ξ(Bin → Ase) × acc(Bin → Ase)
{
ξ(Bin → Ase) ≡ Prob. do processo ocorrer
acc(Bin → Ase) ≡ Prob. de aceitação
Tomando as probabilidades dos estados Ain e Bse definidas na distribuição de Boltzmann,
temos
PA(R
ν
i ) =
e−βVν(Ri)
ν∑
ν=1
e−βVν(Ri)
, PB(R
ν′
i ) =
e−βVν′(Ri)
ν∑
ν=1
e−βVν(Ri)
. (9)
Portanto, a aplicação do produto das probabilidades de encontrar os respectivos estados
e suas transições, pode ser dada por
PApi(Ain → Bse) = PBpi(Bin → Ase) (10)
pi(Ain → Bse)
pi(Bin → Ase) =
PB
PA
pi(Ain → Bse)
pi(Bin → Ase) =
e−βVν′(Ri)
��
��
��
�ν∑
ν=1
e−βVν(Ri)
×
 e
−βVν(Ri)
��
��
��
�ν∑
ν=1
e−βVν(Ri)

−1
pi(Ain → Bse)
pi(Bin → Ase) = e
−βHν′(Ri) × eβVν(Ri) = e−β[Vν′(Ri)−Vν(Ri)]
pi(Ain → Bse)
pi(Bin → Ase) = e
−β[Vν′(Ri)−Vν(Ri)]. (11)
6
Para o lado esquerdo, aplicando a decomposição na probabilidade de transição:
((((
(((ξ(Ain → Bse)
((((
(((ξ(Bin → Ase) ×
acc(Ain → Bse)
ξ(Bin → Ase) = e
−β[Vν′(Ri)−Vν(Ri)]
acc(Ain → Bse)
ξ(Bin → Ase) = e
−β[Vν′(Ri)−Vν(Ri)] (12)
A partir da equação 12 se aplica o critério de Metrópolis, de modo que:
acc(Ain → Bse) ⇒

PB
PA
, se PB < PA
1, se PB ≥ PA
A partir dos critérios de Metrópolis definidos, o protocolo seguido pelo algoritmo é
disposto da seguinte forma:
• Gerar uma configuração aleatória c2 a partir de uma configuração prévia c1.
• Se a energia de c2 for menor do que c1, a configuração é incluída no conjunto de
amostragem.
• Caso contrário, um número aleatório entre 0 e 1 é gerado. Se esse número for menor
que PB
PA
, aceita-se a configuração c2 que será utilizada como configuração c1 para um
novo movimento.
• Caso contrário, a configuração c1 permanece para uma outra tentativa.
• O procedimento se repete em função de um critério de parada, de modo que, a partir
do conjunto de configurações formado as propriedades termodinâmicas podem ser
calculadas.
3 Fundamentos de Dinâmica Molecular: Formalismo Clássico
Ao contrário do Monte Carlo, no método de Dinâmica Molecular (MD) o tempo é
um parâmetro explícito no formalismo envolvido, uma vez que o espaço configuracional é
determinado pela variação da posição de cada partícula do sistema prevista por equações
clássicas de movimento. A ideia central no método de MD consiste em determinar a
evolução temporal através da força que as partículas aplicam umas sobre as outras. Desse
modo, uma vez que
{~RNi } =

~R1 = ~X1 + ~Y1 + ~Z1
~R2 = ~X2 + ~Y2 + ~Z2
~R3 = ~X3 + ~Y3 + ~Z3
~R4 = ~X4 + ~Y4 + ~Z4
...
~RN = ~XN + ~YN + ~ZN
{~vNi } =

~v1 =
d~X1
dt
+
d~Y1
dt
+
d~Z1
dt
~v2 =
d~X2
dt
+
d~Y2
dt
+
d~Z2
dt
~v3 =
d~X3
dt
+
d~Y3
dt
+
d~Z3
dt
~v4 =
d~X4
dt
+
d~Y4
dt
+
d~Z4
dt
...
~vN =
d~XN
dt
+
d~YN
dt
+
d~ZN
dt
7
com {~RNi } e {~vNi } sendo o conjunto de posições e velocidades, respectivamente, das
partículas do sistema, é possível determinar a evolução para mínimos de energia no seu
espaço de fase correspondente. Portanto, permitindo calcular as médias no ensemble.
A base fundamental do método de MD consiste na resolução das equações de movimento
a partir do formalismo da mecânica newtoniana. A segunda lei de Newton prevê que
N∑
i=1
~Fi = −∇~Ri
N∑
i=1
V ({~RNi }) = −∇~Ri
N∑
j=1
V (|~Ri − ~Rj|), (13)
de modo que a força de interação entre as partículas é dada pelo gradiente negativo
do potencial de interação clássico entre elas∗∗. Desse modo, tendo obtido a força entre
as partículas por meio do potencial, é possível relacionar com o conjunto de massas e
acelerações da seguinte forma:
N∑
i=1
~Fi =
N∑
i=1
mi~ai (14)
N∑
i=1
~Fi =
N∑
i=1
mi
d2{~RNi }
dt2
= m1
d2 ~R1
dt2
+m2
d2 ~R2
dt2
+m3
d2 ~R3
dt2
+ . . .+mN
d2 ~RN
dt2
−∇~Ri
N∑
j=1
V (|~Ri − ~Rj|) = m1d
2 ~R1
dt2
+m2
d2 ~R2
dt2
+m3
d2 ~R3
dt2
+ . . .+mN
d2 ~RN
dt2
Logo,
−∇~Ri
N∑
j=1
V (|~Ri − ~Rj|) =
N∑
i=1
mi
d2{~RNi }
dt2
(15)
A parte esquerda da equação 15 pode ser obtida a partir do potencial de interação
previamente parametrizado, como por exemplo, utilizando aproximações empíricas. Para
o termo da direita, a expansão das posições e velocidades no tempo pode ser obtida
utilizando métodos de aproximação numérica, definindo pequenos intervalos de tempo δt
no qual o termo será integrado. [2,4,11,15] A próxima seção tratará desse tema de modo
introdutório.
3.1 Integração Numérica: Método das Diferenças Finitas
Para integração do termo d
2{~RNi }
dt2
diversos algoritmos foram propostos com base no
método de diferenças finitas, o qual utiliza uma expansão em séries de Taylor para solução
de equações diferenciais. Tomando uma função g(x), sua expansão em intervalos ∆x na
∗∗O potencial de interação clássico é determinado a partir da parametrização das constantes de força
associadas à interação de pares no sistema. Em geral, são potenciais ajustados empiricamente, ou ainda,
parametrizados a partir de cálculos quânticos. Entretanto, vale ressaltar que nesse potencial podem ser
incluídas todas as interações de diferentes "naturezas"(intermoleculares, interação de Coulomb, estira-
mento de ligações, variações angular e torções de diedro), sendo assim chamado de Campo de Força do
sistema - popularmente chamado de Force Field.
8
série de Taylor é definida como:
g(x+ ∆x) =
N∑
n=0
g(n)
n!
(x−∆x)n
= g +
dg
dx
(x−∆x)1
1!
+
d2g
dx2
(x−∆x)2
2!
+
d3g
dx3
(x−∆x)33!
+ . . .
As difereças finitas podem ser obtidas da expansão e retrocesso em intervalos de h e
−h, respectivamente, ou seja:
g(x+ h) = g +
dg
dx
h+
d2g
dx2
h2
2
+
d3g
dx3
h3
6
+ ϑ(h4) + . . .+ ϑ(hn)
g(x− h) = g − dg
dx
h+
d2g
dx2
h2
2
− d
3g
dx3
h3
6
+ ϑ(h4) + . . .+ ϑ(hn)
onde, ϑ(hn) são os termos de ordens superiores. Desse modo, considerando o incremento
h = t+ ∆t, e aplicando o método de diferenças finitas na expansão de ~Ri na equação 15,
alguns algoritmos foram propostos a partir de algumas técnicas formais.
• Algoritmo de Verlet
O algoritmo de Verlet, proposto por Loup Verlet em 1960, é um dos mais utilizados
entre os métodos de diferenças finitas aplicados em Dinâmica Molecular, no qual expande-
se a posição ~Ri da partícula para os instantes t+∆t e t−∆t nas séries de Taylor, truncando
a série a partir do quarto termo, ou seja,
~Ri(t+Δt) = ~Ri(t) +
d~Ri(t)
dt
Δt+
d2 ~Ri(t)
dt2
Δt2
2
+ . . . ,
~Ri(t−Δt) = ~Ri(t)− d
~Ri(t)
dt
Δt+
d2 ~Ri(t)
dt2
Δt2
2
− . . . ,
~Ri(t+Δt) = ~Ri(t) + ~vi(t)Δt+ ~ai(t)
Δt2
2
+ . . .
~Ri(t−Δt) = ~Ri(t)− ~vi(t)Δt+ ~ai(t)Δt
2
2
− . . .
Somando as equações ~Ri(t+ ∆t) e ~Ri(t−∆t), temos:
~Ri(t+ ∆t) + ~Ri(t−∆t) = 2~Ri(t) + ~ai(t)∆t2 + . . .
⇓
~Ri(t+ ∆t) = 2~Ri(t)− ~Ri(t−∆t) + ~ai(t)∆t2 + . . . , (16)
onde, uma vez que,
~ai(t) =
Fi
mi
, logo, ~Ri(t+ ∆t) = 2~Ri(t)− ~Ri(t−∆t) +
(
Fi
mi
)
∆t2 + . . . , (17)
e assim,
~Ri(t+ ∆t) = 2~Ri(t)− ~Ri(t−∆t)−
[
1
mi
∇~Ri
N∑
j=1
V (|~Ri − ~Rj|)
]
∆t2. (18)
A posição da partícula i em t+∆t depende da posição no tempo anterior a t e na diferencial
do potencial de interação com relação à posição previamente calculado.
9
Para as velocidades, tomamos a diferença em ~Ri(t+ ∆t) - ~Ri(t−∆t), onde,
~Ri(t+ ∆t)− ~Ri(t−∆t) = 2~vi(t)∆t
⇓
~vi(t) =
~Ri(t+ ∆t)− ~Ri(t−∆t)
2∆t
(19)
Observamos a partir das operações realizadas na determinação de ~Ri(t+ ∆t) e ~vi(t) que
o erro das posições com o algoritmo de Verlet é da ordem de (∆t)4, e das velocidades
de (∆t)3. Desse modo, a integração numérica proporciona a evolução do sistema para
mínimos locais no espaço de fase; para sistemas mais simples, diversas propriedades ter-
modinâmicas podem então ser calculadas a partir dessa amostragem das configurações no
tempo, as quais são construídas pelo processo dinâmico.
4 Replica Exchange Method (REM): Parallel Tempering
Com o aumento da complexidade dos sistemas uma infinidade de mínimos locais se
formam no espaço de fase. Nesse contexto, uma das dificuldades dos métodos de Monte
Carlo e Dinâmica Molecular tradicionais se refere ao aprisionamento do sistema em al-
gum desses mínimos locais, impedindo a evolução no procedimento de amostragem para
regiões com maior estabilidade termodinâmica. Em geral, isso ocorre devido às altas bar-
reiras de energia que separam esses mínimos e que não podem ser vencidas com muitas
temperaturas de interesse − como a temperatura ambiente.
O método de troca de réplicas (Replica Exchange Method − REM) foi desenvolvido
primeiro por Swendsen e Wang com o método Monte Carlo, [16] e então, sendo posterior-
mente melhorado por Hukushima e Nemoto, [17] cuja intenção era contornar o problema
de transposição das barreiras de energia entre mínimos locais. Esse procedimento con-
siste no uso de M diferentes cópias, ou réplicas, do sistema que não interagem entre si,
formando assim o que é definido como generalized ensemble. Cada réplica corresponde a
uma cópia original de um ensemble específico − e aqui será utilizado o ensemble canônico
por descrever muitas situações experimentais − onde cada uma delas está imersa num
termostato diferente, definindo M temperaturas Tm(m = 1, 2, . . . ,M) para M réplicas,
ou seja, como temperaturas paralelas (parallel tempering). De modo simbólico:
Generalized ensemble

Coordenadas relativas Ri,pi :
x1(R,p) ⇒ T1
x2(R,p) ⇒ T2
. . .
xi(R,p) ⇒ Tm
⇒

Estados e temperaturas :
X1(x1) ⇒ T1
X2(x2) ⇒ T2
. . .
Xi(xi) ⇒ Tm
Portanto, há uma correspondência entre as réplicas e a temperatura, de modo que existe
uma identificação mútua entre elas que será útil no critério de troca, ou seja, como critério
de permutação entre réplicas com suas temperaturas, onde i(i = 1, 2, . . . ,M) coorresponde
à identificação para as réplicas e m(m = 1, 2, . . . ,M) para as temperaturas.
Logo, podemos representar os estados no generalized ensemble por
X = (x[i(1)]1 ,x
[i(2)]
2 , . . . ,x
[i(M)]
M )︸ ︷︷ ︸
Replica exchange
= (x
[1]
m(1),x
[2]
m(2), . . . ,x
[M ]
m(M))︸ ︷︷ ︸
Temperature exchange
, (20)
10
onde, o estado X é dado por M réplicas, de modo que em cada uma delas as coordenadas
no espaço de fase das N partículas são dadas por R[i] e p[i] em uma das temperaturas em
Tm, identificadas pelo índice [i] como notação para as réplicas.
Definiremos agora a troca de um par de réplicas no generalized ensemble identificadas
por i e j, imersas em termostatos de temperaturas Tm e Tn, respectivamente, o que resulta
na troca de dois estados como
X = (. . . ,x[i]m, . . . ,x[j]n , . . .) → X′ = (. . . ,x[j]′m , . . . ,x[i]′n , . . .). (21)
Logo, podemos escrever as trocas sob a perspectiva das réplicas e das temperaturas da
seguinte forma:
Replicas
{
x
[i]
m ≡ (R[i],p[i])m Swap−−−−−−−→ x[j]′m ≡ (R[j]′,p[j]′)m
x
[j]
n ≡ (R[j],p[j])n Swap
′−−−−−−−→ x[i]′n ≡ (R[i]′,p[i]′)n
Uma vez que não há interação entre as M réplicas, o peso estatístico do estado X
no ensemble generalizado é dado pelo produto do fator de Boltzmann para cada réplica,
dependendo de como a troca será efetuada. Desse modo, temos que
W (X) =
M∏
i=1
e−βmH(R
[i],p[i]) = exp
{
M∑
i=1
−βmH(R[i],p[i])
}
. (22)
Omitindo os argumentos R[i],p[i] de H(R[i],p[i]), consideraremos a troca das réplicas
expressas da seguinte forma:
exp
−β1H [1(1)] + . . .− βmH [i(m)] − βnH [j(n)]︸ ︷︷ ︸
troca
+ . . .− βnH [M(M)]

...
exp
−β1H [1(1)] + . . .− βmH [j(n)] − βnH [i(m)]︸ ︷︷ ︸
Réplicas trocadas
+ . . .− βnH [M(M)]

Desse modo, temos a troca de temperatura de duas réplicas i e j pela transição de estado
dada por
W (X) ⇒ W (X′), Com probabilidade de troca w(X→ X′)
Tanto no método Monte Carlo quanto em Dinâmica Molecular, o critério para troca
das temperaturas das réplicas será dado pela condição de balanço detalhado, uma vez que
a distribuição de equilíbrio nessa troca será dada pelo produto do fator de Boltzmann
com a probabilidade de troca. Logo,
W (X)w(X→ X′) = W (X′)w(X′ → X)
w(X→ X′)
w(X′ → X) =
W (X′)
W (X)
, (23)
Todavia, o formalismo entre ambos os métodos se diferem por alguns detalhes. Veremos
esses pontos a seguir.
11
• Parallel Tempering Monte Carlo (PTMC)
A primeira característica no uso do PTMC segundo a proposta de Swendsen e Hu-
kushima, está na troca das réplicas, onde somente a parte configuracional da equação 3
será considerada. Isso porque a integral referente à energia cinética possui solução ana-
lítica para cada temperatura. Desse modo, o hamiltoniano do produtório da equação 22
para o fator de Boltzmann considera apenas a contribuição da interação entre as partícu-
las, de modo que se aplica a probabilidade Pβm(R) (sistema com conjunto k de partículas)
para ocorrência do estado X em uma seleção aleatória,ou seja,
W (X) =
M∏
i=1
e−βmV (R
[i]) onde Pβm(R) =
e−βVk(R)∑
k
e−βVk(R)
(24)
A partir dessas definições, a distribuição de probabilidades das M réplicas (Ωβ) é
dada a partir das temperaturas elencadas em m(m = 1, 2, . . . ,M) pelo produtório das
probabilidades de cada estado onde propomos a troca:
Ωβ =
M∏
i=1
Pβm(R) = [P
[1]
β1
(R)× · · · × P [i]βm(R)× P
[j]
βn
(R)︸ ︷︷ ︸
Troca
× . . .× P [M ]βM (R)] (25)
Portanto, no uso da condição de balanço detalhado, considerandoas transições entre
réplicas em equilíbrio dadas por X → X′ e X′ → X com suas respectivas probabilidades
de troca já definidas anteriormente, temos que
M∏
i=1
Pβm(R)w(X→ X′) =
M∏
i=1
Pβn(R)w(X′ → X)[
. . .× P [i]βm(R)× P
[j]
βn
(R)× . . .
]
w(X→ X′) = [. . .× P [j]βm(R)× P
[i]
βn
(R)× . . .]w(X′ → X).
De modo similar ao método Monte Carlo, onde a probabilidade de transição configuracio-
nal é dada pelo produto da probabilidade de ocorrência do processo com a probabilidade
de aceitação, em PTMC a probabilidade de troca é decomposta no produto da probabi-
lidade da troca ocorrer (k(X → X′)) com a probabilidade de aceitação (ACC(X → X′).
Portanto, pela condição de balanço detalhado temos:
M∏
i=1
Pβm(R)k(X→ X′)ACC(X→ X′) =
M∏
i=1
Pβn(R)k(Xn → X′m)ACC(Xn → X′m)
k(X→ X′)ACC(X→ X′)
k(X′ → X′)ACC(X′ → X) =
M∏
i=1
Pβn(R)
M∏
i=1
Pβm(R)
���
���k(X→ X′)ACC(X→ X′)
���
���k(X′ → X′)ACC(X′ → X′) =
���
���
�
P
[1]
β1
(R)× . . . P [i]βn(R) × P
[j]
βm
(R)���
���
��
. . .× P [M ]βM (R)
���
���
�
P
[1]
β1
(R)× . . . P [i]βm(R) × P
[j]
βn
(R)���
���
��
. . .× P [M ]βM (R)
ACC(Xm → X′n)
ACC(Xn → X′m)
=
P
[i]
βn
(R) × P [j]βm(R)
P
[i]
βm
(R)× P [j]βm(R)
, (26)
12
onde, os termos do lado direito correspondem às únicas réplicas onde ocorreram as trocas.
Portanto, pela definição da probabilidade de ocorrência de dado estado a partir do fator
de Boltzmann, temos que:
ACC(X→ X′)
ACC(X′ → X) =
e−βnV
[i](R)∑
k
���
�
e−βkVk
×
 e−βmV [j](R)∑
k
���
��
e−βkVk(R)

−1
× e
−βmV [i](R)∑
k
���
��
e−βkV (R)
×
 e−βnV [j](R)∑
k
���
��
e−βkV (R)

−1
= e−βnV
[i](R) × e−βmV [j](R) × eβmV [i](R) × eβnV [j](R)︸ ︷︷ ︸
Agrupando
= exp
[−βnV [i](R)− βmV [j](R) + βmV [i](R) + βnV [j](R)]
= exp
{−βn [V [i](R)− V [j](R)]+ βm [V [i](R)− V [j](R)}
= exp
{
(βm − βn)
[
V [i](R)− V [j](R)]}
Considerando (βm − βn)
[
V [i](R)− V [j](R)] = ∆, temos
ACC(X→ X′)
ACC(X′ → X) = e
Δ (27)
A equação 27 segue os mesmos critérios de troca que a equação 13 para aceitação
de um movimento segundo o critério de Metrópolis. Após a discussão do método de
parallel tempering no contexto da Dinâmica Molecular, alguns pontos importantes serão
detalhados uma vez que uma equação similar será obtida.
• Parallel Tempering Molecular Dynamics (PTMD)
O uso do método de Dinâmica Molecular no contexto do cálculo da trajetória em
réplicas paralelas foi proposto de Sugita e Okamoto. [18] Nessa abordagem, os momen-
tos lineares das partículas nas réplicas trocadas devems er considerados, sendo, então,
ajustados pela seguinte relação:
p[i]
′
=
√
Tn
Tm
p[i], p[j]
′
=
√
Tm
Tn
p[j]. (28)
No que tange ao protocolo de troca das réplicas, à relação entre a probabilidade de
ocorrência de um determinado estado e a condição de balanço detalhado, o procedimento
é o mesmo daquele apresentado no PTMC. Vale também que, exceto para a transição de
estados envolvida na troca das réplicas envolvidas, as demais réplicas permanecem nas
mesmas temperaturas. No entanto, a principal diferença com relação ao PTMC, além da
dinâmica em cada réplica, na qual se aplica o algoritmo de integração na propagação das
velocidades e posições das partículas no tempo, consiste no uso do hamiltoniano no fator
de Boltzmann e na função de partição da função de probabilidade. Ou seja, a distribuição
de probabilidades e a probabilidade de obtenção de um estado são expressas em função
do hamiltoniano do sistema:
W (X) =
M∏
i=1
e−βmH(R
[i],p[i]) onde Pβm(R,p) =
e−βHk(R,p
[i])∑
k
e−βHk(R,p
[i])
. (29)
13
O protocolo de desenvolvimento com base na troca de duas réplicas segue o mesmo forma-
lismo aplicado ao PTMC. Portanto, seguiremos a partir da aplicação do balanço detalhado
no agrupamento das funções exponenciais:
ACC(X→ X′)
ACC(X′ → X) =
P
[i]
βn
(R,p) × P [j]βm(R,p)
P
[i]
βm
(R,p)× P [j]βm(R,p)
=
e−βnH(R
[i],p[i])∑
k
((((
((
e−βkHk(R,p)
×
e−βmH(R[j],p[j]
′
)(R)∑
k
((((
((
e−βkHk(R,p)

−1
× e
−βmH(R[i],p[i])∑
k
((((
((
e−βkHk(R,p)
×
 e−βnH(R[j],p[j])∑
k
((((
((
e−βkHk(R,p)

−1
= exp
[
−βnH(R[i],p[i]′)− βmH(R[j],p[j]′) + βmH(R[i],p[i]) + βnH(R[j],p[j])
]
Uma vez que H(R[i],p[i]) = K(p[i]) + V (R[i]),
ACC(X→ X′)
ACC(X′ → X) = exp
{
−βn
[
K(p[i]
′
) + V (R[i])
]
− βm
[
K(p[j]
′
) + V (R[j])
]}
+ exp
{
βm
[
K(p[i]) + V (R[i])
]
+ βn
[
K(p[j]) + V (R[j])
]}
.
Organizando em blocos de K(p) e V (R), temos:
ACC(X→ X′)
ACC(X′ → X) =

exp
[
−βnK(p[i]′)− βmK(p[j]′) + βmK(p[i]) + βnK(p[j])
]
+
exp
[−βnV (R[i])− βmV (R[j]) + βmV (R[i]) + βnV (R[j])] .
O primeiro bloco de energias cinéticas é resolvido utilizando as relações definidas nas
equações em 28, e ainda, a partir das definições de energia cinética em termos de momento
linear, considerando βm = 1/kBTm, temos:
exp
[
−βnK(p[i]′)− βmK(p[j]′) + βmK(p[i]) + βnK(p[j])
]
exp
[
−βn (p
[i]′)2
2m
− βm (p
[j]′)2
2m
+ βmK(p
[i]) + βnK(p
[j])
]
exp
[
−βn
(
Tn
Tm
)
(p[i])2
2m
− βm
(
Tm
Tn
)
(p[j])2
2m
+ βmK(p
[i]) + βnK(p
[j])
]
exp
[
− 1
kB��Tn
(
��Tn
Tm
)
(p[i])2
2m
− 1
kB��Tm
(
��Tm
Tn
)
(p[j])2
2m
+ βmK(p
[i]) + βnK(p
[j])
]
exp
[
− 1
kBTm
(p[i])2
2m
− 1
kBTn
(p[j])2
2m
+ βmK(p
[i]) + βnK(p
[j])
]
Observemos a conversão de K(p[i]′)→ K(p[i]) e K(p[j]′)→ K(p[j]), de modo que:
exp
[
���
����−βmK(p[i])−������βnK(p[j]) +�����
�
βmK(p
[i]) +���
���βnK(p
[j])
]
,
14
e assim, o primeiro termo se anula
exp
[
−βnK(p[i]′)− βmK(p[j]′) + βmK(p[i]) + βnK(p[j])
]
= 0, (30)
de modo que:
ACC(X→ X′)
ACC(X′ → X) = exp
[−βnV (R[i])− βmV (R[j]) + βmV (R[i]) + βnV (R[j])]
= exp
{−βn [V [i](R)− V [j](R)]+ βm [V [i](R)− V [j](R)}
= exp
{
(βm − βn)
[
V [i](R)− V [j](R)]}
ACC(X→ X′)
ACC(X′ → X) = e
Δ (31)
Desse modo, as equações em função da exponencial definida por ∆ para PTMC e
PTMD, equações 27 e 31, respectivamente, podem ser satisfeitas pelo critério de Metró-
polis:
ACC(X→ X′) ⇒

eΔ, se ∆ > 0
1, se ∆ ≤ 0
(32)
Às equações definidas em 27 e 31, nos métodos de PTMC e PTMD, respectivamente,
se aplica o critério de Metrópolis definido em 32, no qual se estabelece as operações de
aceitação ou rejeição para a troca de duas réplicas. Cada réplica possui uma temperatura
fixa (e diferentes entre si) nos seus respectivos ensembles canônicos, onde os sistemas
são simulados de modo simultaneo e independente. Desse modo, um par de réplicas
vizinhas são trocadas pelo critério de aceitação a partir do agoritmo de Metrópolis. Uma
vez que o método otimiza a amostragem do espaço configuracional do sistema, uma das
principais vantagens dos métodos consiste no cálculo das observáveis termodinâmicas em
uma dada temperatura intermediária, que necessariamente não faz parte do conjunto de
temperaturas definidas para as réplicas. O weighted histogram analysis method (WHAM)
é dos mais aceitos nesse propósito, e maiores detalhes podem ser acessados na publicação
de S. Kumar [19].
Referências
[1] J.-P Hansen. Molecular-Dynamics Simulation os Statistical-Mechanical Systems.
Amsterdam, 1986.
[2] Terrel L. Hill. Statistical Mechanics: Principles and Selected Applications. New York,
Addison-Wesley Publishing Company, 1987.
[3] Donald Allan McQuarrie. Physical Chemistry: a Molecular Approach. Sausalito,
University Science Books, 1997.
[4] Donald Allan McQuarrie. Statistical Thermodynamics. Sausalito, University Science
Books, 2000.
15
[5] Y. Prigogine. O Fim das Certezas: O Tempo, Caos e as Leis da Natureza. Editora
Unesp, São Paulo, 2011.
[6] I.Müller. A History of Thermodynamics: A Doctrine of Energy and Entropy.
Springer-Verlag Berlin Heidelberg, 2007.
[7] R. E. Dickerson. Molecular Thermodynamics. W. A. Benjamin, Inc., New York,
1969.
[8] W. Reiher. Hammersley, j. m., d. c. handscomb: Monte carlo methods. Biometrische
Zeitschrift, 8(3), 1966.
[9] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.
Equation of state calculations by fast computing machines. The Journal of Chemical
Physics, 21(6), 1953.
[10] B. J. Alder and T. E. Wainwright. Studies in molecular dynamics. i. general method.
The Journal of Chemical Physics, 31(2), 1959.
[11] Michael P. Allen and D. J. Tildesley. Compute Simulation of Liquids. Oxford, Oxford
University Press, 1987. cap 3.
[12] Daan Frenkel and Berend Smith. Understanding Molecular Simulation: From Algo-
rithms to Applications. 2a ed. San Diego, Academic Press, 2001.
[13] J. M. Haile. Molecular Dynamics Simulation: Elementary Methods. 1a ed. New York,
John Wiley & Sons Inc., 1992.
[14] Loup Verlet. Computer "experiments"on classical fluids. i. thermodynamical proper-
ties of lennard-jones molecules. Physical Review, 159:98–103, 1967.
[15] D.C. Rapaport. The Art of Molecular Dynamics Simulation. 2a ed. Cambridge,
Cambridge University Press, 2004.
[16] R. H. Swendsen and J.-S. Wang. Replica monte carlo simulation of spin-glasses.
Physical Review Letters, 57:2607–2609, 1986.
[17] K. Hukushima and K. Nemoto. Exchange monte carlo method and application to
spin glass simulations. Journal of the Physical Society of Japan, 65(6):1604–1608,
1996.
[18] Y. Sugita and Y. Okamoto. Replica-exchange molecular dynamics method for protein
folding. Chemical Physics Letters, 314(1–2):141 – 151, 1999.
[19] S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman. The
weighted histogram analysis method for free-energy calculations on biomolecules. i.
the method. Journal of Computational Chemistry, 13(8):1011–1021, 1992.
16