Termodinâmica estatística de líquidos com o método de Monte Carlo. I. Metodologia

Statistical thermodynamics of liquids using the Monte Carlo method. I. Methodology

Monte Carlo simulation; free energy simulation; electrostactic potential; CHELPG charges

DIVULGAÇÃO

Termodinâmica estatística de líquidos com o método de Monte Carlo. I. Metodologia

Vania Elisabeth Barlette

Instituto de Química - Universidade Estadual de Campinas - CP 6154 - 13083-970 - Campinas - SP

Luiz Carlos Gomide Freitas

Departamento de Química - Universidade Federal de São Carlos - CP 676 - 13565-905 - São Carlos - SP

Recebido em 18/2/97; aceito em 14/4/98

Statistical thermodynamics of liquids using the Monte Carlo method. I. Methodology. Statistical mechanics Monte Carlo simulation is reviewed as a formalism to study thermodynamic properties of liquids. Considering the importance of free energy changes in chemical processes, the thermodynamic perturbation theory implemented in the Monte Carlo method is discussed. The representation of molecular interaction by the Lennard-Jones and Coulomb potential functions is also discussed. Charges derived from quantum molecular electrostatic potential are also discussed as an useful methodology to generate an adequate set of partial charges to be used in liquid simulation.

Keywords: Monte Carlo simulation; free energy simulation; electrostactic potential; CHELPG charges.

1. INTRODUÇÃO

Os líquidos estão presentes em uma quantidade inumerável de sistemas na natureza e desempenham papel de grande importância em ciência básica, tecnologia e em processos vitais. Avanços no conhecimento do estado líquido tem sido obtidos utilizando teorias analíticas a partir da formulação estatística de Boltzmann e Gibbs para estados microscópicos da matéria^1-4. Entretanto, devido a dificuldades no tratamento matemático das equações envolvidas, esses estudos restringem-se a abordagem de líquidos formados por moléculas simples. Em décadas recentes, o desenvolvimento de novas técnicas computacionais e experimentais permitiu uma nova abordagem qualitativa e quantitativa para o estudo de líquidos^5-8. Historicamente, podemos relacionar o início dessa nova etapa com a publicação dos trabalhos pioneiros de Metropolis e col.⁹ com o método de Monte Carlo (1953), e Alder e Wainwright^10-11 com o método de Dinâmica Molecular (1957,1959). Essas técnicas computacionais permitem a investigação de líquidos moleculares complexos.

A utilização de métodos computacionais específicos para estudar sistemas líquidos está associada à complexidade das interações intermoleculares e à ausência de simetria espacial que caracterizam estes sistemas. Com isso, o desenvolvimento de teorias analíticas está restrito a modelos moleculares simplificados. Estudos analíticos pressupõem o auxílio de aproximações matemáticas para a solução de equações e, em geral, os resultados obtidos nem sempre refletem a qualidade do modelo proposto para a representação molecular dos componentes do sistema estudado. Contrastando a essas dificuldades, as técnicas de simulação são capazes de tratar exatamente o modelo e testar sua qualidade a partir da comparação com resultados experimentais.

Os modelos moleculares para o estudo de líquidos com simulação molecular podem levar em consideração características e propriedades das moléculas que constituem o sistema, e o modo como estas interagem. A metodologia de mecânica estatística é, então, utilizada para reproduzir as relações entre interações intermoleculares e propriedades termodinâmicas^12,13. Pode-se afirmar que os potenciais intermoleculares empregados são cruciais para a adequada descrição das relações entre propriedades moleculares e características observadas experimentalmente.

Neste trabalho, apresentamos o formalismo básico do método de Monte Carlo, discutindo sua implementação para calcular propriedades termodinâmicas de líquidos. A utilização dessa metodologia para desenvolver um modelo molecular para o clorofórmio e estudar propriedades termodinâmicas deste líquido será apresentada em trabalho subsequente.

2. SIMULAÇÕES DE MONTE CARLO

2.1. Considerações Gerais

Simulações de Monte Carlo (MC), bem como simulações de Dinâmica Molecular (MD), utilizam sistemas protótipos, os quais interagem através de um campo de força conhecido, e fornecem valiosas informações sobre termodinâmica, estrutura e, no caso de MD, dinâmica de líquidos. Esses métodos modelam os sistemas a nível molecular. Interações intermoleculares específicas entre soluto-solvente e solvente-solvente, necessárias para fornecer detalhes sobre a estrutura e a origem dos efeitos de solvente, também podem ser tratadas. O sucesso desses métodos reside, principalmente, na qualidade dos modelos de potencial empregados, sem comprometer a viabilidade computacional dos cálculos. A aplicação desses modelos a sistemas de tamanho reduzido, comparados aos sistemas macroscópicos, são suficientes para permitir medidas de propriedades observáveis em muito bom acordo com dados experimentais.

As relações entre propriedades termodinâmicas e níveis de energia mecânico-quânticos de uma molécula são fornecidas através da função partição^4,12,13. Métodos de mecânica quântica ab initio e semiempíricos permitem a determinação da função partição para moléculas constituídas por poucos átomos, através da determinação dos seus níveis de energia. Para sistemas compostos por muitas moléculas, no entanto, um tratamento clássico é adequado. Nessa aproximação, a energia do sistema é representada por uma Hamiltoniana , que depende somente das posições e dos momenta das partículas constituintes. A função partição é, então, transformada de uma soma sobre os níveis de energia para uma integral multidimensional sobre todas as possíveis posições e momenta de todas as partículas. O cálculo de funções termodinâmicas, como entalpia, entropia e demais propriedades relacionadas às suas flutuações, as quais dependem diretamente da função partição, requer o conhecimento de todas as possíveis configurações espaciais do sistema, e a resolução de integrais multidimensionais. Entretanto, a despeito da redefinição do sistema em termos clássicos, ainda assim, a função partição não pode ser determinada quando potenciais realísticos são utilizados para representar a interação intermolecular. Para sistemas complexos e de alta entropia, como sistemas moleculares em fase líquida, muitas configurações espaciais são necessárias para a sua adequada caracterização. O algoritmo de Metropolis implementado no método de MC consiste, basicamente, de um processo de evolução probabilística com a finalidade de calcular valores médios de observáveis mecânicas sem a necessidade de determinar a função partição. Esse método foi proposto por Metropolis e colaboradores⁹, e é muito eficiente para explorar o comportamento termodinâmico de sistemas em fase líquida.

