Replica Exchange
Introdução ao Método de : 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.
[email protected];
[email protected] 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 desenvolvimento 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 vapor. 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 Henderson é 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:
[email protected];
[email protected] 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 formalí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,S pode ser admitida como rebento direto da evolução do conceito de calor, se aproximando das protodefiniçõ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 segunda 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 observaçã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 declara: 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 cientí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 Dickerson 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 boltzmanniana 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. S 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 𝑛𝑖 partículas com energia 𝜖𝑖 dispostas de uma determinada forma (estado), podendo ser contadas em um ou mais estados de energia definidos pelo número de degenerecência 𝑔𝑖 . 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 no ensemble, no qual cada microestado é representativo do sistema macroscópico em questão. Portanto, a partir de determinadas operações (previamente definidas no formalismo da termodinâmica estatística) aplicadas sobre uma função que descreve o ensemble, chamada função de partição 𝒵 , 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 microestado 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 partículas distinguíveis (boltzons) definidos pela estatística de Maxwell-Boltzmann. Nesses sistemas, o grau de degenerecência 𝑔𝑖 para cada estado é muito grande com relação ao número de partículas 𝑛𝑖 (distinguíveis), de modo que, 𝑔𝑖 >> 𝑛𝑖 , possibilitando a determinação exata do ensemble mais simples conhecido como microcanônico (𝑁 𝑉 𝐸 )¶ . 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 possibilidade 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 colaboradores (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 consideraçã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, 𝑁 , o volume do sistema, 𝑉 , e a energia, 𝐸, 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 (𝑁 𝑉 𝑇 ), em que o número de partículas (𝑁 ), o volume (𝑉 ) e a temperatura (𝑇 ) são constantes, o qual também é conhecido como ensemble 𝑁 𝑉 𝑇 , a energia do sistema pode ser calculada com base no Hamiltoniano 𝐻(𝑞𝑖 , 𝑝𝑖 ). O Hamiltoniano é definido pelos momentos e posições das partículas do sistema, em que ⎧ ⎪ 𝑁 ⎨R𝑖 ≡ R𝑖 (R1 , R2 , . . . , R𝑁 ) ∑︁ p2𝑖 (1) + 𝑉 (R𝑖 ), ⇒ p𝑖 ≡ p𝑖 (p1 , p2 , . . . , p𝑁 ) 𝐻(R𝑖 , p𝑖 ) = 𝐾(p𝑖 ) + 𝑉 (R𝑖 ) = ⎪ 2𝑚 ⎩ 𝑖=1 𝑚 ≡ 𝑚𝑘 for 𝑘 = 1, 2, . . . , 𝑁 onde, as partículas possuem posições R𝑖 (R1 , R2 , . . . , R𝑁 ), momentos p𝑖 (p1 , p2 , . . . , p𝑁 ) e massas 𝑚𝑘 para 𝑘 = 1, 2, . . . , 𝑁 , os quais definem os termos de energia cinética (𝐾(p𝑖 )) e potencial de interação das partículas (𝑉 (R𝑖 )). 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 ) 𝑒−𝛽𝐻(R𝑖 ,p𝑖 ) 𝑑R𝑖 𝑑p𝑖 1 ⟨𝐴(R𝑖 , p𝑖 )⟩𝑇 = 1 2 ∫︁ ∫︁ 𝑁 ∫︁ , onde 𝛽 = . (2) 𝑘𝐵 𝑇 ... 𝑒−𝛽𝐻(R𝑖 ,p𝑖 ) 𝑑R𝑖 𝑑p𝑖 1
2
𝑁
Utilizando o Hamiltoniano definido na equação 1, em que 𝐾(p𝑖 ) e 𝑉 (R𝑖 ) são respectivamente dependentes apenas dos momentos e das posições, a energia média do sistema ⟨𝐸(R𝑖 , p𝑖 )⟩ = ⟨𝐴(R𝑖 , p𝑖 )⟩𝑇 pode ser calculada a partir das contribuições individuais de 𝐾(p𝑖 ) e 𝑉 (R𝑖 ) utilizando a equação :
⟨𝐻(R𝑖 , p𝑖 )⟩𝑇 = ⟨𝐾(p𝑖 )⟩𝑇 + ⟨𝑉 (R𝑖 )⟩𝑇 ∫︁ ∫︁ ∫︁ ∫︁ ∫︁ ∫︁ −𝛽𝐾(p𝑖 ) ... 𝐾(p𝑖 ) 𝑒 𝑑p𝑖 ... 𝑉 (R𝑖 ) 𝑒−𝛽𝑉 (R𝑖 ) 𝑑R𝑖 + 1 ∫︁2 ∫︁ 𝑁 ∫︁ .(3) = 1 ∫︁2 ∫︁ 𝑁 ∫︁ −𝛽𝐾(p𝑖 ) −𝛽𝑉 (R𝑖 ) ... 𝑒 𝑑p𝑖 ... 𝑒 𝑑R𝑖 1
2
𝑁
1
4
2
𝑁
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 (𝒵 ) pode ser obtida de forma exata e, assim, levando ao cálculo de qualquer propriedade termodinâmica. Logo, uma vez que ∫︁ ∫︁ ∫︁ ∫︁ ∫︁ ∫︁ ... 𝐻(R𝑖 ) 𝑒−𝛽𝑉 (R𝑖 ) 𝑑R𝑖 1 ∫︁2 ∫︁ 𝑁 ∫︁ ⟨𝐻(R𝑖 )⟩𝑇 = . ⇒ 𝒵= ... 𝑒−𝛽𝐻(R𝑖 ) 𝑑R𝑖 , (4) 1 2 𝑁 ... 𝑒−𝛽𝑉 (R𝑖 ) 𝑑R𝑖 1
2
𝑁
temos que 𝒵 é 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 𝒵 ) em um sistema de 𝑁 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 mecanismo 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 significa 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 configuraçõ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:
∫︁ ∫︁ 𝒵=
∫︁ ...
1
2
𝑒−𝛽𝑉 (R𝑖 ) 𝑑R𝑖 ≈
𝑁
𝜈 ∑︁
𝑒−𝛽𝑉𝜈 (R𝑖 ) .
(5)
𝜈=1
A partir dessa aproximação, é possível determinar a probabilidade de observação para cada configuração por 𝑒−𝛽𝑉𝜈 (R𝑖 ) . (6) 𝑃 (R𝜈𝑖 ) = 𝜈 ∑︁ 𝑒−𝛽𝑉𝜈 (R𝑖 ) 𝜈=1
Portanto, cada configuração 𝜈 é selecionada pelo algoritmo Metrópolis utilizando a amostragem 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 A𝑖𝑛 inicial para o estado B𝑠𝑒 seguinte:
A𝑖𝑛 ≡ {R𝜈𝑖 }
Prob. de Transição
−−−−−−−−−−−−−−→ 𝜋(A𝑖𝑛 →B𝑠𝑒 )
B𝑠𝑒 ≡ {R𝜈′ 𝑖 },
(7)
onde, a mudança de estado ocorre a partir da mudança da configuração espacial das partículas de {R𝜈𝑖 } para {R𝜈′ 𝑖 } com probabilidade dada por 𝜋(A𝑖𝑛 → B𝑠𝑒 ). A probabilidade de transição pode ser decomposta em {︃ 𝜋(A𝑖𝑛 → B𝑠𝑒 ) = 𝜉(A𝑖𝑛 → B𝑠𝑒 ) × acc(A𝑖𝑛 → B𝑠𝑒 )
𝜉(A𝑖𝑛 → B𝑠𝑒 ) ≡ Prob. do processo ocorrer acc(A𝑖𝑛 → B𝑠𝑒 ) ≡ 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: A𝑠𝑒 ≡ {R𝜈𝑖 }
Prob. de Transição
←−−−−−−−−−−−−−−−− 𝜋(B𝑖𝑛 →A𝑠𝑒 )
B𝑖𝑛 ≡ {R𝜈′ },
(8)
e assim, {︃ 𝜋(B𝑖𝑛 → A𝑠𝑒 ) = 𝜉(B𝑖𝑛 → A𝑠𝑒 ) × acc(B𝑖𝑛 → A𝑠𝑒 )
𝜉(B𝑖𝑛 → A𝑠𝑒 ) ≡ Prob. do processo ocorrer acc(B𝑖𝑛 → A𝑠𝑒 ) ≡ Prob. de aceitação
Tomando as probabilidades dos estados A𝑖𝑛 e B𝑠𝑒 definidas na distribuição de Boltzmann, temos 𝑒−𝛽𝑉𝜈 (R𝑖 ) 𝑒−𝛽𝑉𝜈′ (R𝑖 ) 𝑃A (R𝜈𝑖 ) = 𝜈 , 𝑃B (R𝜈′ . (9) )= 𝜈 𝑖 ∑︁ ∑︁ −𝛽𝑉𝜈 (R𝑖 ) −𝛽𝑉𝜈 (R𝑖 ) 𝑒 𝑒 𝜈=1
𝜈=1
Portanto, a aplicação do produto das probabilidades de encontrar os respectivos estados e suas transições, pode ser dada por (10)
𝑃A 𝜋(A𝑖𝑛 → B𝑠𝑒 ) = 𝑃B 𝜋(B𝑖𝑛 → A𝑠𝑒 ) 𝑃B 𝜋(A𝑖𝑛 → B𝑠𝑒 ) = 𝜋(B𝑖𝑛 → A𝑠𝑒 ) 𝑃A ⎤−1
⎡
⎢ −𝛽𝑉𝜈 (R ) ⎥ 𝑖 ⎢ 𝑒 ⎥ 𝜋(A𝑖𝑛 → B𝑠𝑒 ) 𝑒−𝛽𝑉𝜈′ (R𝑖 ) ⎢ 𝜈 ⎥ = 𝜈 × ∑︁ ⎢ ∑︁ ⎥ 𝜋(B𝑖𝑛 → A𝑠𝑒 ) ⎣ 𝜈 (R𝑖 ) 𝜈 (R𝑖 ) ⎦ 𝑒−𝛽𝑉 𝑒−𝛽𝑉 𝜈=1
𝜈=1
𝜋(A𝑖𝑛 → B𝑠𝑒 ) = 𝑒−𝛽𝐻𝜈′ (R𝑖 ) × 𝑒𝛽𝑉𝜈 (R𝑖 ) = 𝑒−𝛽[𝑉𝜈′ (R𝑖 )−𝑉𝜈 (R𝑖 )] 𝜋(B𝑖𝑛 → A𝑠𝑒 ) 𝜋(A𝑖𝑛 → B𝑠𝑒 ) = 𝑒−𝛽[𝑉𝜈′ (R𝑖 )−𝑉𝜈 (R𝑖 )] . 𝜋(B𝑖𝑛 → A𝑠𝑒 ) 6
(11)
Para o lado esquerdo, aplicando a decomposição na probabilidade de transição: ( (( 𝜉(A →(B(𝑠𝑒 ) (𝑖𝑛 ( ( (( 𝜉(B →(A(𝑠𝑒 ) (𝑖𝑛 (
×
acc(A𝑖𝑛 → B𝑠𝑒 ) = 𝑒−𝛽[𝑉𝜈′ (R𝑖 )−𝑉𝜈 (R𝑖 )] 𝜉(B𝑖𝑛 → A𝑠𝑒 ) acc(A𝑖𝑛 → B𝑠𝑒 ) = 𝑒−𝛽[𝑉𝜈′ (R𝑖 )−𝑉𝜈 (R𝑖 )] 𝜉(B𝑖𝑛 → A𝑠𝑒 )
(12)
A partir da equação 12 se aplica o critério de Metrópolis, de modo que: ⎧ 𝑃B ⎪ se 𝑃B < 𝑃A ⎨ 𝑃A , acc(A𝑖𝑛 → B𝑠𝑒 ) ⇒ ⎪ ⎩ 1, se 𝑃B ≥ 𝑃A A partir dos critérios de Metrópolis definidos, o protocolo seguido pelo algoritmo é disposto da seguinte forma: Gerar uma configuração aleatória 𝑐2 a partir de uma configuração prévia 𝑐1 . Se a energia de 𝑐2 for menor do que 𝑐1 , 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 𝑃𝑃AB , aceita-se a configuração 𝑐2 que será utilizada como configuração 𝑐1 para um novo movimento. Caso contrário, a configuração 𝑐1 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 ⎧ ⃗ 1 𝑑Z ⃗1 ⃗ 1 𝑑Y ⎪ 𝑑X ⎪ ⎪ ⃗v1 = + + ⎪ ⎪ 𝑑𝑡 𝑑𝑡 𝑑𝑡 ⎪ ⎧ ⎪ ⎪ ⃗ ⃗ ⃗2 ⎪ 𝑑 X 𝑑 Y 𝑑 Z ⃗ ⃗ ⃗ ⃗ ⎪ ⎪ 2 2 R1 = X1 + Y1 + Z1 ⎪ ⎪ ⃗ v = + + ⎪ ⎪ 2 ⎪ ⎪ ⎪ ⎪ 𝑑𝑡 𝑑𝑡 𝑑𝑡 ⃗2=X ⃗2+Y ⃗2+Z ⃗2 ⎪ ⎪ R ⎪ ⎪ ⎪ ⎪ ⃗ ⃗ ⃗3 ⎪ ⎪ 𝑑 Y 𝑑 Z 𝑑 X 3 3 ⎨R ⎨ ⃗v = ⃗3=X ⃗3+Y ⃗3+Z ⃗3 + + 3 ⃗ 𝑁} = 𝑑𝑡 𝑑𝑡 𝑑𝑡 {R {⃗v𝑖𝑁 } = 𝑖 ⃗4=X ⃗4+Y ⃗4+Z ⃗4 R ⎪ ⎪ ⃗ ⃗ ⃗4 ⎪ ⎪ 𝑑 X4 𝑑 Y 4 𝑑 Z ⎪ ⎪ ⎪ ⎪ . ⎪ ⎪ ⃗ v = + + . 4 ⎪ ⎪ . ⎪ ⎪ 𝑑𝑡 𝑑𝑡 𝑑𝑡 ⎪ ⎪ ⎪ ⎪ ⎩R ⎪ .. ⃗𝑁 =X ⃗𝑁 +Y ⃗𝑁 +Z ⃗𝑁 ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ ⃗𝑁 ⃗𝑁 ⃗𝑁 ⎪ 𝑑X 𝑑Y 𝑑Z ⎪ ⎩ ⃗v𝑁 = + + 𝑑𝑡 𝑑𝑡 𝑑𝑡 7
⃗ 𝑁 } e {⃗v𝑁 } sendo o conjunto de posições e velocidades, respectivamente, das com {R 𝑖 𝑖 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 𝑁 ∑︁
𝐹⃗𝑖 = −∇R ⃗𝑖
𝑖=1
𝑁 ∑︁
⃗ 𝑁 }) = −∇ ⃗ 𝑉 ({R 𝑖 R𝑖
𝑖=1
𝑁 ∑︁
⃗𝑖−R ⃗ 𝑗 |), 𝑉 (|R
(13)
𝑗=1
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: 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁
𝐹⃗𝑖 = 𝐹⃗𝑖 =
𝑖=1
−∇R ⃗𝑖
𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑖=1
𝑁 ∑︁
(14)
𝑚𝑖⃗a𝑖 𝑚𝑖
⃗ 𝑁} ⃗1 ⃗2 ⃗3 ⃗𝑁 𝑑2 { R 𝑑2 𝑅 𝑑2 𝑅 𝑑2 𝑅 𝑑2 𝑅 𝑖 = 𝑚 + 𝑚 + 𝑚 + . . . + 𝑚 1 2 3 𝑁 𝑑𝑡2 𝑑𝑡2 𝑑𝑡2 𝑑𝑡2 𝑑𝑡2
⃗𝑖−R ⃗ 𝑗 |) = 𝑚1 𝑉 (|R
𝑗=1
⃗1 ⃗2 ⃗3 ⃗𝑁 𝑑2 𝑅 𝑑2 𝑅 𝑑2 𝑅 𝑑2 𝑅 + 𝑚 + 𝑚 + . . . + 𝑚 2 3 𝑁 𝑑𝑡2 𝑑𝑡2 𝑑𝑡2 𝑑𝑡2
Logo,
−∇R ⃗𝑖
𝑁 ∑︁
⃗𝑖−R ⃗ 𝑗 |) = 𝑉 (|R
𝑗=1
𝑁 ∑︁ 𝑖=1
𝑚𝑖
⃗ 𝑁} 𝑑2 { R 𝑖 𝑑𝑡2
(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 ⃗ 𝑁} 𝑑2 {R
𝑖 diversos algoritmos foram propostos com base no Para integração do termo 𝑑𝑡2 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 𝑔(𝑥), sua expansão em intervalos ∆𝑥 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, estiramento 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:
𝑔(𝑥 + ∆𝑥) =
𝑁 ∑︁ 𝑔 (𝑛)
𝑛!
𝑛=0
= 𝑔+
(𝑥 − ∆𝑥)𝑛
𝑑𝑔 (𝑥 − ∆𝑥)1 𝑑2 𝑔 (𝑥 − ∆𝑥)2 𝑑3 𝑔 (𝑥 − ∆𝑥)3 + 2 + 3 + ... 𝑑𝑥 1! 𝑑𝑥 2! 𝑑𝑥 3!
As difereças finitas podem ser obtidas da expansão e retrocesso em intervalos de ℎ e −ℎ, respectivamente, ou seja:
𝑑𝑔 ℎ+ 𝑑𝑥 𝑑𝑔 𝑔(𝑥 − ℎ) = 𝑔 − ℎ + 𝑑𝑥 𝑔(𝑥 + ℎ) = 𝑔 +
𝑑2 𝑔 ℎ2 𝑑3 𝑔 ℎ3 + 3 + 𝜗(ℎ4 ) + . . . + 𝜗(ℎ𝑛 ) 𝑑𝑥2 2 𝑑𝑥 6 𝑑2 𝑔 ℎ2 𝑑3 𝑔 ℎ3 − 3 + 𝜗(ℎ4 ) + . . . + 𝜗(ℎ𝑛 ) 2 𝑑𝑥 2 𝑑𝑥 6
onde, 𝜗(ℎ𝑛 ) são os termos de ordens superiores. Desse modo, considerando o incremento ⃗ 𝑖 na equação 15, ℎ = 𝑡 + ∆𝑡, e aplicando o método de diferenças finitas na expansão de R 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⃗ 𝑖 da partícula para os instantes 𝑡+∆𝑡 e 𝑡−∆𝑡 nas séries de Taylor, truncando se a posição R a série a partir do quarto termo, ou seja, 2⃗ 2 2 ⃗ ⃗ 𝑖 (𝑡 + Δ𝑡) = R ⃗ 𝑖 (𝑡) + ⃗v𝑖 (𝑡)Δ𝑡 + ⃗a𝑖 (𝑡) Δ𝑡 + . . . ⃗ 𝑖 (𝑡 + Δ𝑡) = R ⃗ 𝑖 (𝑡) + 𝑑R𝑖 (𝑡) Δ𝑡 + 𝑑 R𝑖 (𝑡) Δ𝑡 + . . . , R R 𝑑𝑡 𝑑𝑡2 2 2 2 2⃗ 2 ⃗ Δ𝑡 𝑑 R (𝑡) 𝑑 R (𝑡) Δ𝑡 𝑖 𝑖 ⃗ 𝑖 (𝑡 − Δ𝑡) = R ⃗ 𝑖 (𝑡) − ⃗v𝑖 (𝑡)Δ𝑡 + ⃗a𝑖 (𝑡) ⃗ 𝑖 (𝑡 − Δ𝑡) = R ⃗ 𝑖 (𝑡) − − ... R R Δ𝑡 + − . . . , 2 𝑑𝑡 𝑑𝑡2 2
⃗ 𝑖 (𝑡 + ∆𝑡) e R ⃗ 𝑖 (𝑡 − ∆𝑡), temos: Somando as equações R ⃗ 𝑖 (𝑡 + ∆𝑡) + R ⃗ 𝑖 (𝑡 − ∆𝑡) = 2R ⃗ 𝑖 (𝑡) + ⃗a𝑖 (𝑡)∆𝑡2 + . . . R ⇓ ⃗ ⃗ ⃗ 𝑖 (𝑡 − ∆𝑡) + ⃗a𝑖 (𝑡)∆𝑡2 + . . . , R𝑖 (𝑡 + ∆𝑡) = 2R𝑖 (𝑡) − R
(16)
onde, uma vez que,
𝐹𝑖 ⃗a𝑖 (𝑡) = , logo, 𝑚𝑖
⃗ 𝑖 (𝑡 + ∆𝑡) = 2R ⃗ 𝑖 (𝑡) − R ⃗ 𝑖 (𝑡 − ∆𝑡) + R
(︂
𝐹𝑖 𝑚𝑖
)︂
∆𝑡2 + . . . ,
(17)
e assim,
[︃
]︃ 𝑁 ∑︁ 1 ⃗ 𝑖 (𝑡 + ∆𝑡) = 2R ⃗ 𝑖 (𝑡) − R ⃗ 𝑖 (𝑡 − ∆𝑡) − ⃗𝑖−R ⃗ 𝑗 |) ∆𝑡2 . R ∇⃗ 𝑉 (|R 𝑚𝑖 R𝑖 𝑗=1
(18)
A posição da partícula 𝑖 em 𝑡+∆𝑡 depende da posição no tempo anterior a 𝑡 e na diferencial do potencial de interação com relação à posição previamente calculado. 9
⃗ 𝑖 (𝑡 + ∆𝑡) - R ⃗ 𝑖 (𝑡 − ∆𝑡), onde, Para as velocidades, tomamos a diferença em R ⃗ 𝑖 (𝑡 + ∆𝑡) − R ⃗ 𝑖 (𝑡 − ∆𝑡) = 2⃗v𝑖 (𝑡)∆𝑡 R ⇓ ⃗ 𝑖 (𝑡 + ∆𝑡) − R ⃗ 𝑖 (𝑡 − ∆𝑡) R ⃗v𝑖 (𝑡) = 2∆𝑡
(19)
⃗ 𝑖 (𝑡 + ∆𝑡) e ⃗v𝑖 (𝑡) que Observamos a partir das operações realizadas na determinação de R o erro das posições com o algoritmo de Verlet é da ordem de (∆𝑡)4 , e das velocidades de (∆𝑡)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 termodinâ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 algum 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 barreiras 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 posteriormente 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 consiste no uso de 𝑀 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 𝑀 temperaturas 𝑇𝑚 (𝑚 = 1, 2, . . . , 𝑀 ) para 𝑀 réplicas, ou seja, como temperaturas paralelas (parallel tempering ). De modo simbólico: ⎧ ⎧ ⎪ ⎪ Coordenadas relativas R𝑖 , p𝑖 : Estados e temperaturas : ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨x1 (R, p) ⇒ 𝑇1 ⎨X1 (x1 ) ⇒ 𝑇1 Generalized ensemble x2 (R, p) ⇒ 𝑇2 ⇒ X2 (x2 ) ⇒ 𝑇2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ . . . ... ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩X (x ) ⇒ 𝑇 ⎩x (R, p) ⇒ 𝑇 𝑖 𝑚 𝑖 𝑖 𝑚 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 𝑖(𝑖 = 1, 2, . . . , 𝑀 ) coorresponde à identificação para as réplicas e 𝑚(𝑚 = 1, 2, . . . , 𝑀 ) para as temperaturas. Logo, podemos representar os estados no generalized ensemble por [𝑖(1)]
X = (x1 ⏟
[𝑖(2)]
, x2
[𝑖(𝑀 )]
, . . . , x𝑀 ⏞
Replica exchange
[1]
[2]
[𝑀 ]
) = (x𝑚(1) , x𝑚(2) , . . . , x𝑚(𝑀 ) ), ⏟ ⏞ Temperature exchange
10
(20)
onde, o estado X é dado por 𝑀 réplicas, de modo que em cada uma delas as coordenadas no espaço de fase das 𝑁 partículas são dadas por R[𝑖] e p[𝑖] em uma das temperaturas em 𝑇𝑚 , identificadas pelo índice [𝑖] como notação para as réplicas. Definiremos agora a troca de um par de réplicas no generalized ensemble identificadas por 𝑖 e 𝑗 , imersas em termostatos de temperaturas 𝑇𝑚 e 𝑇𝑛 , respectivamente, o que resulta na troca de dois estados como [𝑖]′ ′ [𝑗]′ [𝑗] X = (. . . , x[𝑖] 𝑚 , . . . , x𝑛 , . . .) → X = (. . . , x𝑚 , . . . , x𝑛 , . . .).
(21)
Logo, podemos escrever as trocas sob a perspectiva das réplicas e das temperaturas da seguinte forma: {︃ [𝑖] Swap [𝑗]′ x𝑚 ≡ (R[𝑖] , p[𝑖] )𝑚 −−−−−−−→ x𝑚 ≡ (R[𝑗]′ , p[𝑗]′ )𝑚 Replicas ′ Swap [𝑗] [𝑖]′ x𝑛 ≡ (R[𝑗] , p[𝑗] )𝑛 −−−−−−−→ x𝑛 ≡ (R[𝑖]′ , p[𝑖]′ )𝑛 Uma vez que não há interação entre as 𝑀 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 {︃ 𝑀 }︃ 𝑀 ∏︁ ∑︁ [𝑖] [𝑖] 𝑊 (X) = 𝑒−𝛽𝑚 𝐻(R ,p ) = exp −𝛽𝑚 𝐻(R[𝑖] , p[𝑖] ) . (22) 𝑖=1
𝑖=1
Omitindo os argumentos R[𝑖] , p[𝑖] de 𝐻(R[𝑖] , p[𝑖] ), consideraremos a troca das réplicas expressas da seguinte forma: ⎧ ⎫ ⎨ ⎬ [1(1)] [𝑖(𝑚)] [𝑗(𝑛)] [𝑀 (𝑀 )] exp −𝛽1 𝐻 + . . . − 𝛽𝑚 𝐻 − 𝛽𝑛 𝐻 + . . . − 𝛽𝑛 𝐻 ⏟ ⏞ ⎩ ⎭ troca
.. .
⎧ ⎫ ⎨ ⎬ exp −𝛽1 𝐻 [1(1)] + . . . − 𝛽𝑚 𝐻 [𝑗(𝑛)] − 𝛽𝑛 𝐻 [𝑖(𝑚)] + . . . − 𝛽𝑛 𝐻 [𝑀 (𝑀 )] ⏞ ⏟ ⎩ ⎭ Réplicas trocadas
Desse modo, temos a troca de temperatura de duas réplicas 𝑖 e 𝑗 pela transição de estado dada por
𝑊 (X)
⇒
𝑊 (X′ ),
Com probabilidade de troca
𝑤(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,
𝑊 (X)𝑤(X → X′ ) = 𝑊 (X′ )𝑤(X′ → X) 𝑤(X → X′ ) 𝑊 (X′ ) = , 𝑤(X′ → X) 𝑊 (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 Hukushima, 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 analí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ículas, de modo que se aplica a probabilidade 𝑃𝛽𝑚 (R) (sistema com conjunto 𝑘 de partículas) para ocorrência do estado X em uma seleção aleatória,ou seja,
𝑊 (X) =
𝑀 ∏︁
[𝑖] )
𝑒−𝛽𝑉𝑘 (R) 𝑃𝛽𝑚 (R) = ∑︁ 𝑒−𝛽𝑉𝑘 (R)
onde
𝑒−𝛽𝑚 𝑉 (R
𝑖=1
(24)
𝑘
A partir dessas definições, a distribuição de probabilidades das 𝑀 réplicas (Ω𝛽 ) é dada a partir das temperaturas elencadas em 𝑚(𝑚 = 1, 2, . . . , 𝑀 ) pelo produtório das probabilidades de cada estado onde propomos a troca: 𝛽
Ω =
𝑀 ∏︁
[1]
[𝑗]
[𝑖]
[𝑀 ]
𝑃𝛽𝑚 (R) = [𝑃𝛽1 (R) × · · · × 𝑃𝛽𝑚 (R) × 𝑃𝛽𝑛 (R) × . . . × 𝑃𝛽𝑀 (R)] ⏟ ⏞ 𝑖=1
(25)
Troca
Portanto, no uso da condição de balanço detalhado, considerando as 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 𝑀 ∏︁
′
𝑃𝛽𝑚 (R)𝑤(X → X ) =
𝑖=1
[︁
... ×
[𝑖] 𝑃𝛽𝑚 (R)
×
[𝑗] 𝑃𝛽𝑛 (R)
𝑀 ∏︁
𝑃𝛽𝑛 (R)𝑤(X′ → X)
𝑖=1
]︁
[𝑗]
[𝑖]
× . . . 𝑤(X → X′ ) = [. . . × 𝑃𝛽𝑚 (R) × 𝑃𝛽𝑛 (R) × . . .]𝑤(X′ → X).
De modo similar ao método Monte Carlo, onde a probabilidade de transição configuracional é 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 probabilidade da troca ocorrer (k(X → X′ )) com a probabilidade de aceitação (ACC(X → X′ ). Portanto, pela condição de balanço detalhado temos: 𝑀 ∏︁
𝑃𝛽𝑚 (R)k(X → X′ )ACC(X → X′ ) =
𝑖=1
𝑀 ∏︁
𝑃𝛽𝑛 (R)k(X𝑛 → X′𝑚 )ACC(X𝑛 → X′𝑚 )
𝑖=1 𝑀 ∏︁
k(X → X′ )ACC(X → X′ ) = 𝑖=1 ′ ′ ′ 𝑀 k(X → X )ACC(X → X) ∏︁ ′ X k(X → )ACC(X → X′ ) = ′ ′ k(X X )ACC(X′ → X′ ) →
𝑃𝛽𝑛 (R) 𝑃𝛽𝑚 (R)
𝑖=1 [1] 𝑃 𝛽 1 (R) × . . . [1] × ... 𝑃𝛽 (R) 1 [𝑖]
[𝑗]
[𝑀 ]
[𝑖]
[𝑗]
[𝑀 ]
(R) 𝑃𝛽𝑚 (R) × 𝑃𝛽𝑛 (R) . . .× 𝑃 𝛽𝑀 [𝑗]
𝑃𝛽𝑛 (R) × 𝑃𝛽𝑚 (R) ACC(X𝑚 → X′𝑛 ) = , [𝑖] [𝑗] ACC(X𝑛 → X′𝑚 ) 𝑃𝛽𝑚 (R) × 𝑃𝛽𝑚 (R) 12
[𝑖]
(R) 𝑃𝛽𝑛 (R) × 𝑃𝛽𝑚 (R) . . .× 𝑃 𝛽𝑀
(26)
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: ⎡ ⎡ ⎤−1 ⎤−1 [𝑖] [𝑗] ACC(X → X′ ) 𝑒−𝛽𝑛 𝑉 (R) ⎢ 𝑒−𝛽𝑚 𝑉 (R) ⎥ ⎢ ⎥ = ∑︁ × ⎣ ∑︁ ⎦ −𝛽 𝑉 −𝛽 𝑉𝑘 (R) 𝑘 𝑘 𝑘 ACC(X′ → X) 𝑒 𝑒
𝑘
[𝑖] ⎢ 𝑒−𝛽𝑛 𝑉 [𝑗] (R) ⎥ 𝑒−𝛽𝑚 𝑉 (R) ⎢ ⎥ × ∑︁ × ⎣ ∑︁ −𝛽 ⎦ −𝛽 𝑉 (R) 𝑘 𝑘 𝑉 (R) 𝑒 𝑒
𝑘
−𝛽𝑛 𝑉 [𝑖] (R)
= ⏟𝑒
𝑘
−𝛽𝑚 𝑉 [𝑗] (R)
×𝑒
⏞× 𝑒
𝛽𝑚 𝑉 [𝑖] (R)
×𝑒
𝑘 𝛽𝑛 𝑉 [𝑗] (R)
Agrupando
[︀ ]︀ = exp −𝛽𝑛 𝑉 [𝑖] (R) − 𝛽𝑚 𝑉 [𝑗] (R) + 𝛽𝑚 𝑉 [𝑖] (R) + 𝛽𝑛 𝑉 [𝑗] (R) {︀ [︀ ]︀ [︀ }︀ = exp −𝛽𝑛 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) + 𝛽𝑚 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) {︀ [︀ ]︀}︀ = exp (𝛽𝑚 − 𝛽𝑛 ) 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) [︀ ]︀ Considerando (𝛽𝑚 − 𝛽𝑛 ) 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) = ∆, temos ACC(X → X′ ) = 𝑒Δ ACC(X′ → X)
(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 momentos lineares das partículas nas réplicas trocadas devems er considerados, sendo, então, ajustados pela seguinte relação: √︂ √︂ 𝑇𝑛 [𝑖] 𝑇𝑚 [𝑗] [𝑗]′ [𝑖]′ 𝑝 , 𝑝 = 𝑝 . (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:
𝑊 (X) =
𝑀 ∏︁
[𝑖]
−𝛽𝑚 𝐻(R[𝑖] ,p[𝑖] )
𝑒
onde
𝑖=1
𝑒−𝛽𝐻𝑘 (R,p ) 𝑃𝛽𝑚 (R, p) = ∑︁ . [𝑖] 𝑒−𝛽𝐻𝑘 (R,p ) 𝑘
13
(29)
O protocolo de desenvolvimento com base na troca de duas réplicas segue o mesmo formalismo 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)
[𝑖]
=
[𝑗]
𝑃𝛽𝑛 (R, p) × 𝑃𝛽𝑚 (R, p) [𝑖]
[𝑗]
𝑃𝛽𝑚 (R, p) × 𝑃𝛽𝑚 (R, p) ⎡ [𝑖]
=
⎤−1 [𝑗]
[𝑖]
[𝑗]′
⎢ 𝑒−𝛽𝑚 𝐻(R ,p ) (R) ⎥ 𝑒−𝛽𝑛 𝐻(R ,p ) ⎢ ∑︁ ∑︁ × ( ( ⎥ ( ( ⎣ ⎦ 𝐻𝑘( (R,p) 𝐻𝑘( (R,p) (𝑘( (𝑘( ( ( 𝑒−𝛽 𝑒−𝛽 𝑘
⎤−1
⎡ [𝑖]
[𝑖]
[𝑗]
[𝑗]
⎢ 𝑒−𝛽𝑛 𝐻(R ,p ) ⎥ 𝑒−𝛽𝑚 𝐻(R ,p ) ⎢ ∑︁ × ∑︁ × ( (⎥ ( ( ⎣ ⎦ 𝐻𝑘( (R,p) 𝐻𝑘( (R,p) (𝑘( (𝑘( ( ( 𝑒−𝛽 𝑒−𝛽
𝑘
𝑘
𝑘
[︁ ]︁ ′ ′ = exp −𝛽𝑛 𝐻(R[𝑖] , p[𝑖] ) − 𝛽𝑚 𝐻(R[𝑗] , p[𝑗] ) + 𝛽𝑚 𝐻(R[𝑖] , p[𝑖] ) + 𝛽𝑛 𝐻(R[𝑗] , p[𝑗] )
Uma vez que 𝐻(R[𝑖] , p[𝑖] ) = 𝐾(p[𝑖] ) + 𝑉 (R[𝑖] ),
{︁ [︁ ]︁ [︁ ]︁}︁ ACC(X → X′ ) [𝑖]′ [𝑖] [𝑗]′ [𝑗] = exp −𝛽 𝐾(p ) + 𝑉 (R ) − 𝛽 𝐾(p ) + 𝑉 (R ) 𝑛 𝑚 ACC(X′ → X) ]︀ [︀ ]︀}︀ {︀ [︀ + exp 𝛽𝑚 𝐾(p[𝑖] ) + 𝑉 (R[𝑖] ) + 𝛽𝑛 𝐾(p[𝑗] ) + 𝑉 (R[𝑗] ) . Organizando em blocos de 𝐾(p) e 𝑉 (R), temos: ⎧ [︁ ]︁ [𝑖]′ [𝑗]′ [𝑖] [𝑗] ⎪ ⎪ exp −𝛽𝑛 𝐾(p ) − 𝛽𝑚 𝐾(p ) + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p ) ACC(X → X′ ) ⎨ = + ACC(X′ → X) ⎪ ⎪ [︀ ]︀ ⎩ [𝑖] exp −𝛽𝑛 𝑉 (R ) − 𝛽𝑚 𝑉 (R[𝑗] ) + 𝛽𝑚 𝑉 (R[𝑖] ) + 𝛽𝑛 𝑉 (R[𝑗] ) . 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 𝛽𝑚 = 1/𝑘𝐵 𝑇𝑚 , temos: [︁ ]︁ [𝑖]′ [𝑗]′ [𝑖] [𝑗] exp −𝛽𝑛 𝐾(p ) − 𝛽𝑚 𝐾(p ) + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p )
exp exp exp exp
]︂ [︂ ′ ′ (p[𝑗] )2 (p[𝑖] )2 [𝑖] [𝑗] − 𝛽𝑚 + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p ) −𝛽𝑛 2𝑚 2𝑚 (︂ )︂ [𝑗] 2 ]︂ [︂ (︂ )︂ [𝑖] 2 𝑇𝑛 (p ) 𝑇𝑚 (p ) [𝑖] [𝑗] −𝛽𝑛 − 𝛽𝑚 + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p ) 𝑇𝑚 2𝑚 𝑇𝑛 2𝑚 [︂ (︂ )︂ [𝑖] 2 (︂ )︂ [𝑗] 2 ]︂ 1 𝑇 (p ) 1 𝑇 (p ) 𝑛 𝑚 [𝑖] [𝑗] − − + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p ) 𝑘𝐵 𝑇 𝑇𝑚 2𝑚 𝑘𝐵 𝑇 𝑇𝑛 2𝑚 𝑛 𝑚 [︂ ]︂ 1 (p[𝑖] )2 1 (p[𝑗] )2 [𝑖] [𝑗] − − + 𝛽𝑚 𝐾(p ) + 𝛽𝑛 𝐾(p ) 𝑘𝐵 𝑇𝑚 2𝑚 𝑘𝐵 𝑇𝑛 2𝑚 ′
′
Observemos a conversão de 𝐾(p[𝑖] ) → 𝐾(p[𝑖] ) e 𝐾(p[𝑗] ) → 𝐾(p[𝑗] ), de modo que: [︁ ]︁ [𝑖] [𝑗] [𝑗] [𝑖] exp −𝛽 𝛽𝑛 𝐾(p ) + 𝛽𝑚𝐾(p ) + 𝛽𝑛 𝐾(p ) , 𝑚 𝐾(p ) − 14
e assim, o primeiro termo se anula [︁ ]︁ ′ ′ exp −𝛽𝑛 𝐾(p[𝑖] ) − 𝛽𝑚 𝐾(p[𝑗] ) + 𝛽𝑚 𝐾(p[𝑖] ) + 𝛽𝑛 𝐾(p[𝑗] ) = 0,
(30)
de modo que:
[︀ ]︀ ACC(X → X′ ) = exp −𝛽𝑛 𝑉 (R[𝑖] ) − 𝛽𝑚 𝑉 (R[𝑗] ) + 𝛽𝑚 𝑉 (R[𝑖] ) + 𝛽𝑛 𝑉 (R[𝑗] ) ′ ACC(X → X) {︀ [︀ ]︀ [︀ }︀ = exp −𝛽𝑛 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) + 𝛽𝑚 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) {︀ [︀ ]︀}︀ = exp (𝛽𝑚 − 𝛽𝑛 ) 𝑉 [𝑖] (R) − 𝑉 [𝑗] (R) ACC(X → X′ ) = 𝑒Δ (31) ACC(X′ → X) 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: ⎧ Δ ⎪ se ∆>0 ⎨𝑒 , ′ ACC(X → X ) ⇒ (32) ⎪ ⎩ 1, se ∆≤0 À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) é o algoritmo entre os 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 Algorithms 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 properties 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
