UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE INSTITUTO DE QUÍMICA CURSO DE QUÍMICA CAROLINE RIBEIRO CARNEIRO SIMULAÇÃO COMPUTACIONAL DOS COEFICIENTES DE AUTODIFUSÃO E DE VISCOSIDADE DO TETRACLOROALUMINATO DE 1-ETIL-3-METILIMIDAZÓLIO Natal / RN, 2015. 2 SIMULAÇÃO COMPUTACIONAL DOS COEFICIENTES DE AUTODIFUSÃO E DE VISCOSIDADE DO TETRACLOROALUMINATO DE 1-ETIL-3-METILIMIDAZÓLIO Por Caroline Ribeiro Carneiro Trabalho de conclusão de curso apresentado à Coordenação do curso de Química da Universidade Federal do Rio Grande do Norte, como Requisito Final à obtenção do título de Bacharel em Química. Orientador: Prof. Dr. Jones de Andrade. Natal / RN, 2015. 3 4 Caroline Ribeiro Carneiro Simulação Computacional dos Coeficientes de Autodifusão e de Viscosidade do Tetracloroaluminato de 1-Etil-3-Metilimidazólio Trabalho de conclusão de curso apresentado à Coordenação do curso de Química da Universidade Federal do Rio Grande do Norte, como Requisito Final à obtenção do título de Bacharel em Química. Orientador: Prof. Dr. Jones de Andrade. Monografia aprovada em: 17 / 12 / 2015. COMISSÃO EXAMINADORA ___________________________________________ Prof. Dr. Jones de Andrade (Orientador) Universidade Federal do Rio Grande do Norte _____________________________________________ Prof.ª Drª. Sibele Berenice Castella Pergher (Membro) Universidade Federal do Rio Grande do Norte ___________________________________________ Prof. Dr. Davi Serradella Vieira (Membro) Universidade Federal do Rio Grande do Norte 5 AGRADECIMENTOS Agradeço, em primeiro lugar, a Deus por me dar forças para superar as dificuldades. Aos meus pais, Albanisa e Francisco, por sempre me apoiarem, acreditarem nos meus sonhos e pelo seu amor incondicional, vocês foram fundamentais para que eu chagasse até aqui. A minha irmã e melhor amiga Thaise Nayane, obrigada por ser meu referencial e estar sempre ao meu lado. A minha tia Albinha por todos os conselhos e apoio. E a Alice, Eva e Nuvem, as melhores coelhas do mundo, pelo seu amor mais puro e genuíno. Ao professor Jones, por toda sua paciência, dedicação e ensinamentos. Aos colegas do laboratório LAMMSA. E a todos que de alguma maneira, direta ou indiretamente, fizeram parte da minha formação, meu muito obrigado. “Seja quem você for, [...] tenha sempre como meta muita força, muita determinação e sempre faça tudo com muito amor e com muita fé em Deus, que um dia você chega lá. De alguma maneira você chega lá”. Ayrton Senna. 6 RESUMO Fazendo uso da simulação computacional por dinâmica molecular, foram estudas propriedades de transporte (viscosidade e coeficiente de difusão), para o líquido iônico tetracloroaluminato de 1-etil-3-metilimidazólio. Avaliou-se o efeito do ajuste de cargas no sistema, sobre essas propriedades, por meio de três modelos: sem transferência de carga, com transferência de carga mas sem rearranjo, e também tanto com transferência de carga como com rearranjo. As simulações computacionais foram executadas com o ensemble isotérmico-isobárico (NpT), utilizando o campo de força AMBER e o pacote de programas GROMACS. Dentre os métodos, modelos e configurações estudados, o modelo apenas com transferência de cargas apresentou o comportamento difusivo mais adequado em tempos longos, enquanto que o modelo com transferência e rearranjo apresentou os melhores resultados no que diz respeito ao coeficiente de autodifusão. Já os valores de viscosidade obtidos indicam que o modelo sem transferência ou rearranjo de cargas é o que mais se aproximou do valor apresentado na literatura. As verificações realizadas também indicam uma forte relação entre as propriedades macroscópicas estudadas e as interações eletrostáticas provenientes do cátion e ânion para os líquidos iônicos. Palavras-chave: Líquidos iônicos; dinâmica molecular; viscosidade; coeficiente de difusão; transferência de carga; 7 ABSTRACT Using computer simulation by molecular dynamics, the ionic liquid 1-ethyl-3- methylimidazolium tetrachloroaluminate transport properties (viscosity and diffusion coefficient) were studied. The charge adjustment effects on these properties of the system were evaluated by three models: without charge transference; with charge transference but without rearrangement; and also with both charge transference and rearrangement. Computer simulations were performed within the isothermal-isobaric ensemble (NpT), using AMBER force field and GROMACS program package. Among the methods, designs and configurations studied, the model only with charge transference proved to be most adequate in relation to diffusion behavior in long times, while the model with both charge transfer and rearrangement showed the best results regarding the self-diffusion coefficient. Finally, the shear viscosity obtained values indicate that the model with neither charge transfer or rearrangement was most closely related to the value presented in the literature. The performed verifications also indicate a strong relationship between the studied macroscopic properties and the ionic liquids cation and anion electrostatic interactions. Keywords: Ionic liquids; molecular dynamics; viscosity; diffusion coefficient; charge transference. 8 ÍNDICE LISTA DE FIGURAS LISTA DE TABELAS LISTA DE ABREVIATURAS E SIGLAS 1 INTRODUÇÃO ..................................................................... 12 1.1 LÍQUIDOS IÔNICOS ........................................................... 12 1.1.1 Tetracloroaluminato de 1-etil-3-metilimidazólio.................... 14 1.2 DINÂMICA MOLECULAR E ALGORITMOS ........................ 16 1.3 POTENCIAIS E CAMPOS DE FORÇA ................................ 18 1.3.1 Polarização............................................................................ 20 1.4 PROPRIEDADES DE TRANSPORTE ................................. 21 1.4.1 Difusão ................................................................................. 21 1.4.2 Viscosidade........................................................................... 23 2 METODOLOGIA................................................................... 26 2.1 GERAÇÃO DA CONFIGURAÇÃO INICIAL DAS MOLÉCULAS ....................................................................... 26 2.2 EQUILIBRAÇÃO DO SISTEMA ........................................... 26 2.3 MODELOS DE TRANSFERÊNCIA ESTÁTICA DE CARGAS PERTURBANDO O CAMPO DE FORÇA ............ 30 2.4 AQUISIÇÃO DE TRAJETÓRIA E ANÁLISES ...................... 33 3 RESULTADOS E DISCUSSÃO .......................................... 34 3.1 Coeficiente de Autodifusão .................................................. 34 9 3.2 Coeficiente de Viscosidade .................................................. 42 4 CONCLUSÃO ...................................................................... 49 5 REFERÊNCIAS ................................................................... 51 ANEXOS .............................................................................. 54 10 LISTA DE FIGURAS Figura 1 – Estrutura molecular do cátion 1-etil-3-metilimidazólio............................. 14 Figura 2 – Estrutura molecular do ânion tetracloreto de alumínio.............................15 Figura 3 - Ilustração das condições de contorno periódicas.....................................17 Figura 4 – Fluxograma executado no ajuste da densidade do sistema....................27 Figura 5 - Variação da densidade do sistema ..........................................................28 Figura 6 - Número associado a cada átomo da estrutura molecular ........................31 Figura 7 – MSD (a) do cátion EMI+ e (b) do ânion AlCl -4 nos 3 modelos .................35 Figura 8 – MSD (a) do cátion EMI+ e (b) do ânion AlCl -4 em escala logarítmica para os 3 modelos .............................................................................................................36 Figura 9 – VACFs para os 3 modelos estudados do cátion EMI+, apresentadas no intervalo de tempo de correlação (a) de 0 à 2,5 ps e (b) de 2 à 25 ps ......................40 Figura 10 – VACFs para os 3 modelos testados do ânion AlCl -4 , apresentadas no intervalo de tempo de correlação (a) de 0 à 2,5 ps e (b) de 2 à 25 ps.......................41 Figura 11 – Relação de Einstein para viscosidade calculada a partir das flutuações dos termos fora da diagonal do tensor de pressão nos 3 modelos estudados..........43 Figura 12 – TCAFs global e seccionados para (a) o modelo 0, (b) o modelo 1 e (c) modelo 2.....................................................................................................................45 Figura 13 – TCAFs globais dos 3 modelos estudados..............................................48 11 LISTA DE TABELAS Tabela 1 – Barostato, termostato e número de passos empregados nas 3 subetapas da segunda fase de equilibração do sistema ............................................................29 Tabela 2 – Distribuição de carga parcial nos átomos................................................32 Tabela 3 – Resultados obtidos das análises de MSD do cátion ...............................37 Tabela 4 - Resultados obtidos das análises de MSD do ânion ................................38 Tabela 5 – Parâmetros do ajuste das TCAFs global e seccionadas do Modelo 0....46 Tabela 6 – Parâmetros do ajuste das TCAFs global e seccionadas do Modelo 1....46 Tabela 7 – Parâmetros do ajuste das TCAFs global e seccionadas do Modelo 2....47 Tabela 8 – Coeficientes de viscosidade obtidos pelos 3 modelos estudados, calculados pelas suas TCAFs....................................................................................47 12 LISTA DE ABREVIATURA E SIGLAS AMBER – campo de força utilizado, do inglês Assisted Model Building with Energy Refinement. EMI+ - 1-etil-3-metilimidazólio. MD – Dinâmica molecular, do inglês Molecular Dinamics. PBC - Condições periódicas de contorno, do inglês Periodic Boundary Conditions. PIMs – modelo de íons polarizáveis, do inglês Polarizable Ions Model. IL – líquido iônico, do inglês Ionic Liquid. RESP – Método para cálculo de cargas pontuais com base no potencial eletrostático molecular, do inglês Restrained ElectroStatic Potencial. RIMs – modelo de íons rígidos, do inglês Rigid Ions Model. TCAF – função de autocorrelação de correntes transversais, do inglês Transverse- current Correlation Functions. MSD – deslocamento médio quadrático, do inglês Mean Square Displacement. VACF – função de autocorrelação da velocidade, do inglês Velocity Autocorrelation Function. 13 1. INTRODUÇÃO 1.1 Líquidos Iônicos Tradicionalmente, uma substância é líquida na temperatura ambiente quando as interações entre seus átomos constituintes são mais fortes, e, portanto, mais favoráveis, que aquelas que podem existir no estado gasoso, todavia sendo ainda mais fracas que as interações que constituem o estado sólido. Moléculas constituídas por íons possuem entre si uma forte atração: desse modo, costumam se agregar no estado sólido da matéria. Logo, pode-se dizer que grande parte das substâncias líquidas são constituídas de partículas neutras. (CONSORTI et al., 2001). Porém, averiguações mais profundas sobre a natureza das interações moleculares e ligações químicas constataram a existência de sais líquidos: substâncias formadas por íons negativos e positivos mas em estado líquido, nos quais a força de coesão entre os íons é predominantemente eletrônica. (FOGAÇA et al., 2010) Em meados da década de 1940, foi descoberto que ao combinar-se o cloreto de alquilpiridínio e o tricloreto de alumínio, esses adquiriam a capacidade de compor um sistema iônico com baixa temperatura de fusão. Já em 1963, Yoke e colaboradores produziram um composto com base no ânion diclorocuprato (CuCl -2 ) e no cátion trietilamônio (NH(C +2H5)3 ), que por sua vez eram obtidos pela mistura de cloreto cuproso (CuCl) e cloreto de trietilamônio (NH(C2H5)3.Cl): tratava-se de um composto iônico que era oleoso e não volátil a temperatura ambiente. Em seguida, foi divulgada a aplicação do benzoato de tetra-hexilamônio (N(C6H9)4.C6H5-COO) como um tipo diferente de solvente que foi verificado como sal líquido a 25ºC. Por este ser miscível em solventes orgânicos variados e apresentar condutividade inerente, foi utilizado para estudos eletroquímicos. No começo dos anos 80 adentrou-se os líquidos iônicos obtidos pela reação do cloreto de 1,3- dialquilimidazólios com cloreto de alumínio (AlCl3), que possuem temperaturas de fusão mais baixas e maior janela eletroquímica de potencial elétrico que seus análogos contendo o cátion alquilpiridínio. (FRANZOI et al., 2011) Após estas e outras descobertas na área, surgiu uma nova classe de solventes: sais de natureza iônica (constituída somente por íons) formados por ânions inorgânicos e cátions orgânicos, caracterizados por fraca interação interiônica (possivelmente em função da dispersão de carga por esses compostos), baixa 14 energia de retículo cristalino, baixo ponto de fusão (inferiores à 100°C), baixa volatilidade (em função da sua pressão de vapor ser desprezível a temperatura ambiente), praticamente não inflamáveis e apresentando boa estabilidade térmica e eletroquímica, capacidade catalítica (posto que têm sítios capazes de se coordenarem aos catalisadores) e de solvatar em uma grande gama de moléculas polares e apolares, dentre outros atributos. (FRANZOI et al., 2011) Além disso, as mais diversas combinações de cátions e ânions podem levar a propriedades químicas específicas (termodinâmicas, de solvatação, de transporte, etc.), e com isso ampliando os campos de atuação dos líquidos iônicos (síntese catalítica, lubrificantes, processos de separação, etc.). (BORODIN, 2009) Estes foram então nomeados de sais orgânicos líquidos à temperatura ambiente, ou apenas líquidos iônicos (ILs, do inglês Ionic Liquid), como são popularmente conhecidos hoje em dia. Fundamentalmente, são eletrólitos desenvolvidos pelo contato de um cátion assimétrico volumoso com um ânion fracamente coordenante, sendo líquidos não só a temperatura ambiente, mas numa ampla faixa de temperaturas. (PAULA, F., 2006) O estado líquido alcançado por esses compostos é resultante da sua entalpia de fusão ser relativamente baixa enquanto sua entropia de fusão é elevada. (FRANZOI et al., 2011) Essa nova classe de solventes encontra-se em crescente investimento e é particularmente interessante por causar menor dano ao ambiente quando comparada com a maior parte dos solventes usados pela indústria química, além da sua baixa toxicidade e de serem potencialmente recicláveis. Com isso, ocorre minimização de geração de resíduos ao ambiente, sendo desse modo coerentes com os princípios da química limpa. (BRASIL, 2013) Além disso, os ILs podem ser “modulados” para serem aplicados a um determinado processo específico na indústria química, ou seja, há a possibilidade de ser efetuado o design de um solvente “perfeito” para o processo desejado, modificando as propriedades de interesse por meio da escolha do cátion e ânion adequados. Estima-se que podem existir da ordem de 109 líquidos iônicos possíveis. (Da SILVA, 2004; BRASIL, 2013) Os cátions e ânions mais comumente encontrados nos ILs são mais complexos do que aqueles apresentados por simples sais inorgânicos em geral (como por exemplo Na+ e Cl- que formam o cloreto de sódio). No que diz respeito aos cátions, os mais populares são os cátions dialquilimidazólio. O imidazólio trata- se de um composto aromático heterocíclico composto por um anel de cinco 15 membros contendo dois átomos de nitrogênio, onde é substituído com duas cadeias alquil diferentes. Já os ânions mais comuns a esse grupo são preferencialmente fracamente básicos, apresentam uma distribuição distorcida de carga negativa e são, em sua maioria, menores que seu cátion correspondente. (SCHRÖDER et al., 2011) 1.1.1 Tetracloroaluminato de 1-etil-3-metilimidazólio No presente trabalho foi empregado o cátion 1-etil-3-metilimidazólio, composto por um cátion imidazólio com dois grupos substituintes nas posições 1 e 3, metil e etil, de sigla EMI+ e com massa molar de 111,1650 g.mol -1. O anel aromático dispersa a carga positiva na superfície do cátion e com isso reduz a concentração da densidade de carga. (Da SILVA, 2004) Figura 1 – Estrutura molecular do cátion 1-etil-3-metilimidazólio. Fonte: Autora, 2015. Como ânion foi utilizado o tetracloroaluminato, de fórmula química [AlCl ]-4 , um composto de geometria tetraédrica, solúvel em solventes orgânicos, que pode tanto decompor-se em sais inorgânicos e haletos simples, como têm a capacidade de formação e interconversão de “clusters” com composição variada. (de ANDRADE J., 2008) Possui massa molar 168,7918 g.mol-1. A alta simetria do ânion dispersa a carga negativa sobre os quatro átomos de cloro ligados ao átomo central (Alumínio). Devido a estas características, tanto o cátion quanto o ânion não exibem grande polarização superficial de carga, com isso não ocorrendo uma forte agregação de íons. Consequentemente, não há formação de cristal, ou seja, o meio permanece líquido até temperaturas bastante baixas. Para o IL em estudo o ponto de fusão relatado na literatura é 7 ºC. (Da SILVA, 2004) 16 Figura 2 – Estrutura molecular do ânion tetracloreto de alumínio. Fonte: Autora, 2015. Uma propriedade de transporte conveniente para o estudo dos ILs e para medir o desempenho de um fluxo é o coeficiente de difusão ou autodifusão. Também chamado de difusividade de massa, corresponde a facilidade com que determinado soluto move-se em um solvente. Em um processo dinâmico, diversos elementos individuais, como momento linear e as próprias moléculas, exibem uma trajetória randômica, e como consequência desses movimentos individuais aleatórios o conjunto se difunde. Em geral, a difusão acontece enquanto o sistema é conduzido para o estado de equilíbrio. Para um processo de difusão usual, há uma dependência linear para com o crescimento temporal do deslocamento médio quadrático (MSD – do inglês, mean squared displacement). (PEDRON; MENDES, 2005) Já a viscosidade é uma propriedade dos fluídos correspondente ao transporte microscópico de quantidade de movimento em função da difusão das moléculas e em função disso é considerada uma propriedade chave para tecnologia e engenharia química, visto que avalia o desempenho de um fluxo líquido afetando tanto a transferência de calor e de massa como também o transporte, a agitação e separação da mistura. Logo, a alta viscosidade típica dos líquidos iônicos costuma ser uma desvantagem para sua utilização na indústria química, uma vez que o fluído passa a se movimentar com menor velocidade. Desse modo se faz necessário o estudo microscópico desses compostos por meio de simulações computacionais (especialmente investigando as ligações de hidrogênio e energias de interações intermoleculares), para que se possa projetar ILs com baixa viscosidade e, portanto, mais vantajosos industrialmente. Como o tamanho e a forma do composto afetam 17 essa propriedade, geralmente os líquidos iônicos que possuem baixa viscosidade são aqueles com estruturas relativamente pequenas e planares. (ZHANG et al., 2015) Objetivando a investigação de propriedades físico-químicas do ponto de vista microscópico dessa classe de compostos, introduziram-se os métodos computacionais para o estudo dos ILs. Nesses métodos, é feita uma abordagem da Química Teórica com o auxílio de computadores, os quais possibilitam o avanço de cálculos numéricos suficientemente complexos que se tornariam impossibilitados de serem realizados sem automação. (GU; YAN et al., 2013) 1.2. Dinâmica molecular e algoritmos As simulações computacionais cumprem o princípio da conservação da energia. Sendo assim, a quantidade total de energia em um sistema isolado permanece constante. A energia cinética é dada pela expressão da mecânica clássica, enquanto que a energia potencial (U) é caracterizada pelo campo de força empregado. (De ANDRADE J., 2008) O sistema simulado é representado, fundamentalmente, por uma “caixa” de simulação na qual estão contidos os átomos/moléculas de interesse. Esse sistema é muito menor que um sistema real, ou seja, a quantidade de átomos utilizados na simulação é significativamente inferior ao número de Avogadro. Como resultado, a superfície das fronteiras da “caixa” tem grande influência no resultado final da simulação: é o chamado “efeito de superfície”. Para tanto, a fim de corrigir esse efeito utilizam-se as “Condições Periódicas de Contorno” (PBC, do inglês “Periodic Boundary Conditions”). Trata-se de uma abstração, em que se impõe que o sistema encontra-se ladeado de réplicas idênticas de si mesmo, infinitamente. Desse modo, cada partícula presente na caixa de simulação irá possuir imagens em cada uma de suas réplicas e em todas as direções. Por um lado, assim assegura-se que se determinada partícula sair da caixa, por qualquer que seja dos lados, sua imagem entrará pelo lado oposto, preservando o número total de partículas presentes, como representa a figura 3. Por outro lado, esse fato exige a aplicação da correção chamada “Convenção de Imagem Mínima” (do inglês, “Minimum Image Convention”), que emprega um raio de corte a fim de excluir as interações artificiais de uma molécula com suas “réplicas”. (NEYERTZ, et al. 1996) 18 Figura 3 – Representação visual das Condições de Contorno Periódicas. Fonte: (ADRIANE; SEGATTO, 2013) A partir das posições e velocidades iniciais das partículas presentes no meio, é possível obter as coordenadas dessas partículas em um tempo posterior. A resolução de cálculos como esse só é possível por aproximações de derivadas por diferenças finitas, é o chamado método de diferenças finitas, no qual ocorre a integração das equações diferenciais que descrevem o movimento. Para esse estudo a forma de integração utilizada se deu através do algoritmo Verlet. (NEYERTZ, et al. 1996) Nesse algoritmo, a força resultante (Fi) sobre uma partícula (i) em um instante de tempo (t) é dada por: (NEYERTZ, et al. 1996) d ²ri Fi (t)ai   ri (t)  dt² mi 19 Onde, ai representa a aceleração da partícula, ri a posição de acordo com o tempo e mi a massa. Para obter as posições calculadas no tempo t  t , se faz necessário a utilização de uma expansão de Taylor: (NEYERTZ, et al. 1996)  t²  t³ ri (t  t)  ri (t)  tri (t)  ri (t)  r 4 i (t)  n( t ) 2! 3!  t²  t³ ri (t  t)  ri (t)  tri (t)  ri (t) ri (t)  n( t 4) 2! 3! Somando-se ambas as expressões, resulta em: (NEYERTZ, et al. 1996)  t² r 4i (t  t)   ri (t  t) 2ri (t)  Fi (t) ( t ) m Com a qual são calculadas as posições das partículas em t  t . Já para calcular a velocidade da partícula é necessário subtrair as equações, originando a expressão: (NEYERTZ, et al. 1996) 1 ri (t)  [ri (t  t)  ri (t  t)] 2 t 1.3. Potenciais e campos de força Os métodos de cálculo utilizados na dinâmica molecular ou MD (do inglês, Molecular Dynamics) clássica baseiam-se num aspecto clássico da estrutura molecular como um conjunto de esferas ligadas por molas que possuem constantes de força características. O campo de força pode ser descrito como uma função matemática que possui parâmetros ajustáveis e que objetiva descrever a energia potencial de um determinado sistema ou conjunto de partículas, sendo então constituído do conjunto de parâmetros que descrevem as interações componentes da função potencial do sistema. Através dos campos de força, se faz possível descrever vários fenômenos moleculares, como mudanças de conformação e interações, de diversas classes de moléculas. A confiabilidade do campo de força varia com a qualidade dos parâmetros que são incorporados na função do potencial energético. Estes são divididos em duas categorias: os campos não polarizáveis ou modelos de íons rígidos (RIMs, do inglês, rigid ion model) e modelos de íons polarizáveis (PIMs, do inglês polarizable ion model). A inserção de componentes 20 energéticos polarizáveis mostra-se auxiliar para uma dinâmica coletiva apropriada, de sistemas fortemente carregados, corrigindo o desvio negativo nos coeficientes de difusão, uma vez que neste modelo os íons exibem uma maior mobilidade na simulação computacional. (SCHRÖDER et al., 2011) O campo de força AMBER (do inglês, Assisted Model Building with Energy Refinement) é a reunião de parâmetros que descrevem as interações presentes entre as espécies e a energia potencial intramolecular. Na última década, observou- se uma crescente diversidade de modelos clássicos de dinâmica molecular em líquidos iônicos, baseados em campos de força com cargas permanentes. (TAYLOR et al., 2013) Para o campo de força AMBER, a energia potencial é descrita pela seguinte equação: (BARREIRO et al., 1997) Utotal Uligação U U UÂngulos diedros nãoligados Nesta, os dois primeiros elementos (U e U ) correspondem ligação Ângulos aos estiramentos e deformações angulares respectivamente, descritas com potenciais harmônicos. Já o termo Udiedros refere-se ao potencial de torção dos ângulos diedros próprios e impróprios da molécula, e é representado por um potencial periódico. O último termo, U , compreende às interações de van der nãoligados Waals (modeladas de acordo com o potencial de Lennard-Jones) e as interações eletrostáticas (regidas pelo potencial de Coulomb entre as cargas atômicas pontuais). (de ANDRADE J. et al.,2008) Principalmente no que diz respeito aos ILs, as interações eletrostáticas são de extrema importância, uma vez que o sistema se encontra altamente carregado com cátions e ânions, e as interações entre eles interferem diretamente na organização do sistema. Sendo assim, torna-se imprescindível o cálculo apropriado de cargas pontuais para o resultado final. No presente estudo foi utilizado o método RESP para o cálculo de cargas pontuais. Esse método baseia-se em um mapa de potencial eletrostático que circunda as moléculas. Nele, restringe-se as cargas para aproximadamente uma carga atômica ideal (zero) de acordo com sua magnitude relativa. As cargas parciais são oriundas do ajuste à esse mapa de potencial eletrostático que é calculado ao redor do arranjo tridimensional dos átomos em uma 21 rede de pontos no espaço, sendo também capaz de reproduzir de maneira adequada diversos momentos multipolares das moléculas. (BOURSCHEIDT et al. , 2005) 1.3.1 Polarização Assim, em campos de força (do tipo RIM), cada átomo é imaginado particularmente como uma esfera rígida, que interage com os outros átomos presentes no meio conforme as leis da mecânica Newtoniana, e a cada átomo atribui-se apenas uma única carga elétrica que permanece inalterada ao longo da simulação. Desse modo tornam-se limitados os efeitos eletrônicos no meio, bem como impossibilitando-se a presença de dipolos induzidos. Para íons pequenos, como o AlCl -4 , a polarização iônica é, basicamente, aditiva. Campos de força do tipo PIM baseiam-se no modelo clássico de Drude, entre outros (PPD e FQ). Esse modelo simula os efeitos da polarizabilidade eletrônica no campo de força por meio da inserção de partículas Drude. (ANISIMOV et al., 2007) Desse modo, para obter-se respostas mais realistas aos átomos da simulação, incluem-se partículas de Drude: tratam-se de sítios virtuais, que não apresentam massa, possuem apenas uma carga parcial e acoplam-se aos átomos por meio de uma mola harmônica. Sendo assim, tanto a constante da mola, quanto as cargas parciais sobre os átomos e partículas Drude, determinam a resposta para o campo eletrostático local aplicado. Logo, as partículas Drude servem como um mediador para alterações na distribuição de carga eletrônica do sistema, possibilitando a polarização induzida. (ANISIMOV et al., 2007) Conforme o modelo Drude, cada átomo é apresentado por duas cargas pontuais, qcore e  , que estão ligadas através de uma mola harmônica de constante de força KD . Quando somadas, essas duas cargas pontuais geram uma carga atômica parcial associada ( qA ) a determinado átomo A. (ANISIMOV et al., 2007) qA  qcore  Quando se aplica um campo elétrico externo (E), ocorre um deslocamento (Δx) da partícula Drude que está carregada, que por sua vez dá origem a um dipolo induzido (  ).(ANISIMOV et al., 2007)  E 22 O deslocamento (Δx) é oposto, já que a mola é uma força restauradora: F harm  KD.x E, finalmente, a polarizabilidade de um átomo é dada por: (ANISIMOV et al., 2007)  ²   KD 1.4. Propriedades de transporte 1.4.1 Difusão A difusão pode ser definida como o processo pelo qual um conjunto de partículas é transportada em função de um gradiente de concentração. Quando essas partículas estão sujeitas a um campo, seja gravitacional ou elétrico, além do fluxo difusivo, as partículas também experimentam um movimento coletivo, em uma mesma direção guiada pelo campo de força. A esse movimento dá-se o nome de arraste ou arrasto. (WEISS, 1869) Um processo usual de difusão (comportamento difusivo), de maneira geral, tem como característica básica o deslocamento quadrático médio ou MSD crescente linearmente com o tempo. Porém, experimentos constataram desvios na difusão normal, que poderia vir a se tornar mais lenta ou mais rápida, o que foi denominado de “difusão anômala”. A escala das trajetórias de difusão pode ser, de maneira simplificada, classificada da seguinte forma:(ASSIS, 2006)  r²(t)  t (1) Na qual  r²  representa o movimento aleatório das partículas na dimensões x,y e z ao longo do tempo de simulação. A difusão é então classificada pelo índice  : quando  1 trata-se de uma difusão “normal”, já caso  1 temos um comportamento anômalo “superdifusivo”, (incluindo a hipótese particular na qual   2 : nesse caso, é chamado comportamento “balístico”). Para  1 são os processos chamados de “subdifusivos”. Sendo assim, em um gráfico de log  r²  vs log t torna-se possível 23 determinar a classificação do comportamento difusivo presente no sistema em estudo. A fim de calcular o coeficiente de autodifusão das moléculas do solvente, é possível utilizar a análise do MSD, na qual obtém-se o coeficiente angular (m) da reta por regressão na região do gráfico onde é observado o comportamento difusivo. Dado esse coeficiente angular do MSD, o coeficiente de autodifusão da molécula (D) é dado por: m D  (2) 6 Outra forma de calcular o coeficiente de autodifusão surge do fato que a dinâmica molecular clássica baseia-se na evolução de coordenadas (r) e momento (p) de um sistema ao longo do tempo, quando aplicadas leis clássicas de movimento sobre uma partícula (i). Sendo assim, a função das coordenadas no espaço de fase pode ser genericamente definida por: (TAVARES; COSTA, 2007) A{p(t),r(t)} A(p,r;t)  A(t) Sendo A(t) qualquer propriedade dependente do tempo. Então, a função de correlação no tempo de A(t) pode ser descrita como: (TAVARES; COSTA, 2007) N CA(t ')   Ai (t ')Ai (0) i1 Onde Ai (t ') indica o valor da propriedade em questão em determinado tempo e Ai (0) o valor de A em diversas origens de tempo. Assim, a função pode por exemplo ser reescrita em termos de velocidade de uma determinada partícula i presente no sistema: (TAVARES; COSTA, 2007) CV (t ')  Vi (t').Vi (0) E para esse caso específico, CV (t ') é dita função de autocorrelação da velocidade, ou VACF (do inglês, Velocity Autocorrelation Function), das moléculas do sistema. A partir dessa função, torna-se possível calcular o coeficiente de 24 autodifusão através da relação de Green-Kubo para a difusão: (TAVARES; COSTA, 2007)  1  1 D   Vi (t ').Vi (0)dt '   CV (t ')dt ' (3) 3 0 3 0 1.4.2 Viscosidade A viscosidade de um determinado fluxo representa a resistência do mesmo ao escoamento, ou deformação progressiva. Essa propriedade é resultado das interações entre partículas vizinhas que estão em movimento no meio, a diferentes velocidades e direção. A viscosidade nos líquidos iônicos, é influenciada pelas interações eletrostáticas e, em sua maioria, pela tendência de formação de ligações de hidrogênio, bem como pelas interações de van de Walls. Quanto maior for a intensidade das ligações de hidrogênio presentes, mais fortemente as espécies estarão atraídas entre si. Dessa forma eleva-se a dificuldade de difusão das partículas ao longo do líquido e com isso aumentará a viscosidade do IL. (CONSORTI et al., 2001) No que diz respeito aos ânions, sua influência na viscosidade se dá por meio da interação dessas moléculas com os cátions, também através de ligações de hidrogênio e/ou interações de van der Walls. Dessa forma, quanto maior for a distribuição de cargas, mais fracas serão as ligações de hidrogênio formadas, o que acarretará em uma menor viscosidade do sistema. Desse modo, o ânion em estudo, AlCl -4 , devido a sua alta simetria dispersa a carga negativa sobre os 4 átomos de cloro, o que acarreta em uma maior dispersão das cargas sobre a superfície da molécula, com isso influenciando na redução da viscosidade do IL. Em relação ao cátion EMI+, trata-se de um dos cátions com menor viscosidade já estudados na literatura, pois possui cadeia alquílica curta que, associada com uma relativa baixa massa molar, promove aos ILs formados por esse cátion uma mobilidade satisfatória, gerando uma baixa viscosidade. Propriedades individuais das partículas, como por exemplo o coeficiente de difusão, podem ser calculadas com mais precisão, por MD, do que propriedades 25 coletivas, como a viscosidade, uma vez que no primeiro caso, a média estatística é realizada para as partículas ao longo do tempo, enquanto que para a viscosidade perde-se a possibilidade de média sobre as partículas, reduzindo a precisão estatística. Para aumentar a precisão, se faz necessário corridas mais longas para propriedades coletivas. (MEDINA et al., 2011) O coeficiente de viscosidade pode ser calculado por variados métodos partindo de simulações. Aqueles mais simples baseiam-se em MD de equilíbrio e empregam como propriedade dependente do tempo as flutuações (dos termos fora da diagonal do tensor) de pressão, Pxz. Assim, a relação de Green-Kubo para a viscosidade de cisalhamento é dada por: 1     Pxz (t ') P xz (0)  dt ' (4) VkBT 0 Onde  é a viscosidade, V o volume e kB a constante de Boltzmann. Essa pode ainda ser reescrita como relação de Einstein: 1V d t '   lim ( Pxz (t ')dt ')² t (5) t ' 2k T dt ' 0 0B A viscosidade ainda pode ser determinada utilizando a TCAFs (do inglês Transverse-current Correlation Functions) A TCAF é calculada no equilíbrio e consequentemente requer simulações relativamente longas. Em comparação com o método mais comumente utilizado para esse tipo de cálculo, a relação de Green-Kubo, a TCAF é mais vantajosa por identificar facilmente um comportamento não hidrodinâmico (influencia na transferência de massa/calor pelo sistema) e fornece um método natural para estimar a magnitude dos chamados efeitos de tamanho finito (decorrentes do confinamento geométrico dos átomos em um volume fixo). (PALMER, 1994) Em escala microscópica, o fluído é constituído de partículas discretas que executam movimentos térmicos aleatórios mesmo quando no equilíbrio. Esses 26 movimentos acarretam na aparição de gradientes de momento que podem então ser analisados para obtenção de coeficientes de transporte. (PALMER, 1994) Para o cálculo da TCAF utiliza-se as equações dos coeficientes da expansão de Fourier para o campo microscópico de momento transversal, na forma: N  micro (k, t)  kˆ j .p (t)sin[k.r j (t)] (6) j1 N  micro (k, t)  kˆ j j .p (t) cos[k.r (t)] (7) j1 Onde  micro são os coeficientes de Fourier, k é o vetor no espaço de Fourier, t  é o tempo e p é o momento microscópico. A função de correlação dos diferentes  micro é que é denominada TCAF. Para sistemas infinitos, a TCAF deve depender  somente da magnitude do vetor k. (PALMER, 1994) Calculam-se então os coeficientes de Fourier. Mas, devido a uma inversão de simetria do sistema, o  (k) obtido é uma função quadrática de k, e assim: (k)   ak² (8) Onde  é o limite infinito do sistema (k nulo), cujo valor extrapolado é o resultado procurado. 27 2. METODOLOGIA A MD propõe o cálculo do movimento das moléculas em um determinado sistema por meio do potencial de interação entre suas partículas e das equações que conduzem seu movimento. Desse modo, torna-se possível a análise da evolução das posições relativas e energia de cada átomo ao longo do tempo. Assim, é possível obter-se as propriedades macroscópicas desse sistema, como os coeficiente de difusão e viscosidade, pois o conjunto de coordenadas de cada átomo ao decorrer da evolução temporal caracteriza a trajetória do mesmo. (ADRIANE; SEGATTO, 2013) Os cálculos computacionais foram executados no sistema operacional Linux, Kernel versão 3.1.10-1.19, distribuição OpenSUSE, versão 12.1. A geração da caixa de simulação foi feita através do programa Packmol, versão 1.1 (MARTÍNEZ et al., 2009), e para as simulações de dinâmica molecular e análise das trajetórias foi utilizado o pacote de programas GROMACS, versão 5.0.2 (PRONK, et al., 2013). (do inglês GROningen MAchine for Chemical Simulations). 2.1 Geração da configuração inicial das moléculas. Dada a massa molar do líquido iônico em estudo, de 279,9568 g/mol, calculou-se o lado da caixa cúbica de simulação, utilizando 1000 moléculas, como sendo 59,93 Å (para a densidade encontrada na literatura de 1,302 g/cm³). Em virtude das grandes intensidades das interações entre os íons que compõem os ILs, não se iniciou a etapa de equilibração na densidade real de 1,302 g/cm³, mas a partir de uma densidade mais baixa. No caso utilizou-se 0,500 g/cm³, que foi aumentada continuamente na sequência. Caso contrário o sistema “explodiria” e aos átomos passariam a não mais interagir. A geração da “caixa” de simulação inicial foi efetuada com o uso do programa Packmol. (MARTÍNEZ et al., 2009) 2.2 Equilibração do sistema. A seguir, procedeu-se a equilibração do sistema. Esta se deu em duas fases. A primeira fase, objetivando a adequação da densidade do sistema, se deu através 28 de nove sucessivas subetapas de minimização da energia e dinâmica molecular, que foram procedidas de acordo com o algoritmo representado pela figura 4: Figura 4 – Fluxograma executado no ajuste da densidade do sistema. Fonte: Autora, 2015. O procedimento se seguiu sucessivamente até a configuração atingir a densidade objetivada de 1,302 g/cm³. De acordo com a representação do fluxograma, a densidade foi sendo acrescida de 0,100 g/cm³, como pode ser observado na figura 5: 29 Figura 5 - Variação da densidade do sistema. Fonte: Autora, 2015. As subetapas de minimização energética, também chamadas de otimização geométrica, ocorreram imediatamente após cada etapa de mudança de densidade do sistema e objetivaram localizar um conjunto de coordenadas que relaxassem adequadamente a energia potencial do sistema, ajustando as posições dos átomos. Assim, relaxando a distorções nas ligações químicas, nos ângulos entre ligações e nas interações de van der Walls, alcança-se uma estrutura de partida na qual as simulações de dinâmica molecular podem ser efetuadas sem que surjam por exemplo estiramentos irreais de ligações químicas (“explosões” no sistema simulado) por consequência da integração numérica das equações de movimento. O algoritmo empregado foi Steepest Descend, com 200 passos no máximo, parâmetros de tolerância de energia 500 e de tamanho de passos 0,0001, ambos em unidades internas. (NAMBA et al., 2008) A cada subetapa de minimização energética seguiu-se uma subetapa de dinâmica molecular. Cada subetapa de MD do processo de ajuste da densidade do sistema foi realizada com 10000 passos temporais de integração de 0,002 ps. A correção do centro de massa foi efetuada a cada passo. Essas subetapas foram executadas sem barostato e empregando termostato Nosé-Hoover. 30 Uma vez ajustada a densidade, a segunda fase de equilibração do sistema constituiu-se de subetapas sucessivas de longas dinâmicas, com pequenas perturbações entre cada uma. Em cada uma utilizaram-se 5 milhões de passos temporais de 0,002 ps. Foram verificadas a estabilização termodinâmica do sistema por diversos parâmetros: como energia total, energia potencial, força coulômbica, ângulos diedros, densidade, tamanho da caixa, volume e produto pressão-volume. Entre estas subetapas de equilibração foi modificado o barostato utilizado, como indicado na tabela 1: Tabela 1 – Barostato, termostato e número de passos empregados nas 3 subetapas da segunda fase de equilibração do sistema. Subetapa Termostato Barostato Número de passos 10 Nosé-Hoover desativado 5.000.000 11 Nosé-Hoover Berendsen 5.000.000 12 Nosé-Hoover Parrinello-Rahman 5.000.000 Fonte: Autora, 2015. O termostato utilizado, Nosé-Hoover, controla o período das flutuações de temperatura no estado de equilíbrio. Os parâmetros empregados para o mesmo foram temperatura constante de 298 K. Na primeira subetapa, o barostato encontra- se desativado, não havendo acoplamento de pressão e sendo assim o tamanho da caixa permanece fixo. Na segunda subetapa, utilizou-se o barostato Berendsen (isotrópico – mantem as mesmas propriedades físicas independente da direção considerada) com parâmetros de tempo de acoplamento 2ps e compressibilidade 4,5x10-5 bar-1, onde o relaxamento do acoplamento de pressão é exponencial ao tempo e a caixa é redimensionada a cada passo. Tem sido argumentado que essa opção não produz um conjunto termodinâmico correto, porém mostra-se a maneira mais eficiente de dimensionar a caixa de simulação no início de uma corrida. Por último, foi empregado o barostato Parrinello-Rahman com parâmetros de tempo de acoplamento 2ps e compressibilidade 4,5x10-5 bar-1, no qual os vetores da caixa estão sujeitos a uma equação de movimento, acoplada à equação de movimento dos átomos. Finalmente, chegou-se a configuração equilibrada para o sistema em questão, empregando um barostato e termostato adequados para o ensemble NpT – 31 isotérmico-isobárico (onde o número de partículas, a temperatura e a pressão são mantidos constantes). 2.3. Perturbando o campo de força: modelos de transferência estática de cargas Nos ILs, a transferência de carga é proveniente de uma deslocalização eletrônica entre os orbitais do cátion e do ânion que compõem o sistema. A introdução dessa característica no modelo, ainda que de maneira estática e não dinâmica, aprimoraria a modelagem do IL e teria como consequência a redução da excessiva forte contribuição artificial das interações coulômbicas na energia potencial do sistema. Dessa forma, espera-se observar melhoria nas propriedades de transporte simuladas, que são sabidamente inibidas pela intensidade das interações coulômbicas modeladas, aproximando-se as encontradas experimentalmente. (WEISSHEIMER, 2013) Para a modificação de cargas, se fez necessário o cálculo das mesmas (previamente calculadas em outros estudos, no nível de teoria MP2//6-311(G(d)) para o par iônico. Esse é o doravante chamado “modelo 2”, que se diferencia do modelo original (doravante, “modelo 0”, cujas cargas foram calculadas empregando o nível de teoria RHF//6-31G(d)) justamente por ser baseado nos cálculos efetuados em pares iônicos e não em íons isolados. Numa proposta alternativa, ainda se empregou o artifício de multiplicar a carga resultante do cátion pareado (0,790623) pelas cargas parciais já existentes no modelo 0, com o intuito de modelar o efeito da transferência de cargas sem redistribuição. Este será doravante chamado de “modelo 1”. As cargas pontuais empregadas podem ser associadas a cada átomo correspondente de acordo com a numeração apresentada na figura 6. Essa mesma numeração é empregada na tabela 2, que apresenta para cada átomo sua respectiva carga parcial nos três modelos testados. 32 Figura 6 - Número associado a cada átomo da estrutura molecular. Fonte: Autora, 2015. 33 Tabela 2 – Distribuição de carga parcial nos átomos. Átomo Modelo 0 Modelo 1 Modelo 2 1 0,032818 0,025947 0,017786 2 - 0,165661 - 0,130975 - 0,122775 3 0,058110 0,045943 0,178978 4 0,079453 0,062817 0,261895 5 - 0,079455 - 0,062819 - 0,207493 6 0,258666 0,204507 0,218216 7 - 0,167115 - 0,132125 - 0,155297 8 0,248174 0,196212 0,203363 Anel imidazólio 0,264990 0,209507 0,394673 9 - 0,192460 - 0,152163 - 0,354022 10 0,129251 0,102189 0,158688 11 0,129251 0,102189 0,158688 12 0,129251 0,102189 0,158688 Metila 0,195293 0,154403 0,122042 13 0,204594 0,161757 0,172050 14 0,089159 0,070491 0,041555 15 0,089159 0,070491 0,041555 16 - 0,010004 - 0,007909 - 0,205976 17 0,055603 0,043961 0,074908 18 0,055603 0,043961 0,074908 19 0,055603 0,043961 0,074908 Etila 0,539717 0,426717 0,273908 Cátion +1,00000 +0,790624 +0,790623 20 0,909248 0,718872 0,726569 21 - 0,477312 - 0,377374 - 0,379298 22 - 0,477312 - 0,377374 - 0,379298 23 - 0,477312 - 0,377374 - 0,379298 24 - 0,477312 - 0,377374 - 0,379298 Ânion -1,00000 -0,790624 -0,790623 Fonte: Autora, 2015. 34 2.4 Aquisição de trajetórias e análises Esta etapa novamente foi dividida em três fases. Na primeira fase, exploratória, foram efetuados 5 milhões de passos de simulação, salvos a cada 5 passos. Essa amostragem se mostrou inadequada, tanto para a análise das VACFs (que claramente exigiram uma maior frequência de gravação de configurações), quanto para o MSD, a relação de Einstein e as TCAFs (que por sua vez apresentaram a necessidade de trajetórias mais longas). Para corrigir essas dificuldades, foram efetuadas duas fases adicionais de aquisição de trajetórias. Na segunda fase, foram efetuadas simulações de apenas 50 000 passos, salvos a cada passo. Essas simulações permitiram a obtenção dos resultados adequados para as VACFs que são apresentados na próxima seção. Finalmente, na terceira e última fase, foi efetuada uma simulação de 25 milhões de passos, salvos a cada 500. Essas trajetórias mostraram-se adequadas para as outras análises, cujos resultados são também apresentados na seção de resultados e discussão. 35 3. RESULTADOS E DISCUSSÃO 3.1 Coeficiente de autodifusão. A fim de encontrar o coeficiente de autodifusão para os íons constituintes do sistema em estudo, foi inicialmente calculado o MSD dos íons, o qual fornece uma maneira simples de calcular o coeficiente de autodifusão. Esse coeficiente de autodifusão é calculado com base um ajuste de mínimos quadrados em uma linha reta traçada pelo MSD a tempos longos, aplicado depois na equação 2. Os resultados desse processo para o cátion e o ânion podem ser verificados nas figuras 7a e 7b sem escala logarítmica, respectivamente, e nas figuras 8a e 8b, em escala logarítmica. 36 Figura 7 – MSD (a) do cátion EMI+ e (b) do ânion AlCl -4 nos 3 modelos. (a) (a) (b) Fonte: Autora,2015. 37 Figura 8 – MSD (a) do cátion EMI+ e (b) do ânion AlCl -4 em escala logarítmica para os 3 modelos (a) (b) Fonte: Autora,2015. Nas figuras 8a e 8b observa-se que a já sabida falha de campos de força do tipo RIM em modelar a difusão de ILs é mantida mesmo com a inclusão de um modelo de transferência de cargas. Isso é visível pois nessas figuras os íons deveriam inicialmente possuir um comportamento dinâmico balístico. Isso se daria pois, uma vez que as partículas movem-se com velocidade (v) (supondo que todas as partículas sob ação do campo de força se movem com a mesma velocidade, sem sofrer colisões nem experimentar força de atrito) constante ao longo do tempo, temos que r  v.t : consequentemente, segundo a equação 1:  r²(t) ~ t , onde   2 . Porém, o que se observa é um comportamento inicial aproximadamente 38 difusivo com o deslocamento quadrático médio avançando linearmente ao tempo. (Para o cátion: modelo 0 -  = 0,930, modelo 1 - = 1,059, modelo 2 - = 1,021; para o ânion: modelo 0 -  = 1,097, modelo 1 -  = 1,204, modelo 2 -  = 1,166). Esse comportamento usual difusivo deveria então manter-se ao longo da simulação molecular, entretanto, percebe-se que as simulações passam a se comportar de maneira subdifusiva (Para o ânion: modelo 0 -  = 0,373, modelo 1 -  = 0,507, modelo 2 -  = 0,391; para o cátion: modelo 0 -  = 0,393, modelo 1 -  = 0,686, modelo 2 -  = 0,533). Esses dados são apresentados nas Tabelas 3 e 4. Tabela 3 – Resultados obtidos das análises de MSD do cátion. MSD Região Balística Região Difusiva Modelo Ti Th  Ti Th  D (cm2/s) Ti Th  D (cm2/s) 0 0 0,4ps 0,930 200ps 600ps 0,397 7,448x10-8 10ns 15ns 0,727 2,813 x10-8 1 0 0,4ps 1,059 200ps 600ps 0,686 7,237x10-8 10ns 15ns 0,978 5,988 x10-8 2 0 0,4ps 1,021 200ps 600ps 0,533 2,142x10-7 10ns 15ns 0,837 1,105 x10-7 Esperado: 2 Esperado: 1 8,170x10-7 Esperado: 1 8,170 x10-7 Fonte: Autora, 2015. 39 Tabela 4 – Resultados obtidos das análises de MSD do ânion. MSD Região Balística Região Difusiva Modelo T 2 2i Th  Ti Th  D (cm /s) Ti Th  D (cm /s) 0 0 0,4ps 1,097 200ps 600ps 0,373 1,523x10-7 10ns 15ns 0,585 1,795x10-8 1 0 0,4ps 1,204 200ps 600ps 0,507 4,797x10-7 10ns 15ns 1,010 4,853x10-8 2 0 0,4ps 1,166 200ps 600ps 0,391 2,152x10-7 10ns 15ns 0,890 1,289x10-7 Esperado: 2 Esperado: 1 Esperado: 1 Fonte: Autora, 2015. Já no que diz respeito às regiões entre 200ps e 600ps para ambos os íons, os resultados no comportamento difusivo anômalo obtidos pelo modelo 0 aproximam-se daquele relatado na literatura em estudos prévios (para o intervalo de 100ps à 150ps). Além disso, o valor de  melhora e se aproxima do esperado nos modelos incorporando a transferência estática de cargas (em especial no modelo 1). Mais importante ainda, pode ser observado que o comportamento difusivo dos íons fica muito próximo ao esperado ( =1) quando analisadas dinâmicas de tempo longo (entre 10 ns e 15 ns). Esse fato é especialmente verdade em relação ao modelo 1, que chegou a ser de 0,978 para o cátion e 1,010 para o ânion. A análise de tempos longos também mostrou melhorias no modelo 0 nesse sentido, ainda que muito inferiores aos modelos 1 e 2. Retornando para as figuras 7a e 7b, as mesmas foram empregadas para o cálculo do coeficiente de difusão pelo uso da equação 2, usando o intervalo de tempo longo (entre 10 ns e 15 ns). A regressão linear permite inferir os seguintes coeficientes de difusão: 2,813x10-8 cm2/s, 5,988x10-8 cm2/s e 1,105x10-7 cm2/s para os cátions, e 1,795x10-8 cm2/s, 4,853x10-8 cm2/s e 1,289x10-7 cm2/s para os ânions, ambos em relação aos modelos 0, 1 e 2 respectivamente. Com isso é possível concluir que, comparado ao valor experimental de 8,170x10-7 cm2/s para o cátion, os resultados obtidos no estudo ainda são inferiores em aproximadamente uma ordem 40 de grandeza. Assim sendo, a transferência estática de carga melhorou os coeficientes de autodifusão de acordo com os coeficientes constatados, mas ainda de forma insuficiente. Este comportamento irregular a nível microscópico é, em sua maioria, provavelmente causado tanto por forças externas ao meio como por forças provenientes de outras partículas especialmente, as chamadas interações coulômbicas (interações de longo alcance entre partículas carregadas). Além do que, no momento em que ocorre um processo difusivo, diversos elementos poderão alterar-se frequentemente, como energia e momento linear. Ainda assim, é possível inferir que, de acordo com o comportamento apresentado nas figuras 7 e 8 e nas tabelas 3 e 4 para cada modelo, os ajustes de cargas reduziram a forte contribuição das interações coulômbicas no sistema. Desse modo, obtém-se comportamentos difusivos mais próximos aqueles que deveriam ser observados (no início balístico e em seguida difusivo). Observa-se ainda que apenas a transferência de cargas (modelo 1) é ainda mais eficiente que o rearranjo dessas (modelo 2). Isso pode ser interpretado em termos de dados na tabela 2: em relação ao anel imidazólio, no qual estão contidos os átomos de hidrogênio disponíveis para formação de ligações de hidrogênio, a carga positiva desse é menor para o modelo 1, de tal forma que reduz para esse modelo a possibilidade de interações do tipo ligação de hidrogênio, aumentando a mobilidade dos íons comparativamente aos modelos 0 e 2 simultanêamente. Uma forma alternativa de calcular do coeficiente de difusão é através das VACFs das moléculas, como já apresentado na equação 3. As funções obtidas podem ser vistas nas figuras 9 e 10. 41 Figura 9 – VACFs para os 3 modelos estudados do cátion EMI+, apresentadas no intervalo de tempo de correlação (a) de 0 à 2,5 ps e (b) de 2 à 25 ps. (a) (b) Fonte: Autora,2015. 42 Figura 10 – VACFs para os 3 modelos testados do ânion AlCl -4 , apresentadas no intervalo de tempo de correlação (a) de 0 à 2,5 ps e (b) de 2 à 25 ps. (a) (b) Fonte: Autora,2015. Nas figuras 9a e 10a, é possível observar que a função de autocorrelação da velocidade decai rapidamente a próximo de zero (em aproximadamente 2 ps). Esse comportamento ocorre, pois no IL um dos íons é muito leve em comparação ao outro (nesse caso, o ânion AlCl -4 é muito mais leve quando comparado ao cátion EMI+) e com isso ocorre a realização de movimentos oscilatórios em uma dinâmica de tempo curto. (TAVARES; COSTA, 2007) 43 Para o cátion (figura 9a) esse padrão oscilatório de curto tempo é visivelmente mais acentuado que para o ânion (figura 10a), dado que estes experimentam uma dinâmica mais impedida/presa. Durante a simulação, os íons são bloqueados, temporariamente, pela sua esfera de coordenação e esse bloqueio é mais efetivo nos cátions, o que torna o período de oscilação para o ânion maior e causa essa curva mais “relaxada” como observado na figura 10. A ampliação das figuras 9a e 10a (figuras 9b e 10b respectivamente) para tempos longos (2 à 25 ps) apresenta uma série de ruídos (da ordem de ±0,004 m²/s²). Assim, a função de autocorrelação da velocidade de ambos os íons decai assimptoticamente a zero, a menos de desvios numéricos e oscilações estatísticas. Infelizmente, tais flutuações acabaram comprometendo a integração da área do gráfico além do aceitável, e consequentemente o cálculo do coeficiente de difusão por esse método foi inviabilizado. 3.2 Coeficiente de viscosidade. Novamente, dois métodos diferentes foram utilizados a fim de analisar a viscosidade do líquido iônico EMImAlCl4. No primeiro é empregada a relação de Einstein, dada pela equação 5, essa deveria convergir, em dinâmica de tempo longo, para o valor da viscosidade. Os resultados obtidos por esse método são apresentados na figura 11 para os 3 modelos em estudo. 44 Figura 11 – Relação de Einstein para viscosidade calculada a partir das flutuações dos termos fora da diagonal do tensor de pressão nos 3 modelos estudados. Fonte: Autora, 2015. É possível constatar que este método utilizado não foi adequado: claramente, as flutuações no tensor de pressão não convergiram ao decorrer do tempo empregado, tornando a análise inviável e com oscilações. O segundo método calcula as TCAFs, utilizadas para estimar o coeficiente de viscosidade de cisalhamento (η). A dificuldade para o cálculo de importantes propriedades de transporte, como a viscosidade, por meio da dinâmica molecular se dá, principalmente, em função das grandes incertezas estatísticas e a complexidade das funções de autocorrelação, com a equação 5 recém empregada. Para minimizar essas e outras dificuldades, e assim calcular a viscosidade de cisalhamento, foi utilizado a TCAF. O cálculo foi realizado em séries com 5 000 ps cada para os 3 modelos testados, além da trajetória adquirida completa. A viscosidade de cisalhamento é proveniente da curva TCAF que é plotada em função de k. A dispersão observada nos dados é possivelmente uma consequência das viscosidades baixas encontradas e dos longos tempos de relaxação para o TCAF. Como resultado, há um número relativamente baixo de intervalos não correlacionados para cada simulação (alguns desses intervalos mais discrepantes foram excluídos no momento da plotagem do 45 gráfico – relação dos valores, incluídos e removidos, no Anexo). O resultado desse procedimento pode ser verificado na figura 12. 46 Figura 12 – TCAFs global e seccionados para (a) o modelo 0, (b) o modelo 1 e (c) modelo 2. (a) (b) (c) Fonte: Autora, 2015. 47 A partir das curvas ajustadas, se fez possível verificar os parâmetros  , a e o  r² provenientes do polinômio de ordem 2, para cada um dos modelos, os valores são apresentados nas tabelas 5, 6 e 7. Os elevados valores do parâmetro r² obtidos para cada seção, bem como para a trajetória global, evidenciam a alta confiabilidade dos dados obtidos, assim como a grande concordância e convergência entre os dados de  . Tabela 5 – Parâmetros do fit das TCAFs global e seccionadas do Modelo 0.  a r²  Série 1 4,89353±2,23982 x10-1 -1,03775±1,05338 x10-1 0,87394 Série 2 4,97126±2,06196 x10-1 -1,03775±9,72665 x10-2 0,89247 Série 3 4,90006±2,12154 x10-1 -1,01914±9,97751 x10-2 0,88169 Série 4 5,01663±2,33119 x10-1 -1,06789±1,10156 x10-1 0,87034 Completo 4,95154±2,16594 x10-1 -1,04685±1,02347 x10-1 0,88198 Fonte: Autora, 2015. Tabela 6 – Parâmetros do fit das TCAFs global e seccionadas do Modelo 1.  a r²  Série 1 2,80296±9,73389 x10-2 -6,29043x10-1±5,63104x10-2 0,91899 Série 2 2,76988±8,46416 x10-2 -6,07837x10-1±4,87203x10-2 0,93399 Série 3 2,75013±8,75092 x10-2 -5,98382x10-1±5,04818x10-2 0,92739 Série 4 2,76232±8,63020 x10-2 -6,09270x10-1±4,97431x10-2 0,93169 Série 5 2,73088±1,02721 x10-1 -5,92642x10-1±5,73848x10-2 0,91428 Completo 2,77208±8,52263 x10-2 -6,12222x10-1±4,91596x10-2 0,93378 Fonte: Autora, 2015 48 Tabela 7 – Parâmetros do fit das TCAFs global e seccionadas do Modelo 2.  a r²  Série 1 3,38132±1,56147x10-1 -6,46143x10-1±7,69229x10-2 0,86513 Série 2 3,39690±1,51409x10-1 -6,53410x10-1±7,48768x10-2 0,87378 Série 3 3,42155±1,58659x10-1 -6,67882x10-1±7,83656x10-2 0,86848 Série 4 3,38445±1,57815x10-1 -6,42574x10-1±7,92596x10-2 0,86795 Série 5 3,36999±1,35509x10-1 -6,36529x10-1±6,78789x10-2 0,89789 Completo 3,38914±1,42880x10-1 -6,55669x10-1±6,94184x10-2 0,88144 Fonte: Autora, 2015. Desse modo, os valores obtidos para os parâmetros  para a série completa correspondem a viscosidade de cisalhamento (em cP) obtidos em cada um dos modelos estudados, de acordo coma tabela 8. Tabela 8 – Coeficientes de viscosidade obtidos pelos 3 modelos estudados, calculados pelas suas TCAFs. Modelo Viscosidade (cP) Desvio (%) Modelo 0 4,95154 ± 2,16594x10-1 -72,5% Modelo 1 2,77208 ± 8,52263x10-2 -84,6% Modelo 2 3,38914 ± 1,42880x10-1 -81,2% Experimental 17,98934 Fonte: Autora, 2015. É possível constatar que o modelo 0 mostrou-se surpreendente melhor frente aos outros quando comparado com o resultado da viscosidade obtido experimentalmente, de acordo com a literatura. Este resultou em um menor desvio (–72,5 %) no valor experimental do que os outros modelos estudados. A figura 13 mostra uma comparação das curvas fitadas para os três modelos, demonstrando a eficiência relativa do modelo 0 frente aos demais testados. 49 Figura 13 – TCAFs globais dos 3 modelos estudados. Fonte: Autora, 2015. 50 3. CONCLUSÕES Por meio do estudo computacional por Dinâmica Molecular realizado, comparando os três modelos distintos de distribuição de carga para o líquido iônico tetracloroaluminato de 1-etil-3-metilimidazólio, conclui-se que, no que diz respeito ao comportamento difusivo anômalo, os resultados alcançados foram favoráveis ao modelo com transferência estática de cargas sem redistribuição (modelo 1). Isso se deve à modelagem de cargas nele empregada, que contribuiu para a redução das interações coulômbicas no sistema, alcançando comportamento difusivo adequado a tempos longos (10 à 15ns) e coeficiente de difusão mais próximo ao observado experimentalmente. Desse modo, essa modelagem favoreceu a mobilidade dos íons, principalmente em função redução da carga no anel imidazólio neste modelo. Já no que diz respeito ao coeficiente de difusão, comparativamente com o valor experimental relatado na literatura para este composto (8,170x10-7 cm2/s), todos coeficientes alcançados no estudo são inferiores a este, indicando que a transferência estática de carga se mostrou insuficiente para corrigir melhorar numericamente os coeficientes de difusão. Os modelos com transferência estática de carga (modelos 1 e 2) porém, em comparação com o modelo original (modelo 0), mostram-se mais próximos ao valor experimental, 5,988x10-8 cm2/s e 1,105x10-7 cm2/s para o cátion respectivamente, sendo o modelo 2 (também com transferência de carga) anda mais próximo do valor almejado, um desvio de -86,5%. Já em relação ao coeficiente de viscosidade, quando empregada a relação de Einstein, as flutuações no tensor de pressão não convergiram da maneira adequada, tornando a análise inviável devido a presença de oscilações. Já quando se empregou o método envolvendo o cálculo das TCAFs, curiosamente o modelo original (sem transferência de cargas, modelo 0) foi então o que mais se aproximou do valor apresentado na literatura, com um desvio de apenas –72,5 %. Novos estudos podem ser realizados a partir deste ponto, tais como: 1 – Identificar e sanar as causas dos pontos divergentes na análise pelo método TCAF; 51 2 – Verificar os resultados (e talvez sanar o item 1) obtidos pelas TCAFs gravando as configurações de trajetórias em intervalos de tempos mais próximos (a cada 1 ou 2 passos?); 3 – Aumentar o tempo das simulações de MD para as análises das VACFs, de modo que a integração adequada das funções seja possível (de 100ps para 25ns?); 4 – Verificar o comportamento difusivo a tempos longos identificado nas análises de MSD, através de estender ainda mais as simulações (de 25ns para 50ns?); 5 – Estudo comparativo com um modelo polarizável (Drude?), para verificar a influência da polarização e talvez alcançar resultados mais próximos aos constatados experimentalmente; 6 – Estudo comparativo das técnicas de análise de MD fora do equilíbrio para propriedades de transporte (em especial viscosidade) dos ILs, frente às técnicas empregadas aqui; 7 – Emprego de um modelo dinâmico de transferência de carga, para comparação com os modelos estáticos aqui empregados. 52 4. REFERÊNCIAS ABRAHAM, M. J.; SPOEL, D. VAN DER; LINDAHL, E.; HESS, B.; TEAM, G. D. GROMACS User Manual version 5.0.2. Www.Gromacs.Org, 2014. ADRIANE, C.; SEGATTO, R. Estudo computacinal da difusão de água em nanotubos de carbono induzida por gradiente de pressão. , p. 32–35, 2013. ANISIMOV, V. M.; VOROBYOV, I. V.; ROUX, B.; MACKERELL, A. D. Polarizable empirical force field for the primary and secondary alcohol series based on the classical drude model. Journal of Chemical Theory and Computation, v. 3, p. 1927–1946, 2007. ASSIS, P. C. DE. Universidade Federal do Rio Grande do Norte - PPGF. Processos Difusivos Generalizados. , 2006. BARREIRO, E. J.; RODRIGUES, C. R.; ALBUQUERQUE, M. G.; SANT’ANNA, C. M. R. DE; ALENCASTRO, R. B. DE. Modelagem Molecular: Uma Ferramenta para o Planejamento Racional de Fármacos em Química Medicinal. Química Nova, v. 20, n. 3, p. 300–310, 1997. BOES, E. S.; U.; Universidade Federal do Rio Grande do Sul. Estudo Computacional de Líquidos Iônicos do Tipo Imidazólio com Substituintes Insaturados. , 2012. BORODIN, O. Polarizable force field development and molecular dynamics simulations of ionic liquids. The journal of physical chemistry. B, v. 113, p. 11463– 11478, 2009. BRASIL, O. Química verde. , 2013. CONSORTI, C. S.; SOUZA, R. F. DE; DUPONT, J.; SUAREZ, P. A Z. Líquidos iônicos contendo o cátion dialquilimidazólio: Estrutura, propriedades físico-químicas e comportamento em solução. Quimica Nova, v. 24, n. 6, p. 830–837, 2001. BOURSCHEIDT, L. Universidade Federal do Rio Grande do Sul. Estudo Computacional de Líquidos Iônicos do Tipo Dialquilimidazólio. , 2005. Da SILVA,T. B., Universidade Federal de Pelotas. Líquidos iônicos - alguns aspectos sobre as propriedades , preparação e aplicações, 2004. De ANDRADE J.,FÍSICO-QUÍMICA. Universidade Federal do Rio Grande do Sul Líquidos Iônicos : Desenvolvimento de Campo de Força e Estudo das Propriedades Físicas e Estruturais. , 2008. FRANZOI, A. C.; BRONDANI, D.; ZAPP, E.; et al. Incorporação de líquidos iônicos e nanopartículas metálicas na construção de sensores eletroquímicos. Quimica Nova, v. 34, n. 6, p. 1042–1050, 2011. 53 GU, Y.; YAN, T. Thole model for ionic liquid polarizability. Journal of Physical Chemistry A, v. 117, p. 219–227, 2013. GUIDO, P. R. V. C.; MOLECULAR, M. FFI0776 Ͳ Modelagem g e Engenharia de Proteínas. . HESS, B. Determining the shear viscosity of model liquids from molecular dynamics simulations. Journal of Chemical Physics, v. 116, n. 1, p. 209–217, 2002. LEVASHOV, V. A. Atomic Level Green-Kubo Stress Correlation Function for a Model Crystal: An Insight into Molecular Dynamics Results on a Model Liquid. , 2014. Disponível em: . MARTÍNEZ L., R. Andrade, E. G. Birgin, J. M. Martínez. Packmol: A package for building initial configurations for molecular dynamics simulations. Journal of Computational Chemistry, 30(13):2157-2164, 2009. MEDINA, J. S.; PROSMITI, R.; VILLARREAL, P.; et al. Molecular dynamics simulations of rigid and flexible water models : Temperature dependence of viscosity. Chemical Physics, v. 388, n. 1-3, p. 9–18, 2011. Elsevier B.V. Disponível em: . NAMBA, A. M.; SILVA, V. B. DA; SILVA, C. H. T. P. DA. Dinâmica molecular: Teoria e aplicações em planejamento de fármacos. Ecletica Quimica, v. 33, n. 4, p. 13–24, 2008. NEYERTZ, S.; BROWN D. Local structure and mobility of ions in polymer electrolytes: A molecular dynamics simulation study of the amorphous PEOxNaI system. J. Chem. Phys. (Melville). 1996; 104: 3797 - 3809. PALMER, B. J. Transverse-current calculations. Physical Review E, v. 49, n. 1, p. 359, 1994. PAULA, F. DE. Introdução : características de uma interação iônica, 2006 . PEDRON, I. T.; MENDES, R. D. S. Difusão anômala e equações generalizadas de difusão. Revista Brasileira de Ensino de Física, v. 27, n. 2, p. 251–258, 2005. PRONK, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., Shirts, M. R., Smith, ´ J. C., Kasson, P. M., van der Spoel, D., Hess, B., Lindahl, E. GROMACS 4.5: a highthroughput and highly parallel open source molecular simulation toolkit. Bioinformatics 29(7):845–854, 2013. SCHRÖDER, C.; SONNLEITNER, T.; BUCHNER, R.; STEINHAUSER, O. The influence of polarizability on the dielectric spectrum of the ionic liquid 1-ethyl-3- methylimidazolium triflate. Physical chemistry chemical physics : PCCP, v. 13, p. 12240–12248, 2011. TAVARES, L.; COSTA, D. A. Simulação computacional de Eletrólitos Poliméricos em Poli (oxietileno) e Líquidos iônicos. Universidade de São Paulo, Brasil, 2007. 54 TAYLOR, T.; SCHMOLLNGRUBER, M.; SCHRÖDER, C.; STEINHAUSER, O. The effect of Thole functions on the simulation of ionic liquids with point induced dipoles at various densities. Journal of Chemical Physics, v. 138, n. 20, p. 0–10, 2013. WANG, F.; JORDAN, K. D. A Drude-model approach to dispersion interactions in dipole-bound anions. Journal of Chemical Physics, v. 114, n. 24, p. 10717–10724, 2001. WEISS, T. F. Estas notas de aula estão fortemente baseadas no livro de T. F. Weiss (2 vols.) indicado na bibliografia. , p. 1–20, 1869. WEISSHEIMER, M. I. K. Efeito da transferencia de carga em propriedades de líquidos ionicos contendo o cátion 1-Butil-3-Metilimidazólio. Journal of Chemical Information and Modeling, v. 53, p. 1689–1699, 2013. ZHANG, X.; HUO, F.; LIU, X.; et al. Influence of Microstructure and Interaction on Viscosity of Ionic Liquids. Industrial & Engineering Chemistry Research, v. 54, n. 13, p. 3505–3514, 2015. Disponível em: . 55 ANEXO 1 Lista de comandos utilizados: (GROMACS) Editconf, comando utilizado para converter arquivos de extensão .pdb em arquivos .gro (editconf –f .pdb –o .gro), bem como para ajustar a densidade da caixa de simulação (editconf –f .gro –o .gro –density -pbc). G_energy, comando utilizado para calcular a viscosidade pela relação de Einstein (g_energy –f .edr –s .tpr –o .xvg –vis .xvg). G_velacc comando utilizado para calcular as VACFs (g_velacc; - f .trr –s .tpr –mol – normalize –o .xvg). G_msd, comando utilizado para calcular o MSD (g_msd; -f .xtc –s .tpr –o .xvg). G_tcaf, comando utilizado para calcular as TCAFs (g_tcaf; - f .trr –s .tpr –ao .xvg –o .xvg –of .xvg –ov .xvg –b –e). 56 ANEXO 2 Para o cálculo da viscosidade através do comando g_tcaf foi feito um ajuste para rodar 25 000 ps em 5 séries (sendo cada uma delas com 5 000 ps), dando origem as seguintes tabelas para cada modelo correspondente. Os pontos removidos para o ajuste quadrático são mostrados com uma linha horizontal. Tabela – Modelo 0 K Série 1 Série 2 Série 3 Série 4 Completo 0,898 4,65393 4,546 4,49903 4,4961 4,55885 0,898 4,48789 4,55151 4,45559 4,62336 4,5347 0,898 4,411 4,58081 4,6018 4,79426 4,56896 1,271 2,85721 3,01837 4,89455 2,91219 2,92662 1,271 2,87905 3,06379 4,95947 3,06639 2,99586 1,271 3,04505 2,92439 3,01865 3,01806 3,00759 1,271 2,91458 2,99537 4,97943 2,97187 2,97116 1,271 2,91124 2,99841 2,96187 2,94664 2,96143 1,271 2,88216 2,96782 2,976 3,01337 2,96407 1,556 2,23281 2,31597 2,28169 2,28896 2,28473 1,556 2,25361 2,30502 2,33942 2,28782 2,30159 1,556 2,24571 2,36818 2,3094 2,30779 2,31237 1,556 2,26737 2,33717 2,29625 2,32346 2,31058 1,797 1,87242 1,88499 1,94058 1,89647 1,90268 1,797 1,83211 1,89817 1,92175 1,91411 1,89495 1,797 1,87782 1,87019 1,87886 1,934 1,89574 Fonte: Autora, 2015. 57 Tabela - Modelo 1 K Série 1 Série 2 Série 3 Série 4 Série 5 Completo 0,886 2,35132 2,36788 2,42047 2,41306 0,415138 2,38453 0,886 2,37254 2,42606 2,36494 2,4118 2,34818 2,38381 0,886 2,52376 2,36962 2,35498 2,31255 2,45327 2,401 1,253 1,71851 1,74156 1,69017 1,72398 1,725287 1,71973 1,253 1,67353 1,69726 1,73917 1,70588 1,67149 1,69612 1,253 1,71783 1,72968 1,70382 1,70385 1,67879 1,70566 1,253 1,74609 1,70071 1,71683 1,689 1,71493 1,71242 1,253 1,70645 1,72584 1,71412 1,66711 1,70953 1,70716 1,253 1,69003 1,68794 1,6768 1,74438 1,74577 1,7034 1,535 1,41117 1,43396 1,40323 1,36715 1,39505 1,40105 1,535 1,39298 1,37421 1,41738 1,4093 1,39332 1,39347 1,535 1,39298 1,38867 1,43694 1,39647 1,41224 1,40507 1,773 1,40699 1,41508 1,39129 1,40189 1,3968 1,40162 1,773 0,104072 16,5982 6,75631 1,17955 1,24919 16,9133 1,773 0,104072 1,66556 0,103859 12,3893 14,5537 15,1081 1,773 0,104072 21,3812 0,103859 0,563867 0,103784 0,10379 Fonte: Autora, 2015. 58 Tabela - Modelo 2 K Série 1 Série 2 Série 3 Série 4 Série 5 Completo 0,894 3,18512 3,16691 3,07304 3,15331 3,14857 3,13813 0,894 3,17546 3,12709 3,27328 3,21389 3,11014 3,17093 0,894 2,15565 2,12603 2,11509 2,22522 2,15597 2,15108 1,264 2,13654 2,19619 2,19463 2,18398 2,20681 2,17179 1,264 2,13059 2,13961 2,15791 2,16264 2,12962 2,13958 1,264 2,13901 2,19994 2,13758 2,11349 2,19571 2,15268 1,264 2,14482 2,15622 2,10511 2,1399 2,20781 2,14571 1,264 2,12772 2,10846 2,18726 2,09284 2,14804 2,12675 1,264 8,3314 7,03037 7,29037 6,98187 7,24324 1,71702 1,548 7,44163 7,17031 0,139428 0,139351 319,566 21,981 1,548 94,2914 1,74526 7,31415 0,139351 89,129 1,7307 1,548 1,73289 7,27712 1,71296 14,9063 0,13907 7,27962 1,548 1,48936 1,4836 1,45412 1,48529 1,44764 1,46803 1,787 1,48503 1,47691 1,47454 1,47433 1,477695 1,47483 1,787 1,46643 1,45913 1,44341 1,46774 1,44558 1,45356 1,787 6,72182 2,06877 4,62713 0,0464503 0,0463566 3,47882 Fonte: Autora, 2015.