As propriedades, ou observáveis, obtidas via simulação de MC, aparecem como médias sobre algum espaço amostral. A fim de obter as configurações espaciais que farão parte da amostragem, utiliza-se um gerador de números aleatórios para modificar as coordenadas espaciais das moléculas que compõem o sistema. é importante ressaltar que novas configurações são escolhidas de modo probabilístico, de acordo com um peso estatístico, proporcional à distribuição de equilíbrio termodinâmico do sistema, a distribuição de Boltzmann. Ou seja, ao processo de gerar configurações não está associada uma lei dinâmica de evolução temporal. A cinética de que trata o método de MC é uma cinética estocástica, em que configurações sucessivas, ou passos de MC, não são, necessariamente, temporalmente consecutivas.

Empregar a distribuição estatística de Boltzmann implica que as configurações que não são energeticamente muito próximas à energia média de equilíbrio termodinâmico tem uma contribuição desprezível para o cálculo dos valores médios. Isso requer um tipo especial de amostragem para acelerar a convergência do cálculo dos valores médios. Esse tipo de amostragem é denominado importance sampling. As médias são, então, calculadas sobre o espaço configuracional amostrado, constituido de um número suficientemente grande de configurações espaciais, e escolhidas de acordo com os pesos estatísticos de Boltzmann. Os valores médios calculados com essa metodologia são exatos para o modelo utilizado para representar a interação intermolecular.

2.2. Valores médios e ensembles

As variáveis para descrever um modelo usualmente são escolhidas de acordo com o sistema a ser estudado. Em geral, em simulação de líquidos as moléculas são consideradas rígidas. Nesse caso, a posição x_a de uma molécula não linear a pode ser especificada por um vetor 6-dimensional composto da posição do centro-de-massa da molécula, ou outro ponto conveniente, e os três ângulos de Euler,

x_a = (x_a,y_a,z _a,f_a,q _a,y_a) (1)

Então, uma configuração espacial inteira de um sistema de N moléculas é dada pelo vetor 6 N-dimensional,

x^N = (x₁,x₂,...,x _N). (2)

No formalismo da Termodinâmica Estatística,^4,12,13 uma propriedade termodinâmica observável é uma média no ensemble de estados microscópicos de um sistema. A escolha apropriada de um ensemble está associada às condições de contorno impostas ao sistema em estudo. Os sistemas comumente estudados em simulação de líquidos estão sob as condições do ensemble canônico, com temperatura T, volume V e número N de moléculas fixo (NVT), e ensemble isobárico-isotérmico (NpT), em que p é a pressão do sistema. Seguiremos, portanto, com a discussão para o ensemble NVT com a sua extensão para o ensemble N_pT.

A integral que define o valor médio de um observável pode ser dividida em duas partes: uma contendo a dependência com relação aos momenta das partículas e, a outra, contendo a dependência em relação às coordenadas espaciais. Independente da complexidade do sistema, a integral sobre os momenta das partículas pode ser resolvida analiticamente^{4-6, 12-14}. Devido a essa particularidade, a complexidade matemática está na resolução da integral que contém o potencial de interação intermolecular, a qual é denominada integral configuracional. A contribuição configuracional para o valor médio de uma propriedade A no ensemble NVT é dada por^4-6,12-14

onde k_B é a constante de Boltzmann e T a temperatura absoluta. A equação acima é uma média no espaço de fase configuracional t de 6 N dimensões. O termo (x^N) é a Hamiltoniana, omitindo os termos de energia cinética, constituída por termos de potencial intermolecular. O denominador na equação (3) é a função partição configuracional Z_NVT para esse ensemble.

Para o ensemble N_pT, em que a pressão do sistema é mantida constante e o volume é permitido variar^4-6,12-14

que contém um peso de Boltzmann modificado, e um espaço de integração de (6N+1) dimensões, incluindo o volume como variável adicional. O denominador na equação (4) é a função partição configuracional Z_NpT no ensemble N_pT.

Propriedades que podem ser expressas como médias configuracionais, tais como energia U, entalpia H e volume V, podem ser obtidas a partir das equações (3) e (4). Propriedades relacionadas às flutuações estatísticas de energia, entalpia e volume, também podem ser obtidas das derivadas dessas equações. Particularmente, podemos escrever a capacidade calorífica molar à volume e à pressão constantes, respectivamente, como^4,12,13

Podemos escrever também, em termos de flutuações^4,12,13, a compressibilidade isotérmica,

e o coeficiente de expansão isobárica,

As relações (5) a (7) requerem que sejam armazenadas durante a simulação quantidades como energia, entalpia e volume, bem como os quadrados destas quantidades.

A incerteza DA para uma média A, sobre N medidas, é calculada a partir da relação,

A incerteza para uma propriedade que já é dependente de média, por exemplo, energia livre, C_p, k_T ou a_P, é calculada como⁵,

onde <A>_b é a média de A calculada em um bloco b de configurações, N_b é o número de blocos e <A> é a média global de A.

2.3. Método de Monte Carlo com Algoritmo de Metropolis

O problema crítico de como calcular as médias, a partir das integrais configuracionais, foi resolvido por Metropolis e col.,⁹ implementando condições periódicas de contorno e um procedimento para a seleção de configurações no método de MC convencional.

Com as condições periódicas de contorno a amostra macroscópica, em nosso caso o líquido a ser modelado, é considerada infinitas réplicas de uma célula básica de referência. Essas condições de contorno são válidas para sistemas homogêneos. Nesses sistemas, podemos considerar que os movimentos feitos pelas moléculas nas células réplicas são os mesmos que ocorrem na célula de referência. A célula de referência contém um número finito N de moléculas, usualmente da ordem de algumas centenas, cuja densidade média é finita e a mesma da amostra macroscópica. A vantagem desse procedimento é reduzir um problema de alta dimensionalidade a um problema tratável computacionalmente, ao mesmo tempo que permite alcançar o limite termodinâmico.

A idéia básica do método de MC convencional é calcular numericamente as integrais envolvidas na equação (3), ou (4), aproximando-as por somas sobre um número n de pontos escolhidos aleatoriamente no espaço de fase t^5,6,14,

onde (x_n^N) corresponde a energia potencial total de N moléculas na n-ésima observação.

Na equação acima, cada configuração n tem a probabilidade P(x_n^N)de ser observada,^5,6,14

Entretanto, esse método de amostragem aleatória não é de muita utilidade para muitos problemas de interesse, já que o número n de observações é proibitivamente grande no espaço multidimensional. Ainda, para N moléculas na densidade do líquido, uma amostragem aleatória dos estados muitas vezes levaria a grande aproximidade ou mesmo superposição de duas moléculas. Nesse caso, o potencial de interação seria grande e positivo. Desse modo, estariamos incorrendo em problemas de convergência estatística, uma vez que com grande probabilidade escolheríamos configurações em que exp[- (x_n^N)/k _BT] é pequeno⁹.

Um procedimento mais eficiente para a seleção das configurações é restrito às regiões mais densas do espaço de fase e termodinamicamente mais prováveis. O algoritmo de amostragem modificado para o método de MC convencional, proposto por Metropolis e col., é denominado importance sampling e consiste em evitar a amostragem excessiva de configurações que tem pequena contribuição para o cálculo das médias. Vejamos, a seguir, como isso pode ser feito, examinando o algoritmo de Metropolis.

O algoritmo de Metropolis envolve calcular a probabilidade de que um movimento pode ser feito, caracterizando a transição de uma configuração npara uma configuração n^', através de um processo probabilístico, o pocesso de Markov. Em um processo de Markov é possível construir uma trajetória aleatória de pontos {x_n}, através do espaço de fase, independente de uma lei dinâmica que dite o comportamento físico do sistema. Nesse sentido, o estado a ser gerado depende unicamente do estado atual do sistema, sem informação a respeito da história passada do mesmo.

Este processo é definido por especificar uma probabilidade de transição,W (x^N_n ® x^N_n^, ), de um ponto n para outro ponto n^' no espaço de fase. O processo de Markov tem a propriedade de que a probabilidade P ({x^N_n}) converge para a probabilidade de equilíbrio P_eq ({x^N_n }), para um número muito grande de observações x^N_n, quando,^5,6,14

que significa que a razão das probabilidades de transição depende somente da mudança de energia entre os dois estados, n e n^', D = (x^N_n') - (x^N_n).

Se as observações incluem configurações em que realizam-se mudanças de volume, em adição a mudanças das posições das moléculas, temos, na representação do ensemble N_pT,^5,6,14

onde DV representa a mudança de volume envolvida na transição.

Assim, os pontos na trajetória evoluirão de acordo com um processo Markoviano segundo a probabilidade de transição (13), ou (14). A ordem seqüencial de estados gerados nesse processo não está, necessariamente, relacionada ao tempo físico. Dito de outra forma, as configurações geradas serão selecionadas de acordo com o peso de Boltzmann (13), ou (14), as quais podem ou não ser temporalmente consecutivas.

Vejamos agora como gerar as configurações tentativas que passarão pelo critério dos pesos de Boltzmann acima. Inicialmente, colocamos as N moléculas em uma configuração arbitrária em uma célula de referência, aqui denominada caixa de simulação. Uma configuração tentativa pode ser gerada, fazendo-se um movimento de translação e rotação combinados, para uma única molécula escolhida aleatoriamente. Se os sistemas estão sob as condições do ensemble N_pT, movimentos de volume também levam a novas configurações.

Podemos transladar o centro-de-massa da a-ésima molécula, ao longo de cada um dos três eixos cartezianos fixos no espaço, de acordo com a seguinte prescrição,⁵

onde os x_k são os números aleatórios, sorteados entre [0,1], e Dr_max é o máximo deslocamento permitido. Com essa prescrição, a molécula poderá deslocar-se de -Dr_max a +Dr_max na caixa de referência. De acordo com as condições periódicas de contorno, se a molécula, como um resultado dos movimentos acima, sair da caixa de referência, uma outra molécula da caixa réplica entrará pelo lado oposto, de modo a permanecer constante o número N de moléculas.

A rotação para o centro-de-massa da a-ésima molécula é feita em um dos três eixos cartezianos fixos no espaço e escolhido aleatoriamente. Uma vez escolhido um eixo para a rotação, por exemplo x, os deslocamentos angulares para os outros dois eixos cartezianos, y e z, podem variar entre -DQ_max e +DQ_max de acordo com a prescrição,⁵

A energia da configuração resultante dos dois movimentos conjuntos, de translação e rotação, será testada em relação a energia anterior, calculando-se a diferença de energia,

entre as duas configurações. Se a energia da nova configuração é menor do que a energia da configuração anterior, ou seja,

D < 0 (21)

o movimento é aceito pelo algoritmo de Metropolis. Caso contrário, um número aleatório p entre [0,1] é sorteado e comparado ao peso de Boltzmann (13). Se,

o movimento é aceito pelo algoritmo de Metropolis. Caso contrário, esse movimento é proibido pelo algoritmo, e retornamos a configuração anterior. A configuração anterior é, então, contada novamente para o cálculo das médias. A seguir, um novo movimento é tentado.

A repetição do procedimento acima fornecerá uma cadeia de configurações, distribuídas ao longo do espaço de fase configuracional de acordo com o peso de Boltzmann (13). A diferença fundamental nesse procedimento, em relação ao método de MC convencional, é que os pesos de Boltzmann são atribuídos às configurações antes das mesmas serem selecionadas. Assim, para as M configurações aceitas, a equação (11) reduz-se a,^5,6,9,14

onde A_n é o valor da propriedade A para o sistema, após o n-ésimo movimento ter sido aceito.

Podemos escrever o algoritmo de Metropolis, resumidamente, como,⁶

1. Especificar uma configuração inicial.

2. Gerar uma nova configuração x^N

3. Calcular a mudança de energia D.

4. Se D < 0, aceitar a nova configuração e retornar ao passo 2.

5. Se não, calcular exp[-D /k_BT]

6. Gerar um número aleatório p Î [0,1].

7. Se p < exp[-D /k_BT], aceitar a nova configuração e retornar ao passo 2.

8. Caso contrário, a configuração anterior passa ser a nova configuração e retornar ao passo 2.

Se o ensemble é o N_pT, o sistema também pode expandir ou contrair seu volume. Para uma caixa cúbica de referência as dimensões L_x, L_ye L_z são iguais, de modo que L³ = V. Serão tentados movimentos de -DV_max a +DV_max, de acordo com,⁵

Após o movimento de volume, todas as coordenadas do centro-de-massa das moléculas na caixa de referência são escalonadas pelo fator a = L'/L, de modo que o volume pode flutuar. Com issso, calcula-se a nova energia e testa-se com relação a anterior, seguindo-se o algoritmo acima, mas com o peso de Boltzmann modificado (14).

Os deslocamentos máximos, Dr_max e DQ_max, assim como os movimentos máximos de volume, DV_max, são parâmetros ajustáveis na simulação. Se esses parâmetros forem muito grandes, muitos movimentos serão proibidos. Por outro lado, se forem muito pequenos, as configurações não mudarão significativamente. Valores típicos para esses parâmetros são 0,15Å para Dr_max, 15^o para DQ_max e 200-600ų para DV_max.

Na grande maioria das simulações, tentativas de mover o volume não são feitas para cada nova configuração gerada. Esses movimentos são tentados com uma freqüência controlada, ao longo da simulação. Essa providência implica em uma redução em tempo computacional, uma vez que em um movimento de volume, todos os 1/2[N(N-1)] elementos da matriz de interação entre as N moléculas devem ser recalculados, ao passo que quando só uma molécula é movimentada por translação ou rotação somente uma linha desta matriz deve ser recalculada.

2.4. Amostragem Preferencial Próxima ao Soluto

Quando tratamos de soluções diluídas, é importante considerar uma amostragem diferenciada para as moléculas de solvente. Na amostragem de Metropolis, as moléculas de solvente são movimentadas com a mesma freqüência em toda a região da solução. Devemos considerar como regiões especialmente importantes das soluções aquelas próximas ao soluto, onde as interações soluto-solvente são mais fortes. No entanto, a maioria das moléculas de solvente encontram-se fora dessa região e interagem mais fracamente com o soluto. Se as moléculas de solvente são selecionadas de acordo com o procedimento de Metropolis, problemas estatísticos podem ocorrer devido a baixa concentração de soluto em relação a de solvente¹⁵. Para contornar esse problema, Owicki e Scheraga¹⁶ propuseram uma técnica de amostragem para o soluto e as moléculas de solvente a ele mais próximas. Essa técnica é denominada amostragem preferencial (preferencial sampling). Segundo essa técnica, além dos pesos de Boltzmann, as configurações são escolhidas segundo o peso multiplicativo adicional,

onde R é a distância entre o soluto e a molécula de solvente, e C é uma constante que assume um valor apropriado. Para um soluto poliatômico, R é calculado em relação a um sítio central da molécula.

A amostragem preferencial, definida pela equação (25), permite que as moléculas de solvente mais próximas ao soluto sejam amostradas mais freqüentemente em relação às moléculas de solvente mais distantes deste.

2.5. Função Distribuição Radial de Pares

Características da estrutura de um líquido podem ser obtidas através das funções distribuição radial de pares.

Como estrutura de um líquido entendemos a disposição média relativa entre as moléculas. A função distribuição radial, g_ij(r), entre os átomos i e j de um par de moléculas distintas, pode ser definida por⁵,

onde r é a separação entre os átomos, N_ij (r,r + Dr) é o número médio de átomos j encontrados na camada esférica entre r e r + Dr centrada no átomo i, 4pr²Dr é o elemento de volume da camada esférica, e r_j é a densidade numérica média de átomos j no líquido.

Pode-se dizer que g_ij (r) descreve a variação na distribuição dos átomos j em relação à distribuição que seria encontrada no líquido, para estes átomos, se o líquido fosse uniforme. Se a distribuição dos átomos j é a mesma que seria esperada de uma distribuição uniforme, N_ij(r,r + Dr)/4pr²Dr ® r_j, de modo que g_ij (r)®1. Esse comportamento é verificado para distâncias de separação r grandes, tipicamente da ordem de 10Å.

Experimentalmente, as funções distribuição radial podem ser determinadas a partir de técnicas de difração de raios-X, espalhamento de nêutrons e difração de elétrons^7,8. Em uma simulação computacional, o cálculo de g_ij (r) envolve, basicamente, o cálculo das distâncias r entre os átomos i e j. O número de vezes que ocorre o par de átomos em função da distância r é armazenado em um histograma. Esse histograma é incrementado ao longo da simulação, e por fim, é normalizado para o número N de moléculas utilizadas na caixa de simulação e para o número de passos usados no cálculo de g_ij (r).

2.6. Números de Coordenação

Os máximos das funções distribuição radial representam as distâncias médias entre vizinhos e, no caso de líquidos, são associados às camadas de solvatação. A integração da equação (26), até os mínimos correspondentes, fornece estimativas do número de vizinhos, ou números de coordenação C_ij. Por exemplo, para a primeira camada de solvatação com um mínimo em R_c temos,

onde C_ij(R_c ) é o número de primeiros vizinhos.

2.7. Função Potencial de Interação

Em vista do procedimento de amostragem, é essencial conhecer a energia de uma configuração. Primeiramente, em simulação de líquidos, as moléculas são representadas por um conjunto de sítios, distribuidos sobre a estrutura molecular, e a energia de interação é obtida como uma soma de interações sítio-sítio. Para a grande maioria das aplicações na literatura, as distâancias de ligação e os ângulos de ligação são mantidos fixos, de forma que a geometria das moléculas permanece rígida ao longo da simulação. Portanto, a representação comumente usada para a energia de uma configuração consiste da contribuição intermolecular, desprezando-se contribuições intramoleculares devido a vibrações e rotações internas das moléculas. Essas contribuições internas para a energia podem ser incluídas sem modificar o esquema geral de simulação, apresentado acima. O termo intermolecular, usualmente, é representado por uma soma sobre todos os pares de moléculas a e b,

A energia de interação, E_ab, é dada como uma soma de potenciais de Lennard-Jonnes e Coulomb entre o sítio i da molécula a e o sítio j da molécula b,

onde e₀ é a permissividade do vácuo e r_ij é a distância entre os sítios i e j, em Angstrom. Os parâmetros A_kk e B_kk são expressos, respectivamente, como

onde o conjunto de parâmetros (s_kk, e_kk) representa os parâmetros de Lennard-Jones para o sítio k, e_kk, em unidades de kcal mol^-1, e s_kk, em Angstrom. Em geral, para economizar tempo na dedução de parâmetros, os termos A_ij e B_ij são calculados através das regras de cruzamento,⁵

e usados para obter os parâmetros de Lennard-Jones não diagonais. As cargas q_k, na equação (29), são cargas parciais fixas centradas nos sítios, em unidades de elétron, e em geral coincidem com as posições dos átomos nas moléculas. A distribuição de cargas pontuais sobre a estrutura molecular é usada para levar em conta a contribuição da assimetria da densidade de carga molecular para a energia intermolecular. Na seção 4, apresentaremos uma metodologia eficiente para obter cargas parciais adequadas para serem utilizadas em simulação molecular em fase condensada.

Para sistemas sujeitos a condições periódicas de contorno, há dois procedimentos para gerar novas configurações a partir do cálculo da energia potencial de uma dada configuração⁵, os quais são: (a) a convenção de imagem mínima, e (b) a convenção de imagem mínima com um raio de corte esférico adicional.

O procedimento (b) emprega um raio de corte esférico, r_c, de modo que a interação entre duas moléculas separadas por uma distância maior do que r_c é desprezada. Em geral, r_c é a metade da caixa cúbica de simulação. Então, para uma dada molécula, somente interações com moléculas na esfera vizinha com raio igual a metade do cubo são calculadas. Para distâncias além do raio de corte, r > r_c, as contribuições para a energia são calculadas através da equação⁵,

Na equação (32), assumimos que g_ij(r) = 1 além do raio de corte, e que somente os termos de Lennard-Jones, u_ij, contribuem para o potencial (28). Segundo Keesom¹⁷, a interação dipolo-dipolo, u_d-d, após a média sobre a orientação angular, comporta-se assintoticamente como uma função radial, ou seja,

onde m_i é o momento de dipolo permanente da molécula i e r é a distância entre os centros dos dipolos. Logo, a correção de longo alcance para a interação dipolo-dipolo pode ser introduzida na expressão analítica utilizada para a correção de longo alcance do potencial de Lennard-Jones, apresentada na equação (32). O comportamento assintótico encontrado na equação acima está em acordo com resultados obtidos em estudos de efeitos de tamanho em sistemas compostos por moléculas neutras^18-20, os quais indicam que a contribuição das interações Coulombicas além do raio de corte é muito pequena.

Para maior estabilidade numérica, o potencial de interação é multiplicado por uma função de ajuste²¹ no intervalo [(r_c-0,5)-r_c] Å. Com essa função,o potencial de interação é suavemente reduzido a zero nos limites da esfera de interação definida pelo raio de corte. Dessa forma, são evitados saltos abruptos no valor de energia de interação quando uma molécula entra ou sai da região de interação definida pelo raio de corte.

A convenção de imagem mínima para calcular o potencial de interação é utilizada em ambos os procedimentos, (a) e (b). Com a convenção de imagem mínima, as interações entre uma dada molécula e as imagens mais próximas das outras N - 1 moléculas na caixa de simulação são calculadas. Com isso, o cálculo da energia potencial, devido às interações entre todos os pares de moléculas, envolve 1/2N (N - 1) termos. Para cada molécula, N - 1 interações são calculadas. Assim, para uma caixa de simulação contendo N=100 moléculas e um total de 10⁸ configurações, um total de 10⁸ x 99 »10¹⁰ interações entre pares de moléculas seriam calculadas, uma vez que somente uma molécula é movimentada por configuração. A aproximação de imagem mínima deixa a simulação mais lenta por um fator de 2 em relação a aproximação que utiliza um raio corte esférico adicional. Isso porque para cada molécula são calculadas N - 1 interações, no caso de imagem mínima, e aproximadamente N/2 interações quando um raio de corte esférico é empregado.

3. SIMULAÇÕES DE ENERGIA LIVRE

Para um sistema de N moléculas, descrito pela Hamiltoniana H(x^N), a energia livre de Gibbs pode ser escrita como^4,12,13,

G = ¾ k_BTlnZ_NpT, (34)

onde Z_NpT está definida na equação (4). Uma expressão equivalente pode ser obtida para a energia livre de Helmholtz, com a função partição dada em (3).

Como já exposto, propriedades configuracionais de equilíbrio podem ser obtidas seguindo-se os passos contidos na seção (2.3),mas perdemos informação no denominador das equações (3) e (4). Então, não temos conhecimento das funções partição envolvidas. Como conseqüência, as energias livres não são prontamente obtidas em um cálculo de MC com o método de Metropolis, o qual amostra somente uma região restrita do espaço de fase, a região de mais baixa energia.

O problema de estimar a energia livre, através de simulação de MC com o método de Metropolis, tem sido discutido na literatura^22-28. Foram introduzidas metodologias especiais que fornecem meios para se obter estimativas para as diferenças de energia livre entre dois sistemas previamente escolhidos. As diferenças de energia livre, assim obtidas, podem ser diretamente comparadas com resultados experimentais. Dentre as metodologias existentes, passaremos a discutir a teoria de perturbação termodinâmica, introduzida por Zwanzig²⁹, a qual tem sido largamente utilizada^20,30-37.

3.1. Teoria de Perturbação Termodinâmica

O método perturbativo pode ser usado para calcular a diferença de energia livre de Gibbs, DG, entre dois estados A e B, bem definidos, relacionados pela perturbação DH da forma,

_B =

_A + D. (35)

Assim, a expressão para DG, que acompanha a transformação entre o estado inicial não perturbado A (ou de referência) e o estado final B, pode ser escrita como,

ou, usando as equações (34) e (35),

O aspecto interessante dessa equação é que a primeira exponencial do numerador, dividida pelo denominador, corresponde à distribuição de probabilidade para a amostragem baseada na Hamiltoniana do estado inicial A. Com isso, a equação (37) pode ser expressa como uma média no ensemble das configurações representativas do estado inicial,

A equação (38) é a essência do método perturbativo. Uma vez que somente diferenças de energia livre são envolvidas, os estados inicial e final são intercambiáveis. Assim, podemos escrever,

onde a média amostral <...>_B é baseada na Hamiltoniana do estado final B.

Uma vez definidos

_A e _B, o método de MC pode ser usado para fazer a média dos termos exp (-D /k_BT), calculados para cada configuração. Um procedimento análogo pode ser seguido com respeito à equação (39). Se a equação (38) é aplicada, o cálculo da mudança de energia livre, acompanhando o processo A ® B, é baseado em uma amostragem do sistema inicial A. Se, alternativamente, a equação (39) é usada para determinar a diferença de energia livre na transformação B ® A, a fase de equilíbrio e de média na simulação deve ser sobre estados característicos do estado B, representado pela Hamiltoniana _B.

Estados intermediários entre os estados A e B podem ser definidos, fazendo-se a Hamiltoniana do sistema uma função de um parâmetro de acoplamento l, que toma valores discretos no intervalo [0,1], tal que,

=

(l = 0) =

_A, e, =

(l = 1) =

_B. (40)

Assim, o i-ésimo estado intermediário pode ser definido pela Hamiltoniana (l_i), de acordo com a equação,

(l_i) = l_i

_B + (1 - l_i)_A. (41)

Estados como (l_i+ Dl_i) = (l_i+1) e (l_iDl_i) = (l_i-1) são aqueles imediatamente posterior e anterior ao estado (l_i), respectivamente.

Quando l variar de forma adequada no intervalo [0,1], o sistema inicial, caracterizado por

_A, é gradualmente transformado no sistema final, caracterizado por _B. Com esse procedimento, a diferença de energia livre entre o subintervalo l_i e l_i+1 pode ser obtida como,

DG(l_i ® l_i+1) = -k_BT ln < exp [-D (l_i ® l_i+1)/k_B T]>_li, (42)

onde, D (l_i ® l_i+1) = (x^N, l_i+1) - (x^N, l_i), e em sentido oposto,

DG(l_i ® l_i-1) = -k_BT ln < exp [D (l_i ® l_i-1)/k_B T]>_li, (43)

onde, D (l_i ® l_i-1) = H (x^N, l_i-1)- H (x^N, l_i). É importante observar que nas equações (42) e (43) as amostragens são baseadas no estado intermediário l_i. Com essas equações, as diferenças de energia livre para quaisquer dois subintervalos subseqüentes, interiores ao intervalo [0,1], representados por (x^N, l_i-1) e (x^N, l_{i +1}), podem ser obtidas a partir de uma única simulação, equilibrada para o estado intermediário representado por (x^N, l_i).As duas médias resultantes desses cálculos para obter DG(l_i ® l_i+1) e DG(l_i ® l_i-1), por meio de uma única amostragem para ambas as direções, é conhecida na literatura como amostragem de passo duplo^22,38 (double-wide sampling). Quando a amostragem é realizada em uma única direção, dizemos que a amostragem é de passo único.

A variação total de energia livre é obtida somando-se as contribuições de cada subintervalo na variável l_i,

Como a energia livre é uma função de estado, DG_total não depende da escolha do parâmetro l. A dependência funcional de H(x^N), sobre l, é obtida fazendo-se todos os termos do potencial de interação uma função de l. Assim, a função potencial varia, durante a perturbação, quando um sistema é transformado em outro. Se z é um conjunto de parâmetros de potencial para um dado sítio molecular, então z varia, quando l_i varia de 0 a 1, de acordo com,

z(l_i) = l_iz_B+ (1 - l_i)z_A, (45)

onde i, A e B são os estados de referência, inicial e final, respectivamente.Parâmetros geométricos, como comprimento de ligação e ângulo de ligação, também variam de acordo com a equação acima, se, na transformação A ® B, a geometria dos componentes também é modificada.

4. FUNÇÃO POTENCIAL DE INTERAÇÃO: CARGAS ATÔMICAS DERIVADAS DO POTENCIAL ELETROSTÁTICO MOLECULAR - O MÉTODO CHELPG

A representação usual e mais simples das propriedades eletrostáticas de líquidos é através do potencial gerado por cargas pontuais, centradas em sítios moleculares. Como o valor médio de uma carga sobre um átomo não é uma quantidade rigorosamente definida em mecânica quântica, várias técnicas têm sido propostas e discutidas^39-52, para respresentar a densidade de carga molecular por um conjunto de cargas pontuais.

Um método de obtenção de cargas pontuais concebido para representar adequadamente o potencial eletrostático Coulombico intermolecular é conhecido na literatura como método de cargas derivadas do potencial, ou método CHELPG (Charges from Electrostatic Potential Grid based), proposto por Breneman e Wiberg⁴³. Cargas derivadas do potencial incluem: (i) a determinação de uma função de onda apropriada; (ii) a determinação do potencial eletrostático em vários pontos em torno da molécula; e (iii) o cálculo das cargas, em posições pré-definidas, por meio de um ajuste de mínimos quadrados dos potenciais eletrostáticos, calculados classicamente, aos valores dos potenciais calculados quanticamente.

Formalmente, o potencial eletrostático molecular para um ponto i, na posição x_i, em torno de uma molécula, composta de N_a átomos e N_e elétrons, pode ser definido como a interação entre a distribuição de carga molecular e uma carga positiva unitária nesta posição x_i,^50-53,

onde Z_j é a carga nuclear sobre o átomo j, centrada na posição R_j, e r é a densidade eletrônica. O primeiro termo na equação (46) corresponde à repulsão entre as cargas pontuais Z_j e a carga unitária em r_i, e o segundo termo, à atração eletrostática envolvendo a distribuição de cargas eletrônicas em todo o espaço e a carga unitária positiva em r_i. No formalismo de orbitais moleculares^39-41 (LCAO - MO, Linear Combination of Atomic Orbitals - Molecular Orbitals), o potencial eletrostático molecular é o valor esperado do operador 1/r, e a equação (46) pode ser reescrita como,

Na equação (47), a somatória em k estende-se sobre o número de orbitais moleculares, duplamente ocupados, {X_r} é o conjunto de orbitais atômicos ou funções de base, e {C_rk} são os coeficientes da combinação linear de funções de base {X_r} para o k-ésimo orbital molecular

O potencial eletrostático molecular pode, também, ser representado por uma expansão em multipolos atômicos. Na aproximação de monopolo, o potencial eletrostático no ponto i pode ser calculado, classicamente, a partir de cargas atômicas pontuais q_j,

A densidade de carga molecular pode, então, ser particionada em cargas pontuais q_j entre os átomos constituintes, na condição de que E_i(q₁,q ₂,...,q_Na) reproduza V_i, calculado segundo a equação (47). Para isso, o critério de ajuste de mínimos quadrados é utilizado, de modo a minimizar a função Y,

onde N_p é o número total de pontos em que o potencial eletrostático é calculado. A qualidade do ajuste do potencial eletrostático, gerado classicamente, com aquele determinado quanticamente, é estimada em termos dos parâmetros desvio quadrático médio, RMS = (Y/N_p)^1/2, e desvio quadrático médio relativo,

O método implementado no programa CHELPG para resolver as equações acima usa multiplicadores de Lagrange, com a seguinte condição de vínculo sobre as cargas,

onde q_tot é a carga total da molécula. A função a ser minimizada, com o vínculo dado pela equação (50), é Z = Y + gG, onde g é o multiplicador de Lagrange. O mínimo de Z e os correspondentes q_j's são encontrados resolvendo-se o sistema de (N_a + 1) equações,

que resulta em,

Com a notação r_ik = |r_i - R_k |, e definindo,

tem-se a equação matricial,

Em notação matricial,

B = Aq. (55)

Resolvendo-se para q,

q = A^-1B, (56)

tem-se o conjunto de cargas {q_j }.

As cargas {q_j} são cargas efetivas com o único vínculo de conservar a carga molecular, representado pela equação (50). Generalizações podem ser feitas no procedimento de ajuste, por exemplo, forçar as cargas a reproduzir não só a carga molecular como também o momento de dipolo⁵¹, o que representaria introduzir mais três equações de vínculo, uma para cada componente do momento de dipolo da molécula.

O conjunto de cargas CHELPG obtido acima é o que melhor reproduz o potencial eletrostático molecular em regiões além da superfície de van der Waals da molécula. Essa característica, aliada à constatação de que o potencial eletrostático calculado nessas regiões apresenta pouca dependência com o conjunto de base e o nível de cálculo empregados, confere às cargas derivadas do potencial uma técnica eficaz e viável computacionalmente para representar potenciais eletrostáticos intermoleculares. Esse método tem produzido resultados bastante promissores quando aplicado à simulação molecular em fase condensada, em particular, à simulação de líquidos^20,54,55. Isso, em parte, devido à representação efetiva que o conjunto de cargas CHELPG fornece para os efeitos de polarização existentes em fase condensada, os quais podem ser parcialmente incorporados ao conjunto de cargas fixas no campo de força, como aquele representado pela equação (29). De certa forma, esses efeitos de polarização estão incluídos nos modelos OPLS (Optimized Potential for Liquid Simulation), desenvolvidos por Jorgensen e col.^30,56-60. para estudar propriedades termodinâmicas de líquidos.

5. CONCLUSÕES

Neste trabalho, apresentamos as equações básicas da mecânica estatística e sua implementação através do método de Monte Carlo para estudar propriedades termodinâmicas de líquidos. A implementação da teoria de perturbação termodinâmica com o método de Monte Carlo para calcular variações de energia livre em processos envolvendo interação soluto-solvente foi também discutida. Neste trabalho, as interações intermoleculares são descritas por potenciais clássicos, obtidos pela composição de funções de Lennard-Jones e Coulomb. Uma metodologia para obter o potencial eletrostático efetivo a partir da função de onda molecular foi apresentada. Esse potencial efetivo é gerado por um conjunto de cargas pontuais distribuidas sobre os sítios moleculares. A potencialidade das metodologias apresentadas para calcular propriedades termodinâmicas e distribuições radiais de pares para líquidos complexos é também abordada.

AGRADECIMENTOS

À FAPESP e ao CNPq pelo apoio financeiro.

REFERÊNCIAS

1. Barke, J. A.; Hendersen, D.; Rev. Mod. Phys. 1976, 48, 587.

2. Croxton, C. A.; ed.; Progress in Liquids Physics, Wiley, New York 1978.

3. Ailawadi, N. K.; Phys. Rep.; Rev. Sec. Phys. Lett. 1980, 57, 241.

4. Hansen, J. P.; McDonald, I. R.; Theory of Simple Liquids, Academic Press, London 1986.

5. Allen, M. P.; Tildesley, D. J.; Computer Simulation of Liquids, Clarendon Press, Oxford 1987.

6. Heerman, D. W.; Computer Simulations Methods, Springer-Verlag, Berlin 1986.

7. Teixeira-Dias, J. J.; ed.; Molecular Liquids: New Perspectives in Physics and Chemistry, Kluwer Acad. Pub.; Dordrecht 1992.

8. Evans M. W.; Dynamical Process in Condensed Matter, Wiley, New York 1985.

9. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E.; J. Chem. Phys. 1953, 21, 1087.

10. Alder, B. J.; Wainwright, T. E.; J. Chem. Phys. 1957, 27, 1208.

11. Alder, B. J.; Wainwright, T. E.; J. Chem. Phys. 1959, 31, 459.

12. Hill, T. E.; An Introduction to Statistical Mechanics, Dover, New York 1986.

13. McQuarrie, D. A.; Statistical Mechanics, Harper and Row, New York 1976.

14. Binder, K.; ed.; Monte Carlo Methods in Statistical Physics, Springer-Verlag, Berlin, 2^nd edition 1986.

15. Jorgensen, W. L.; Bigot, B.; Chandrasekhar, J.; J. Am. Chem. Soc. 1982, p. 4584.

16. Owicki, J. C.; Scheraga, H. A.; Chem. Phys. Lett. 1977, 47, 600.

17. Israelachivili, J. N.; Intermolecular and Surface Forces, Academic Press, London 1985.

18. Freitas, L. C. G.; Cordeiro, J. M. M.; J. Mol. Struct. (THEOCHEM) 1995, 335, 189.

19. Freitas, L. C. G.; J. Mol. Struct. (THEOCHEM) 1993, 282, 151.

20. Barlette, V. E.; Garbujo, F. L. L.; Freitas, L. C. G.; J. Mol. Engineering, special issue in Molecular Modeling 1997, 7, 439.

21. Andrea, T. A.; Swope, W. C.; Andersen, H. C.; J. Chem. Phys. 1983, 4576.

22. Gunsteren, W. F. V.; Weiner, P. K., ed.; Computer Simulation of Biomolecular Systems: Theoretical and Experimental Aplications, ESCOM, Leiden 1989.

23. Beveridge, D. L.; DiCapua, F. M.; Annu. Rev. Biophysics 1989, 18, 431.

24. Kollman, P.; Chem. Rev. 1993, 93, 2395

25. Bash, P. A.; Singh, U. C.; Langridge, R.; Kolman, P.; Science 1987, 236, 564.

26. Reynolds, C. A.; King, P. M.; in Richards, W. G.; ed.; Computer-Aided Molecular Design, Cap. 4, IBC Technical Services Ltd.; Great Britain.

27. Mezei, M.; Beveridge, D. L.; Ann. N. Y. Acad. Sci. 1986, 482, 1.

28. Bennet, C. H.; J. Comp. Phys. 1976, 22, 245

29. Zwanzig, R. W.; J. Chem. Phys. 1954, 22, 1420.

30. Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L.; J. Phys. Chem. 1990, 94, 1683.

31. Jorgensen, W. L.; Buckner, J. K.; J. Phys. Chem. 1987, 91, 6083.

32. King, P.; Reynolds, C. A., Essex, J. W.; Worth, G. A.; Richards, W. G.; Mol. Simul. 1990, 5, 265.

33. Jorgensen, W. L.; Ravimohan, C.; J. Chem. Phys. 1985, 83, 3050.

34. Jorgensen, W. L.; Blake, J. F.; Bukner, K.; Chem. Phys. 1989, 129, 193.

35. Jorgensen, W. L.; Severance, D. L.; J. Am. Chem. Soc. 1990, 112, 4768.

36. Freitas, L. C. G.; Silva, L. B.; Botelho, L. F.; Quím. Nova 1996, 19, 166.

37. Freitas, L. C. G.; Botelho, L. F.; Quím. Nova 1994, 17, 489

38. Swaminathan, S.; Whitehead, R. J.; Guth, E.; Beveridge, D. L.; J. Am. Chem. Soc. 1997, 99, 7817.

39. Henre, W. J.; Radom, L.; Schleyer, P. V.; Pople. J. A.; Ab Initio Molecular Orbital Theory, John Wiley & Sons, New York 1986.

40. Szabo, A.; Ostlund, N. S.; Modern Quantum Chemistry, MacMillan, New York, 2^nd edition 1982.

41. Cook, D. B.; An Initio Valence Calculations in Chemistry, Buttherworth, London 1974.

42. Chirlian, L. E.; Francl, M. M.; J. Comp. Chem. 1987, 8, 894

43. Breneman, C. M.; Wiberg, K. B.; J. Comp. Chem. 1990, 11, 361.

44. Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A.; J. Phys. Chem. 1993, 97, 10269.

45. Rappé, A. K.; Goddard III, W. A.; J. Phys. Chem. 1991, 95, 3358.

46. Bader, R. F. W.; Acc. Chem. Res. 1985, 18, 9.

47. Wiberg, K. B.; Rablen, P. R.; J. Comp. Chem. 1993, 14, 1504.

48. Guadagnini, P. H.; Bruns, R. E.; Souza, A. A. D.; Quím. Nova 1996, 19. 148.

49. Merz, Jr.; K. M.; J. Comp. Chem. 1992, 13, 749.

50. Alemán, C.; Luque, F. J.; Orozco, M.; J. Comp. Chem. 1993, 14, 799.

51. Woods, R. J.; Khalil, M.; Pell, W.; Moffat, S. H.; Smith, Jr.; V. H.; J. Comp. Chem. 1990, 11, 297.

52. Besler, B. H.; Merz Jr., K. M.; Kollman, P. A.; J. Comp. Chem. 1990, 11, 431.

53. Szabó, G. N.; Ferenczy, G. G.; Chem. Rev. 1995, 95, 829.

54. Carlson, H. A.; Nguyen, T. B.; Orozco, M.; Jorgensen, W. L.; J. Comp. Chem. 1993, 14, 1240.

55. Sinoti, A. L. L.; Politi, J. R. S.; Freitas, L. C. G.; J. Braz. Chem. Soc. 1996, 7, 133.

56. Jorgensen, W. L.; Chandrasenkhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.; J. Chem. Phys. 1983, 79, 926.

57. Jorgensen, W. L.; Tirado-Rives, J.; J. Am. Chem. Soc. 1988, 110, 1657.

58. Jorgensen, W. L.; Madura, J. D.; Swenson, C. J.; J. Am. Chem. Soc. 6638.

59. Jorgensen, W. L.; J. Phys. Chem. 1986, 90, 1276.

E-mail: [e-mail]

Datas de Publicação

• Publicação nesta coleção
  12 Maio 2000
• Data do Fascículo
  Abr 1999

Histórico

• Aceito
  14 Abr 1998
• Recebido
  18 Fev 1997