UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE CENTRO DE CIÊNCIAS EXATAS E DA TERRA INSTITUTO DE QUÍMICA PROGRAMA DE PÓS-GRADUAÇÃO EM QUÍMICA Estudo de Clusters metálicos de Alumínio-Sódio, Alumínio-Potássio, Alumínio-Lítio e Sódio-Lítio pelas abordagens de Algoritmos Genéticos, Cálculos Quânticos e Análise Topológica Acassio Rocha Santos Dissertação de Mestrado Natal/RN, fevereiro de 2017 ACASSIO ROCHA SANTOS Estudo de Clusters metálicos de Alumínio-Sódio, Alumínio-Potássio, Alumínio-Lítio e Sódio-Lítio pelas abordagens de Algoritmos Genéticos, Cálculos Quânticos e Análise Topológica Dissertação de mestrado apresentada ao Programa de Pós-Graduação em Química, do Centro de Ciências Exatas e da Terra da Universidade Federal do Rio Grande do Norte, como requisito parcial obtenção do título de Mestre em Química. Linha de Pesquisa: Físico-Química e Química Teórica. Orientador: Prof. Dr. Caio Lima Firme Co-orientador: Prof. Dr. Jadson Cláudio Belchior Natal - RN Fevereiro de 2017 Universidade Federal do Rio Grande do Norte - UFRN Sistema de Bibliotecas - SISBI Catalogação de Publicação na Fonte. UFRN - Biblioteca Setorial Prof. Francisco Gurgel De Azevedo – Instituto Química - IQ Santos, Acassio Rocha. Estudo de clusters metálicos de alumínio-sódio, alumínio- potássio, alumínio-lítio e sódio-lítio pelas abordagens de algoritmos genéticos, cálculos quânticos e análise topológica / Acassio Rocha Santos. - 2017. 117 f.: il. Dissertação (mestrado) - Universidade Federal do Rio Grande do Norte, Centro de Ciências Exatas e da Terra, Programa de Pós- Graduação em Química, Natal, 2017. Orientador: Prof. Dr. Caio Lima Firme. Coorientador: Prof. Dr. Jadson Cláudio Belchior. 1. Físico-Química - Dissertação. 2. Algoritmos genéticos - Dissertação. 3. Potencial Gupta - Dissertação. 4. Clusters - Dissertação. 5. QTAIM - Dissertação. 6. Química - Dissertação. I. Firme, Caio Lima. II. Belchior, Jadson Cláudio. III. Título. RN/UF/BS-IQ CDU 544(043.3) DEDICATÓRIA A Deus, fonte de toda inspiração, que me motiva e me auxilia durante todas as etapas da minha vida. AGRADECIMENTOS A Deus, pelo seu amor, fidelidade, compaixão, força e consolo. À minha esposa, Ana Virgínia, pelo seu amor, força, motivação, auxílio, apoio e pelas revisões realizadas neste trabalho. Sem você, eu não teria conseguido chegar ao fim desta dissertação. À minha mãe, Zélia Rocha, pela confiança, força e amor. Ao meu pai, Dilmar Brito, por sempre acreditar em mim e nos meus sonhos. Ao meu orientador, professor Caio Lima Firme, pelas orientações, correções, pelo auxílio e pela disponibilidade durante todo o tempo do mestrado. Ao meu co-orientador, professor Jadson Cláudio Belchior, pela oportunidade de iniciação científica, mesmo sabendo da minha falta de tempo devido à realidade de trabalho e estudo; pelas orientações e contribuições desde o início da minha carreira acadêmica até os dias atuais. Aos professores Gerd Bruno da Rocha, Fabrício Gava Menezes e Robson Fernandes de Farias, por aceitarem o convite para a composição da banca de mestrado. Ao meu colega de laboratório, Fabrício, pelas contribuições para este trabalho. Ao meu amigo Daniel Marcos, pela amizade, pelo apoio e pelos conselhos, tanto em assuntos profissionais quanto pessoais. Fortificai as mãos desfalecidas, fortalecei os joelhos vacilantes. Dizei àqueles que têm o coração perturbado: “Tomai ânimo, não temais! Eis o vosso Deus!” (Isaías 35,3-4) Cada sonho que você deixa para trás é um pedaço do seu futuro que deixa de existir. (Steve Jobs) LISTA DE FIGURAS Figura 1: Fluxograma de operações de um GA padrão ........................................................... 20 Figura 2: Fluxograma do funcionamento do GA modificado ................................................. 23 Figura 3: Representação do VSCC da molécula de Amônia e dos pontos críticos (3,-3) e (3,-1) da molécula de água ....................................................................................................... 46 LISTA DE TABELAS Tabela 1 – Parâmetros de potencial Gupta para os sistemas Al-Li, Al-Na e Al-K ................... 26 Tabela 2 – Acrônimos, sinais dos autovalores e representação dos pontos críticos ................. 41 Tabela 3 – Valores de (b), ( 2b), |1|/3, G(rc)/(rc) de ligações covalentes C-C no etano, benzeno, etileno, das ligações iônicas em LiCl e NaF e do argônio diatômico ....................... 42 LISTA DE ABREVIAÇÕES GA: Algoritmo Genético OP : Operador FF: Fitness Function AN: Operador Aniquilador H: Operador História SEP: Superfície de energia potencial QTAIM: Teoria Quântica de Átomos em Moléculas DI: Índice de deslocalização BCP: Ponto crítico de ligação CCP: Ponto crítico de gaiola RCP: Ponto crítico do anel NACP: Atrator nuclear do ponto crítico NNACP: Atrator não nuclear do ponto crítico VSCC: Concentração de carga na camada de valência HF: Hartree-Fock DFT: Teoria do Funcional de Densidade MP2: Teoria de Perturbação de Møller-Plesset 2 CCSD: Método coupled cluster com substituições simples e duplas CCSD(T): Método coupled cluster com substituições simples, duplas e triplas D3BIA: Índice de aromaticidade baseado na densidade, deslocalização e degenerescência NICS: Deslocamento químico independente do núcleo RESUMO O estudo teórico de clusters metálicos tem despertado um interesse considerável, devido à possibilidade de criar novas ligas de materiais em nanoescala, as chamadas "nanoligas". Pesquisas sobre nanoligas desempenham papel significativo na Ciência de Materiais, pois, entre seus objetivos mais importantes, estão o de prever a estabilidade das estruturas, seus modos de crescimento, bem como o de auxiliar a interpretação de medidas espectroscópicas e outras medições experimentais. Nesse contexto, um grande número de métodos foi relatado nos últimos anos para a otimização do mínimo global de grupos atômicos e moleculares, sendo um dos mais utilizados atualmente o do Algoritmo Genético (doravante, GA), o qual baseia-se em princípios relacionados a processos evolutivos, em operadores inspirados na Teoria da Evolução e na Genética, isto é, na recombinação, mutação e seleção natural. Particularmente, o GA com a implementação do potencial Gupta tem se mostrado eficiente na busca de soluções “ótimas” em problemas de otimização de clusters metálicos. Esta dissertação é composta por capítulos de introdução, de metodologia, de abordagem teórica (Cap. 1, 2 e 3); e também por capítulos que contêm artigos sobre o tema proposto (Cap. 4, 5 e 6). No primeiro artigo (Cap. 4), analisaram-se clusters bimetálicos AlxNay (x+y≤55) por meio da aplicação do GA com a implementação do potencial Gupta. Com base também na aplicação do GA, no segundo capítulo (Cap. 5) foram estudados clusters de AlxLiy e AlxKy (x+y ≤ 55). Em ambos os trabalhos, para elevar a eficiência do GA, introduziu-se mais dois operadores: o Aniquilador e o História. Ao serem comparadas as estruturas obtidas por meio do GA com potencial Gupta para clusters de alumínio puro, lítio puro e alumínio-lítio com resultados recentes da literatura, verificou-se que para os sistemas Al2, Al3, Al6, Al8, Al9, Li5, Li6, Li7, Al1Li5, Al1Li7 e Al1Li8 as geometrias obtidas foram muito semelhantes àquelas resultantes de cálculos de funcional de densidade e ab initio[como CCSD(T)]. No terceiro artigo (Cap. 6), analisou-se um novo algoritmo genético quântico (Q-GA) para pequenos sistemas de clusters NaxLiy com (x+y ≤ 10). Constatou-se que o Q-GA apresenta maior eficiência na busca do mínimo global em relação ao GA com o potencial Gupta. Isso porque o primeiro utiliza método quântico, enquanto o segundo usa um método clássico. Por ser mais preciso, o Q-GA possui uma abrangência menor. Neste artigo, além de cálculos ab inito, também foram realizados cálculos topológicos a partir da Teoria Quântica de Átomos em Moléculas (QTAIM) para as estruturas Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1, obtidas pelo Q-GA. Nessas estruturas, chama a atenção o fato de não haver caminho de ligação envolvendo diretamente os metais, sendo unidos por pseudoátomos, com exceção do Na5Li1. Algumas interações atômicas não foram indicadas pelo caminho de ligação e sua análise foi feita pelo índice de deslocalização (DI). No sistema Na1Li5, os pares atômicos Na1-Li2 e Na1-Li6 têm as interações mais fortes (e equivalentes à do sistema NaLi) de todos os pares Na-Li de todos clusters NaxLiy(x+y=6); ao mesmo tempo, os outros pares Na-Li têm interações dez vezes mais fracas do que aquelas do sistema NaLi. As interações Na-Na dos clusters Na4Li2 e Na5Li1 são as mais fortes quando comparadas com sistemas puros. Por fim, verificou-se que a fórmula do grau de degenerescência do índice de aromaticidade D3BIA e a carga atômica indicaram que os átomos de lítio mais próximo ao átomo de sódio transferem carga para esse último. Palavras-chave: Algoritmo genético. Potencial Gupta. Clusters. QTAIM. ABSTRACT The theoretical study of metal clusters has drawn considerable interest due to the possibility of creating new alloys from materials in nanoscale, the so-called "nanoalloys". Research on nanoalloys has had an important role in materials science, since, among some of its most relevant objectives, we may find the prediction of stability in structures, their manners of growth and further assistance in the interpretation of spectroscopic and other experimental measures. In this context, several methods have been reported in the last few years towards the global minimum optimization of atomic and molecular groups, where the Genetic Algorithm (henceforth GA) is currently considered one of the most used methods, whilst based on principles related to evolutionary processes as well as operators inspired by the Theory of Evolution and Genetics, i. e., by recombination, mutation and natural selection. The GA method in particular, and altogether with the implementation of the Gupta potential, has become efficient in the search for “optimal” solutions for optimization problems in metallic clusters. The present dissertation is composed of chapters consisting of introduction, methodology and theoretical considerations (Chap. 1, 2 and 3), as well as of chapters containing articles on the proposed subject (Chap. 4, 5 and 6). In the first article (Chap. 4), we may find the analysis of AlxNay (x + y ≤ 55) bimetallic clusters through the Genetic Algorithm method with the implementation of the Gupta potential. Also based on the GA application, in the following chapter (Chap. 5) we may find a study regarding AlxLiy e AlxKy (x+y ≤ 55) clusters. In both works, in order to improve GA efficiency, two additional operators have been introduced: Annihilator and History. By being compared to structures obtained by means of GA with Gupta potential for pure aluminum, pure lithium and aluminum-lithium clusters in recent results from literature, it has been verified that, regarding systems Al2, Al3, Al6, Al8, Al9, Li5, Li6, Li7, Al1Li5, Al1Li7 e Al1Li8, the obtained geometries were very similar to those resulting from density functional and ab initio calculations [such as CCSD(T)]. In the third chapter (Chap. 6), we analyzed a new quantum genetic algorithm (Q- GA) for small cluster systems NaxLiy with (x+y ≤ 10). It has been observed that Q-GA presents an improved efficiency towards a global minimum regarding the GA with the Gupta potential. That has been the case since the former uses the quantum method, while the latter uses a classic method. More specifically, the Q-GA has a narrower scope. In this article, besides ab initio calculations, topological calculations were performed as well, grounded on the Quantum Theory of Atoms in Molecules (QTAIM) for the structures Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1 obtained by the Q-GA. In these structures, it is evident that there is no bonding path between the metals, since they are bonded by pseudo atoms, with the exception of the Na5Li1. Some of the atomic interactions have not been suggested by the bonding path, being their analysis performed according to the delocalization index (DI). In the Na5Li1 system, the atomic pairs Na1-Li2 and Na1-Li6 have the strongest interactions (equivalent to the NaLi system) of all Na-Li pairs in all of the NaxLiy (x+y=6) clusters; concurrently, other Na-Li pairs bear interactions ten times weaker than those from the NaLi system. The Na-Na interactions from the clusters Na4Li2 e Na5Li1 are stronger when compared to pure systems. Finally, it has been verified that the degree of degeneracy formula of the aromaticity index D3BIA and the atomic charge suggest that the lithium atoms that are closer to the sodium atom transfer charge to the latter. Key words: Genetic algorithm. Gupta potential. Clusters. QTAIM. SUMÁRIO CAPÍTULO 1 – INTRODUÇÃO E OBJETIVOS ................................................................... 12 1.1. Introdução ......................................................................................................................... 13 1.2. Objetivos ........................................................................................................................... 15 1.2.1. Objetivo Geral .......................................................................................................... 15 1.2.2. Objetivos específicos ................................................................................................ 15 CAPÍTULO 2 - METODOLOGIA........................................................................................... 16 2.1. Tratamento eficiente do sistema e custo computacional .................................................. 17 2.2. O Algoritmo Genético (GA): otimização de estruturas .................................................... 17 2.3. Algoritmo Genético padrão .............................................................................................. 19 2.4. Algoritmo Genético modificado ....................................................................................... 22 2.5. Potencial empírico funcional Gupta ................................................................................. 24 2.6. Energias associadas aos clusters....................................................................................... 26 2.7. O método do Algoritmo Genético quântico Q-GA .......................................................... 27 CAPÍTULO 3 – TEORIA E ABORDAGEM COMPUTACIONAL ....................................... 30 3.1. Teoria de Perturbação Møller-Plesset .............................................................................. 31 3.2. A partição de Møller-Plesset ............................................................................................ 32 3.3. As correções de energia MBPT(2) ou MP2...................................................................... 35 3.4. O método Coupled Cluster ............................................................................................... 36 3.5. Aproximação CCD ........................................................................................................... 37 3.6. As aproximações CCSD e CCSDT .................................................................................. 38 3.7. Teoria Quântica de átomos em moléculas (QTAIM) ....................................................... 39 3.8. Superfície interatômica e bacia atômica ........................................................................... 40 3.9. Pontos críticos .................................................................................................................. 40 3.10. Elipsidade ....................................................................................................................... 41 3.11. Classificação de interações químicas por QTAIM ......................................................... 42 3.12. Os índices de localização (LI) e deslocalização (DI) ..................................................... 44 3.13. Topologia do Laplaciano da densidade de carga ............................................................ 45 CAPÍTULO 4 – RESUMO ESTENDIDO PUBLICADO COMO ARTIGO GERAL NA REVISTA PROCESSOS QUÍMICOS (2015) - QUALIS C .................................................... 47 CAPÍTULO 5 – ARTIGO EM DESENVOLVIMENTO PARA SUBMISSÃO EM UMA REVISTA COM QUALIS B2, A SER SELECIONADA ........................................................ 65 CAPÍTULO 6 – ARTIGO EM DESENVOLVIMENTO PARA SUBMISSÃO EM UMA REVISTA COM QUALIS A1, A SER SELECIONADA ....................................................... 83 CAPÍTULO 7 - CONCLUSÕES GERAIS .......................................................................... 109 REFERÊNCIAS GERAIS ...................................................................................................... 112 12 Capítulo 1 – Introdução e Objetivos 13 1.1. Introdução Materiais nanoestruturados, ou “nanomateriais”, são classificados em: clusters de nanopartículas, nanopartículas, nanotubos, e materiais nanocompósitos. O controle da estrutura dos materiais a nível atômico, molecular ou supramolecular (nanoescala) fornece informações importantes acerca das propriedades macroscópicas, além de auxiliar a obtenção novas funcionalidades para solução de problemas específicos. Por exemplo, na construção civil, a aplicação de materiais nanoestruturados faz com que possam ser melhoradas características como a resistência mecânica, a durabilidade, a leveza e a flexibilidade; e introduzidas funcionalidades como o isolamento térmico, a autolimpeza. Outro exemplo de aplicações refere-se ao desenvolvimento de sensores para controle da segurança de estruturas.[1] Na Ciência dos Materiais, as propriedades de sistemas metálicos podem ser melhoradas por meio da mistura de elementos distintos para gerar compostos intermetálicos e ligas. Em diversas situações, ocorre a melhoria das características específicas sobre ligas devido aos efeitos sinérgicos. A rica diversidade de composições, estruturas e propriedades das ligas metálicas culminou em significativas aplicações em eletrônica, engenharia e catálise.[2] Clusters (nanopartículas) se definem como agregados de milhões de átomos constituídos por átomos idênticos, moléculas, ou ainda duas ou mais espécies diferentes. [2,3] Os clusters podem tanto ser investigados em diversos meios, tais como feixes moleculares, fase de vapor e suspensões coloidais; quanto isolados em matrizes inertes e superfícies. Em se tratando dos clusters metálicos, especificamente, grande parte do interesse no seu estudo decorre da possibilidade de criação de materiais em nanoescala, as chamadas nanoligas. Essas nanoligas constituem um novo tipo de material, cujas propriedades podem ser distintas dos átomos e moléculas ou da matéria macroscópica.[2] Entre essas propriedades, a forma geométrica e a estabilidade energética podem se alterar significativamente com o tamanho dos clusters. Por exemplo, sabe-se que os aglomerados de metais alcalinos, com tamanhos de até milhares de átomos - e clusters também menores de Cu e Ag -, em conformidade com o modelo jellium, são relativamente estáveis em certas nuclearidades - os chamados números mágicos. [4] Por outro lado, aglomerados de metais de transição e alguns metais do grupo principal - (como Al, Ca e Sr), geralmente, apresentam números mágicos que correspondem aos grupos constituídos por 14 poliedros concêntricos de átomos. Nesse sentido, há interesse contínuo em clusters metálicos, em razão do potencial de aplicações em áreas como catálise e nanoeletrônica (por exemplo, em dispositivos de tunelamento de um único elétron).[5] Os clusters constituídos por um pequeno número de átomos (2-100 átomos) são denominados de nano-clusters, os quais são mais adequados para uma abordagem teórica, pois possuem um menor número de partículas. Ao mesmo tempo, há dificuldades na realização de uma abordagem experimental do tema, como na determinação estrutural em relação ao número de partículas. A obtenção de isômeros de maior estabilidade em sistemas moleculares, clusters e estruturas cristalinas é de fundamental importância para o desenvolvimento de novos materiais, em que a estrutura mais estável de um cluster é aquela que possui a menor energia. Na determinação dos clusters, deve-se, a princípio, calcular a conformação de menor energia da superfície de energia potencial (doravante, SEP), ou seja, o mínimo global, o que geralmente pode ser considerado como um problema de otimização. Dada a grande quantidade de isômeros com estruturas distintas, a obtenção do mínimo global na SEP torna- se bastante dispendiosa. Foi constatado um crescimento quase exponencial do número de estruturas de mínimos à medida que o sistema aumenta de tamanho.[6,7] Com a finalidade de contribuir para a resolução das dificuldades mencionadas, diversos métodos de otimização foram desenvolvidos. Nos métodos Basin-Hopping[8,9] e nos Algoritmos Genéticos (GA,)[10] têm sido apresentados bons resultados, [8, 9, 11,12] ao passo que em métodos como Monte Carlo e Molecular Dynamics Simulated Annealing Approaches, as dificuldades em se encontrar os mínimos globais permanecem. [13,14] Neste trabalho, o enfoque se dará sobre as aplicações dos GA no estudo de clusters atômicos metálicos. Além de um método de otimização eficaz, como o GA, faz-se necessário também o uso de um potencial que descreva de forma satisfatória as interações interatômicas nos clusters metálicos. No modelo de sistemas metálicos, tem sido utilizado de modo satisfatório o potencial empírico Gupta,[1,15] que foi construído considerando-se o caráter essencial da banda de ligações metálicas, por meio da aproximação de segundo momento para o modelo tight-binding de ligação química.[16-18] A maior parte das pesquisas acerca de aglomerados metálicos focaliza a análise de clusters de metais alcalinos puros;[19-21] ao mesmo tempo, pesquisas sobre ligas metálicas têm sido pouco exploradas. Estudos realizados por Cheng et al. sobre ligas de alumínio-lítio tem despertado atenção para o tema, uma vez que o cluster de Al1Li5 é uma espécie 15 particularmente estável e desempenha o papel de um bloco de construção para os clusters LinAl.[22] Atualmente, são frequentes pesquisas mais detalhadas com enfoque em clusters de lítio puro,[23] alumínio puro[24] e clusters maiores de ligas de alumínio-lítio, com diversas nuclearidades.[25] Ligas de sódio-lítio possuem alta miscibilidade e, separando-se apenas na fase líquida a 573K, podem tanto solubilizar hidrogênio quanto serem promissoras na construção de anodos.[26–28] Abordagens teóricas e experimentais de clusters NaxLiy podem ser mencionadas. Mohajeri e Mahmoodinia,[29] por exemplo, avaliaram a estrutura dos isômeros para x+y=6. Benichou et. al.[30] relataram o potencial de ionização para diversas composições. Pérez et. al.[31] mediram a polarizabilidade de algumas estruturas selecionadas. Não foram encontradas na literatura análises de clusters de NaxLiy com nove e dez átomos. Na presente dissertação, será aplicado o GA acoplado ao potencial Gupta, no estudo de clusters de alumínio dopados com sódio, potássio e lítio. Além disso, um novo algoritmo genético com abordagem quântica será utilizado, na análise de pequenos clusters para sistema de sódio dopado com lítio. 1.2. Objetivos 1.2.1.Objetivo geral Validar os métodos do Algoritmo Genético com a implementação do potencial Gupta e o Algoritmo Genético Quântico (Q-GA). 1.2.2.Objetivos específicos a) Analisar o processo de crescimento dos sistemas AlxNay, AlxKy e AlxLiy (x + y ≤ 55), por meio do GA com potencial Gupta, de modo a se obter estruturas de menor energia, as quais servirão de condição inicial para os cálculos quânticos posteriores; b) realizar o tratamento quântico de alto nível, como CCSD(T) e MP2; c) realizar estudo de clusters através do Q-GA para o sistema NaxLiy (x + y ≤ 10); d) proceder à análise topológica e ao cálculo dos índices de aromaticidade NICS e D3BIA de alguns clusters do sistema NaxLiy (x + y = 6). 16 Capítulo 2 – Metodologia 17 2.1. Tratamento eficiente do sistema e custo computacional Nos últimos anos, surgiram diversos métodos considerados eficientes para determinar estruturas com mínimo global para clusters atômicos e moleculares.[32,33] Mesmo assim a extrapolação da SEP requer um método que atenda à complexidade de otimização de clusters maiores, em relação ao custo computacional e à elevação da quantidade de mínimos locais que aparecem na SEP com o aumento do cluster. [34,35] Alguns métodos apresentam falhas para otimização de clusters, uma vez que a presença de grande quantidade de mínimos locais mascara o mínimo global de um determinado sistema. O GA tem se revelado eficaz na resolução dessas falhas[32,33,35-37] uma vez que se trata de uma técnica de baixo custo computacional, pois não há a necessidade de cálculo de gradientes em toda a hipersuperfície.[38,39] O tratamento quântico também tem se mostrado eficaz no estudo teórico de pequenos clusters, todavia, o seu custo computacional é bastante elevado, pois a distribuição de possibilidades de estruturas iniciais é bastante ampla; além de o cálculo por meio de métodos ab initio para varrer todas as possibilidades ser muito caro, mesmo para clusters pequenos. Nesse sentido, um método menos robusto como o GA pode apresentar uma boa alternativa como condição inicial para um tratamento quântico posterior. 2.2. O Algoritmo Genético (GA): otimização de estruturas Os Algoritmos Genéticos (GA) são inspirados na Teoria de Darwin e estão relacionados à inteligência artificial no campo de computação evolucionária. Definem-se como um método para resolução de problemas por meio de um processo evolutivo, em que é obtida a solução mais bem adaptada (sobrevivente) dentro de um conjunto de todas as soluções possíveis, conforme a adequação à resposta procurada. Por essa razão, o GA é aplicado com bastante frequência na determinação de valores que correspondem aos mínimos da superfície de energia potencial associadas aos clusters.[33,35] O método GA pode ser comparado com as células presentes em um organismo vivo, no interior das quais os genes, blocos de DNA, compõem um mesmo conjunto de cromossomos e localizam-se em posições determinadas. Esses blocos são responsáveis por codificar certas proteínas, levando à caracterização do ser vivo, como cor de pele, tipo de 18 cabelo, cor dos olhos, etc. Os genes dos pais se combinam durante o processo de reprodução, de modo a formar um novo cromossomo, que, devido a alterações do DNA, pode sofrer mutações, causadas por diversos fatores, como, por exemplo, algum erro na cópia dos genes dos pais.[40] Na aplicação do GA, a fim de otimizar o sistema de uma população teste, inicialmente são usados operadores genéticos (OP) evolucionários, como o operador de Cruzamento (ou acasalamento) entre as espécies, de Mutação e de Seleção Natural. Em seguida, os indivíduos mais bem adaptados são ranqueados, por meio da Fitness Function.[41] Considera-se positivo o conjunto de soluções geradas pelo GA, porém, nem sempre é possível provar ou encontrar uma solução considerada ótima. Nesse sentido, com os objetivos de refinar a análise da busca por estruturas de clusters por meio da minimização da energia e de aumentar a probabilidade de encontrar estruturas de mínimo global na SEP, aborda-se nesta pesquisa o GA modificado, composto pelos operadores Aniquilador (AN) e História (H), além dos operadores já mencionados do GA padrão.[32,42] Em se tratando do funcionamento do GA, inicialmente a população é gerada de forma aleatória, segundo o tamanho e a composição do cluster, sendo as coordenadas espaciais criadas de modo randômico dentro de um determinado intervalo. Depois, o GA faz uso de uma relação de coordenadas cartesianas de cada átomo do cluster [38] para a geração de novos genes, os quais evoluem para a próxima geração através da aplicação dos operadores evolucionários de Cruzamento, Mutação e Seleção Natural. Após esse processo, cada indivíduo da população é deslocado para o seu mínimo local por meio do método BFGS e, assim, os operadores genéticos agem sobre a população atual. [43] Antes de iniciar a aplicação do OP Cruzamento, por meio da técnica de Roleta, os indivíduos pais são escolhidos segundo sua adequação ao sistema. Escolhidos os pais, o OP Mutação pode ou não atuar, uma vez que a frequência desse operador é pequena e ele é bastante lento. Sua aplicação ocorre somente em uma fração de cada população gerada. Desse modo, para os sistemas estudados nesta dissertação, a fração escolhida foi de 15%. Aplicados o OP Cruzamento e o OP Mutação, o OP Aniquilador remove indivíduos e mutantes com energia igual aos membros já existentes na geração; já o OP História “guarda na memória” informações sobre a população dessa geração para uso em uma etapa futura. Posteriormente, os demais indivíduos que sobraram são ranqueados de acordo com seu valor de Fitness Function. Com base nesse ranqueamento, o OP Seleção Natural faz o descarte dos últimos indivíduos da lista para que a população de clusters permaneça constante.[44] 19 O GA padrão constitui um loop que se repete por uma certa quantidade de vezes até atingir a convergência, quando não se conseguem estruturas com energia menor do que as existentes. Nessas circunstâncias, o método fornece o resultado do mínimo global encontrado, porém, isso pode não ser realmente verídico, dado o grande número de isômeros na SEP. O procedimento do GA modificado é semelhante ao do GA padrão, exceto pela inclusão dos OP Aniquilador e História. Nesse sentido, o GA modificado é um aprimoramento do GA padrão, realizando novos ciclos que ocorrem a partir do ponto em que o OP Aniquilador realiza a extinção de todos os indivíduos do sistema refinado pelo método padrão. O procedimento se reinicia, mas a população inicial não é mais totalmente aleatória, e sim promovida pelo OP História. Em seguida, um novo loop completo é repetido até que, após um certo número de vezes, nenhum isômero com energia menor que o último mais estável seja obtido. [32] Em face do exposto, percebe-se a inviabilidade de tratar um problema de clusters por meio de cálculos ab initio robustos com MP2, Coupled-Cluster ou por meio de DFT, dado que é expressivo o número de minimizações locais realizadas pelo GA. Desse modo, é mais interessante utilizar uma abordagem analítica com parâmetros experimentais para a energia potencial entre os átomos do sistema de muitos corpos. Expressões com essa abordagem são denominadas de “potenciais empíricos” de muitos corpos, bastante usados, pois são eficientes na modelagem de ligações entre átomos metálicos. Além disso, esses potenciais empíricos reproduzem com precisão propriedades estruturais e termodinâmicas na maior parte dos metais de transição. [18,32] Nas seções 2.3. e 2.4., serão detalhados o funcionamento do GA padrão e do GA modificado, respectivamente, com aprofundamento da breve descrição apresentada anteriormente. 2.3. Algoritmo Genético padrão Devido ao seu baixo custo computacional, o GA consiste em uma ferramenta largamente utilizada na otimização de clusters atômicos e moleculares. A seguir, é apresentado o fluxograma de funcionamento de um GA padrão. 20 Figura 1: Fluxograma de operações de um GA padrão [45] No método do GA padrão, são aplicados os operadores Cruzamento, Seleção Natural e Mutação, para a resolução de problemas de otimização. Tal aplicação se dá em uma população teste, de modo a obter mínimos globais para um determinado sistema. A população inicial é gerada aleatoriamente com base no número de clusters pré-definidos e na sua composição. Cada composição está associada a uma provável solução correspondente a um indivíduo da população do GA, sendo que a população inicial é obtida por meio das coordenadas cartesianas (x,y,z) que descrevem a posição espacial de cada cluster.[41] Gerada a população inicial e o tipo de indivíduos, ocorre a minimização local das energias dos clusters em relação às suas coordenadas. Em seguida, é aplicado o procedimento de seleção dos indivíduos mais bem adaptados ao final do processo evolucionário.[32] Terminadas essas etapas, a Fitness Function ranqueia os indivíduos conforme suas energias, sendo que o indivíduo com menor energia apresentará maior valor de Fitness Function, a qual verifica se eles são adequados ou não para o sistema. A Fitness Function é dada por:[46]  =  1 − ℎ 2  á  − 1 (), 21 em que Fi corresponde à medida da adequação do i-ésimo indivíduo da população que possui a energia Vi. O valor Vmin corresponde à energia mínima da população da vez, enquanto Vmax corresponde à energia máxima. Por meio da análise da equação, observa-se que o maior o valor de Fi está relacionado com o menor valor da energia do i-ésimo membro da população. Constitutivo do GA padrão, o OP Cruzamento possui a função de cruzar informações entre dois indivíduos dentro de uma população, a qual é avaliada pela Fitness Function para, então, serem selecionados dois membros para acasalamento. Os indivíduos mais bem adaptados(de menor energia) possuem maior probabilidade de serem selecionados, logo, têm maiores chances de repassarem suas informações para as próximas gerações. São geradas através do operador cruzamento uma, duas ou quatro descendências por par de pais. Os cruzamentos podem ser por cortes simples ou duplos. No primeiro, o corte é realizado em uma determinada direção, de modo randômico em cada um dos pais. No segundo, os clusters gerados são formados por duas partes: uma de um dos pais e uma do outro. [32,47] Por sua vez, o OP Seleção Natural possui a função de selecionar os indivíduos para o cruzamento. Ele apresenta diversos modos de atuação, sendo mais usada a técnica Roleta, na qual cada cluster da população é selecionado de acordo com o valor calculado pela Fitness Function. Os indivíduos mais bem adaptados têm maior probabilidade de serem escolhidos; porém, menos frequentemente, há probabilidade de indivíduos menos adaptados serem escolhidos. Apesar de apresentar baixa probabilidade e a ocorrência do fenômeno ser muito lenta, pode surgir após a aplicação do OP Cruzamento a ação do OP Mutação, que desloca, de maneira aleatória, uma fração de átomos em um cluster ou aplicando uma rotação em parte desse cluster. A depender das características do processo avaliativo, o resultado pode ser favorável ou não.[32,41] A principal contribuição do operador de mutação consiste na prevenção da convergência prematura para um mínimo local. Como o processo de mutação é um evento raro, adotou-se a taxa de 15% para a mutação em cada geração de clusters. Desse modo, aplicados os OP Cruzamento e Mutação, aplica-se o OP Seleção Natural, que seleciona indivíduos com menor energia para serem aproveitados na próxima geração e elimina os demais clusters existentes, de modo a manter constante a população. É importante frisar que a simples utilização dos OP Cruzamento e Mutação em indivíduos pertencentes a uma dada população poderia ocasionar a perda de clusters de menor energia. Para que isso não ocorra, o OP Seleção Natural foi desenvolvido para a selecionar 22 indivíduos com maior valor de Fitness Function em uma população intermediária, em que estão presentes os pais e os membros oriundos das OP Mutação e Cruzamento. Nesse sentido, o método GA torna-se mais efetivo na busca de mínimo global para um determinado sistema, pois o operador seleção natural age de forma elitista em vez de estocástica.[32,41] O critério de avaliação do sistema é feito por um operador que realiza a gestão do processo evolucionário. Por meio de condições previamente estabelecidas, caso certa população esteja caminhando no sentido da convergência, o processo sobre essa população é mantido até que aconteça a convergência. Entretanto, ressalte-se que, se após um determinado número de interações uma determinada população, não convergir, é gerada uma nova população no sentido de buscar outra solução para o sistema. Concluídas essas etapas, o processo evolucionário é finalizado, fornecendo a população final obtida para o sistema escolhido.[41] 2.4. Algoritmo Genético modificado Apesar da validade de todos os procedimentos do GA padrão, faz-se necessário neste trabalho utilizar novos operadores na etapa de avaliação, com o objetivo de tornar o método mais robusto[42] e aumentar a eficiência do GA. Propomos, a utilização de dois novos operadores, anteriormente mencionados, que não estão presentes no GA padrão: História (H) e Aniquilação (AN). Ambos possuem a função de elevar a probabilidade da localização do mínimo global na SEP. O processo evolucionário de GA modificado segue o esquema mostrado na Figura 2, a seguir: 23 Figura 2: Fluxograma do funcionamento do GA modificado.[32] A caixa tracejada à esquerda, em forma de retângulo, delimita a ação dos operadores do método do GA padrão. Já a elipse tracejada à direita delimita a ação dos OP Aniquilador e História, do método GA modificado. O OP Aniquilador atua após a aplicação dos OP Cruzamento e Mutação, com a função de remover descendências e mutantes presentes na subpopulação com energias iguais às de clusters pré-existentes na população. O OP História (H) age paralelamente ao AN, com a função de guardar na memória dados da população inicial. A partir desse processo, entra em ação o OP Seleção Natural, com a finalidade de remover os indivíduos menos adaptados, de maneira que a população de clusters permaneça constante. Em geral, após uma certa quantidade de iterações, ocorre a convergência, onde todos os indivíduos da população possuem a mesma energia. Nesse instante, o GA padrão interpreta o resultado como mínimo global obtido. Todavia, como a SEP apresenta um alto índice de mínimos locais, pode ser que a solução encontrada pelo GA padrão não se trate de um mínimo global e sim de um mínimo local. Na tentativa de minimizar esse problema, o OP Aniquilador 24 remove todos indivíduos da população encontrada no final do processo. O OP História age novamente, fornecendo uma nova população inicial, composta por indivíduos guardados de gerações passadas, exceto dos clusters que apresentaram convergência.[32,42] O processo é finalizado quando não é encontrado nenhum isômero mais bem adaptado após um certo número de ciclos. Na abordagem do GA modificado realizada neste trabalho, utiliza-se o máximo de 200 gerações e a quantidade de 30 clusters para a população inicial. Os parâmetros iniciais foram alterados de forma aleatória, com a realização de novos testes, porém o GA apresentou convergência equivalente, o que permite a inferência de que os sistemas analisados são independentes das condições iniciais. No presente trabalho, o sistema estudado através desse GA diz respeito ao crescimento de clusters de alumínio dopados com sódio, potássio e lítio. O sistema metálico é tratado através do potencial empírico clássico Gupta, obtido na literatura.[48] As estruturas mais estáveis e suas respectivas energias obtidas através do GA modificado servirão de condições iniciais para um tratamento quântico, que pode ser desenvolvido posteriormente. 2.5. Potencial empírico funcional Gupta Os potenciais empíricos de muitos corpos têm sido usados com frequência para solucionar diversos problemas, inclusive na Ciência dos Materiais. Eles reproduzem com boa exatidão diversas propriedades, tanto termodinâmicas quanto estruturais de boa parte dos metais de transição, entre as quais se destacam a energia de sublimação, a constante de rede, as constantes elásticas e a energia de formação de vacância.[49,50] Trata-se de um dos meios mais práticos na análise das propriedades de interface e de superfície de ligas metálicas, bem como na abordagem de simulações de defeitos pontuais ou estendidos.[40] As expressões empíricas de energia potencial para muitos corpos apresentam valores positivos quando comparadas com potenciais mais simples que incluem somente interações entre pares. Como a reprodução de características básicas de sistemas metálicos envolve a inclusão do caráter essencial da banda das ligações metálicas, o uso do potencial empírico de muitos corpos é mais apropriado.[18] Esse tipo de potencial analisa o efeito de cada entidade (j) do sistema sobre um dos componentes (i), além de avaliar a interação de cada par (ij). Ressalte-se também que é levado 25 em consideração o efeito causado na componente (i) pela influência que seus vizinhos exercem até certa distância de cada entidade (j). Vale a pena destacar ainda que são somadas todas das contribuições de i e j presentes no sistema. O modelo adotado no presente trabalho é o de potencial Gupta [50] da energia de coesão e das propriedades físico-químicas dos materiais, que está incluído no grupo de potenciais funcionais de pares e baseia-se no modelo do segundo momento de Tight (Tigth- Binding Second-Moment Model). Esse potencial é dado em função dos termos dos pares repulsivos Vr e do termo de acoplamento de muitos-corpos Vm, sendo que a energia de um determinado átomo i, calculada pela interação com todos os átomos vizinhos j, é escrito por:  = !"#( ) − $( )% (&)' , em que Vr descreve a repulsão de pares entre íons metálicos e é definido por: #( ) = !()*+ ,−+ -./.0 − 112 ' /3 (4) . O termo “de muitos-corpos”, Vm, descreve a coesão de muitos corpos da banda de valência, sendo dado por: $( ) = 5∑ 7)*+ 8−29 #:#; − 1<'/3 = >? (@). Os parâmetros A, ρ, r0, A e q foram obtidos através de dados experimentais de energia de coesão, parâmetros de rede e constantes elásticas independentes, da estrutura cristalina a 0 K. Esses dados foram calculados para os metais puros e também para as ligas formadas. O termo rij refere-se à distância entre os átomos i e j; já os parâmetros q e p são adimensionais e fazem uma descrição do alcance dos termos de banda e repulsão, respectivamente. Os parâmetros A e ζ apresentam unidades de energia e determinam a força de Vr e Vm. O termo r0 descreve a distância média dos vizinhos que possuem maior proximidade do átomo no bulk do metal. [51] 26 Os parâmetros para interações homoatômicas Al-Al, Li-Li, Na-Na e K-K presentes neste trabalho foram obtidos da literatura. [52-54] Para as interações heteroatômicas Al-Li, Al- Na e Al-K, efetuamos o cálculo de médias aritméticas simples.[51] Desse modo, o potencial Gupta pode ser utilizado tanto para o modelagem de diversos sistemas sólidos e clusters, quanto para ligas metálicas. Na Tabela 1, são apresentados os valores dos parâmetros utilizados nesta dissertação. Tabela 1: Parâmetros de potencial Gupta para os sistemas Al-Li, Al-Na e Al-K 2.6. Energias associadas aos clusters Uma grandeza de grande aplicabilidade no estudo de clusters é a energia média de ligação, dada por: BC = DEFG(H)H (I). Se os clusters possuírem tamanho considerável, o valor de Eb deve aproximar-se assintoticamente do valor da energia de coesão do bulk relativo aos átomos que compõem o cluster.[32] Outra análise também frequente no estudo de clusters é a segunda diferença de energia média das ligações. [55,56] Seu cálculo é realizado através da seguinte expressão: ∆BC(K) = 2BC(K) − BC(K − 1) − BC(K L 1) (M). 27 A abordagem da segunda diferença de energia permite avaliar o comportamento energético de um dado cluster com N átomos em relação a seus “vizinhos”, ou seja, um cluster com N + 1 e outro com N – 1 átomos. Se ∆BC(K) > 0, encontramos picos de estabilidade, já se ∆BC(K) < 0, aparecem picos de instabilidade. [32] Embora a segunda diferença de energia seja bastante utilizada, a sua eficiência nem sempre descreve bons resultados, pois a comparação da estabilidade está associada somente com os clusters vizinhos. Podem haver situações em que os clusters vizinhos são energeticamente instáveis, ocasionando um pico estabilidade dos aglomerados de N átomos. Nesses casos, a estrutura apresentou pico de estabilidade porque os vizinhos eram muito instáveis energeticamente. Outra forma de avaliar clusters de ligas metálicas é a energia de excesso, a qual fornece informações sobre o quão favorável ou não é a formação da liga quando comparada ao cluster puro correspondente. [2,57] O cálculo da energia de excesso para um cluster de N átomos de uma liga bimetálica AxBN-X é dado por: BRS((TUHT) =  ((TUHT) − V  ((H)K − (K − V) (UH)K (W), em que valores negativos para BRS((TUHT) favorecem a formação da nanoliga correspondente. No nosso trabalho, são realizados cálculos de energia de excesso para os sistemas AlxLiN-x, AlxNaN-x e AlxKN-x com 3 ≤ N ≤ 55. 2.7. O método do Algoritmo Genético quântico Q-GA O Q-GA desenvolvido pelo grupo LDAM1 realiza as otimizações dos clusters através de cálculos quânticos com o uso do pacote GAMESS-US. O tratamento quântico é efetuado usando o método MP2, com o uso de um potencial nuclear efetivo e base LANL2DZ. O referido método é inicializado com uma geração aleatória dos clusters da população inicial, sendo o número de indivíduos definido pelo usuário. Posteriormente, os clusters são otimizados localmente pelo método BFGS, que está presente no pacote GAMESS-US. Como os indivíduos são gerados aleatoriamente, a probabilidade de estarem próximos do mínimo 1 LDAM – Laboratório de Dinâmica Molecular –, do Departamento de Química da Universidade Federal de Minas Gerais. O autor deste trabalho e o seu co-orientador fazem parte desse grupo. 28 global é muito remota. Nesse sentido, com o objetivo de reduzir o custo computacional, os critérios de convergência do gradiente de energia calculado pelo GAMESS-US iniciam-se com altos valores. Nessa etapa inicial, as estruturas com a maior componente do gradiente, com valores menores ou iguais a 0,5, são aceitas como mínimos. Posteriormente, os clusters são classificados em ordem crescente de energia, formando-se um ranking, em que as estruturas que não alcançarem convergência situam-se nas últimas posições. Em uma nova etapa, o OP Predador age, eliminando 25% dos indivíduos com maiores valores da energia potencial da população inicial. Para manter a população constante, três operadores possuem a função de gerar novos indivíduos de formas diferentes. O operador Imigração gera novos indivíduos de forma totalmente aleatória, do mesmo modo como descrito para a geração da população inicial, sendo as estruturas geradas otimizadas localmente pelo método BFGS. Esse operador simula a imigração de indivíduos para o meio submetido ao processo evolucionário, conferindo uma diversidade imparcial à população. Dentre o total de clusters que devem ser incorporados à população para mantê-la constante após a ação do Predador, 5% é incorporado pelo OP Imigração. Se os 5% não corresponderem a um número inteiro de indivíduos, ocorre o arredondamento para o inteiro mais próximo. Caso 5% seja menor do que a unidade, pelo menos um indivíduo é adicionado à população. O operador Mutação é então aplicado sobre os clusters que não foram eliminados pelo Predador; ele realiza deslocamentos aleatórios nas coordenadas dos átomos dos clusters e gera mutantes, mas preserva os originais. Os indivíduos mutantes são otimizados localmente pelo GAMESS-US. Em seguida, o OP Cruzamento atua também de modo aleatório, selecionando pares de indivíduos restantes da população inicial que servirão como genitores de novos clusters. Nesse contexto, a ocorrência do crossover tem como base um método proposto recentemente na literatura, que consiste na utilização de uma esfera para o corte e ligação das estruturas, ao invés de um plano como habitual. [58] Cada genitor é dividido em duas partes, sendo uma na parte externa e outra no interior da esfera. Se a parte interna dos genitores possuírem o mesmo número de átomos, ocorre a permutação entre essas partes.[56] Se a aplicação da esfera não for possível, o crossover se dá pela média aritmética das coordenadas dos pais, sendo que, em seguida, as estruturas geradas são novamente otimizadas pelo método BFGS do GAMESS-US. Nesse momento, ocorre a reintegração da população inicial para 29 cada geração. Do número total de indivíduos que devem ser incorporados à população corrente, após a ação do Predador, 20% vem do OP Mutação e o restante da população é incorporado pelo OP Cruzamento. Depois da ação de todos os operadores descritos, é realizado um novo ranqueamento dos clusters, em que o indivíduo mais bem adaptado é avaliado pelo operador Highlander, cujo objetivo é analisar se o tal indivíduo foi o de menor energia potencial (primeiro do ranking) após um determinado número de gerações consecutivas. Em caso afirmativo, dizemos que o Highlander perdurou. Assim, a convergência do Q-GA é enfim alcançada. Em caso negativo, o Q-GA verifica se foi atingido o número máximo de gerações permitidas definidas pelo usuário. Se tal número não for atingido, o procedimento se reinicia a partir da ação do Predador. Caso contrário, o Q-GA é interrompido e não se assume nenhuma convergência. Por último, conforme o Q-GA avança para números maiores de gerações, os critérios de convergência para a otimização feita pelo GAMESS-US ficam cada vez mais rigorosos, até atingir o valor de 0,0005 como máximo aceito para a maior componente do gradiente de energia da estrutura considerada como mínimo global. 30 Capítulo 3 – Teoria e Abordagem Computacional 31 3.1. Teoria de Perturbação de Møller-Plesset Quando o tratamento de um sistema é realizado para o caso especial de apenas um determinante, a energia de Hartree-Fock(EHF) apresenta seu melhor resultado. No cálculo da energia de um sistema pelo método Hartree-Fock, há um erro associado à energia de correlação, característico de método variacional, que se baseia em um único determinante. A interpretação física é a de que nesse método as interações eletrônicas são tratadas pela média autoconsistente, não havendo mais a contribuição do movimento entre as partículas correlacionadas. A energia de correlação é obtida pela diferença entre a energia do método Hartree-Fock restrito e a energia exata E0 não relativística do sistema:[59] Ecorr = E0 - EHF (8). Na literatura, há uma grande variedade de métodos teóricos desenvolvidos para a obtenção da energia de correlação eletrônica, entre os quais destacamos o método de Interação de Configurações (CI)[60,61] e a Teoria de Perturbação de Muitos Corpos (MBPT – Many Body Perturbation Theory).[62] O prêmio Nobel de Química de 1998, conferido a J. Pople,[63] foi principalmente relacionado à aplicação dessas teorias. A Teoria de Perturbação de Muitos Corpos, desenvolvida por Schrӧdinger em 1926,[64] é considerada ainda nos dias atuais como uma das mais importantes ferramentas para o estudo de física de muitos corpos. Das formulações perturbativas mais simples, destaca-se a Teoria de Perturbação de Rayleigh-Schrӧdinger (RSTP), cujo princípio-chave refere-se à divisão do Hamiltoniano em duas partes: o Hamiltoniano não perturbado, que consiste na parte fundamental com as autofunções conhecidas; a perturbação de pequena escala, sendo o resultado exato não diferenciado significativamente da solução não perturbada. A energia exata é obtida pela soma das infinitas contribuições, denominadas ordens de perturbação.[65] A forma de perturbação mais utilizada para tratamentos de sistemas envolvendo Química Quântica refere-se ao uso da RSTP com a partição do Hamiltoniano realizada por Møller-Plesset.[66] A essa abordagem se dá o nome de Teoria de Perturbação de Møller- Plesset (MP)[67] ou (MBPT).[68] Com a necessidade do avanço para ordens superiores às dos métodos (MP), surgiu uma técnica alternativa, denominada “coupled cluster” (CC), em que determinadas contribuições da série de perturbação podem ser efetivamente somadas até a ordem infinita. [65] 32 3.2. A partição de Møller-Plesset Para o desenvolvimento de uma teoria de perturbação para sistemas moleculares, faz- se necessária a escolha da partição do Hamiltoniano eletrônico, dado por: X = − ! ħ2Z H [ ∇ − !! ]^)4`a0.^ b ^[ H [ L !! 14`a0./ H /c H [ (d), onde m corresponde à massa do elétron, ZA ao número atômico do núcleo A, N ao número de elétrons do sistema e M indica a quantidade de núcleos do sistema. Ao se utilizar o sistema de unidades atômicas em que a distância é dada em raios de Bohr e a energia em hartree, o hamiltoniano pode ser reescrito do seguinte modo: X = − !12 H [ ∇ − !! ]^.^ b ^[ H [ L !! 1./ H /c H [ (e) . Uma alternativa interessante para realização dessa escolha é definir o Hamiltoniano não perturbado como: X(0) = !ℱ(g)H[ (), onde ℱ é o operador de Fock, definido por: ℱ(g) = ℎ(g) L ![ℐC (g) − j(g)] (&), onde: ℎ(g) = − 12 ∇ − ! 1.^ b ^[ (4) 33 O termo ℐ é o operador de coulomb e j é o operador de troca. Observa-se que a soma em g ocorre sobre os elétrons do sistema, processo de seleção denominado de partição de Møller-Plesset.[66] Nesse sentido, a discussão a seguir tratará de spin-orbitais, com a obtenção de dados válidos tanto para RHF[69] quanto para UHF.[70] Aplicando o H0 na função de onda HF, Φ0, tem-se: X0Φ0 = !ϵnΦ0 (@)n A autofunção de H(0) é dada por Φ, com autovalor ∑ ϵn n . Os spin-orbitais encontrados no determinante HF são designados pela soma em c. Utilizado o conceito de determinantes substituídos, o cálculo SCF é feito com um grupo k de funções de base, obtendo-se 2k spin- orbitais moleculares. Desses, são inclusos no determinante HF apenas os N spin-orbitais de menor energia, que são os ocupados. Já os spin-orbitais 2K – N que sobraram são denominados spin-orbitais virtuais. A troca de um ou mais spin-orbitais preenchidos (a,b,c,..,) e spin-orbitais virtuais (r,s,t,...) resulta em um determinante substituído. Quando o spin-orbital a preenchido for substituído pelo spin-orbital virtual r, o determinante simplesmente substituído é designado por Φo# . De maneira análoga, os determinantes duplamente substituídos, indicados por ΦoC# , são aqueles que os spin-orbitais ocupados a e b foram substituídos pelos spin-orbitais r e s virtuais. Ressalte-se ainda que os determinantes substituídos são autofunções de H(0). Seguem os modelos para os determinantes com substituição simples e duplamente substituídos, respectivamente: X0Φo# = ,!a  − ao L a#2 Φo# (I) , X0ΦoC# = ,!a  − ao − aC L a# L a 2 ΦoC# (M). O hamiltoniano não perturbado, obtido pela equação (11), resulta na perturbação: 34  = X − X0 = !ℎ(g) L ! 1./q/ −!(g) (W), ou por:  = ! 1./q/ − !rst (g) (u), onde: rst(g) = ![v(g) − j(g)] (d) . O primeiro termo de V fornece a exata interação elétron-elétron; já o segundo termo fornece o dobro dessa interação na forma de média. É preciso destacar que V possui a função de reparar a dupla contagem da interação elétron-elétron em H0, além de apresentar com detalhes as interações individuais entre os elétrons. Uma maneira de se calcular as correções perturbativas na energia é, inicialmente, obter elementos de matriz dos tipos wΦxyxΦz{ e wΦxyxΦz{, em que Φ e Φz descrevem as autofunções de H0 (determinantes HF e substituídos): y = !rst (g) (&e), e y = ! 1./q/ (&), De modo que: wΦxxΦz{ = wΦxy − yxΦz{ (&&). 35 3.3. As correções de energia em MBPT (2) ou MP2 A obtenção das primeiras correções na energia de um sistema de muitos elétrons pode ser analisada com base na equação (23), em que B0(0) é: B0(0) = !ϵn (&4).n A correção de primeira ordem é realizada pela equação: |Φ0||Φ0~ = |Φ0|y − y|Φ0~ = −12 !| ‖ ~o, (&@), de forma que a energia corrigida até primeira ordem expressa-se por: B0(0) L B0() = !ϵn − 12 !| ‖ ~o, (&I).n Como esse resultado é justamente a energia de Hartree-Fock, isso significa que as correções de energia para o MBPT ou MP só ocorrem a partir da segunda ordem. A correção de segunda ordem é assim representada: B0() = !! || ‖.‚~|ϵo L ϵC − ϵ# − ϵ #q oqC (&M), em que as restrições <  e . < ‚ possuem a função de evitar que uma mesma substituição seja inclusa em duplicidade. Desse modo, o método MP2 ou MBPT(2) é dado pela seguinte expressão: Bƒ„ = B…† L B0() (&W). Tal expressão inclui somente substituições duplas. Como a correção de primeira ordem apenas recupera a energia HF, tem-se que o MP2 ou MBPT(2) é o método de menor ordem na abordagem MBPT. Como o sinal de B0() é sempre negativo, a energia do MP2 é 36 menor do que a energia do método Hartree-Fock. Na maioria das abordagens, B0() superestima a energia de correlação, ocasionando uma energia menor do que a energia exata obtida por MPBT de ordem superior ou coupled cluster. Apesar da teoria de perturbação tratar-se de um método bastante difundido, há situações em que a série de perturbação é divergente,[71] sendo uma opção mais interessante o tratamento pelo método coupled cluster. 3.4. O método Coupled Cluster O método coupled cluster (CC) é uma alternativa eficaz na busca pela energia de correlação eletrônica. Esse método, bastante difundido nos últimos anos, foi formulado por Coster[72] e Kümel[73] e desenvolvido por Čížek.[74,75] Ele consiste em separar um sistema de muitos elétrons em diversos aglomerados (clusters) com uma quantidade menor de elétrons. Em seguida, são calculadas das interações entre os elétrons do mesmo cluster e, posteriormente, é calcula a interação entre clusters distintos. A função de onda do método CC é obtida do seguinte modo: [74,75] Φ = )‡Φ0 (&u), em que T é o operador de cluster, representado por: ˆ = ˆ L ˆ L⋯ ˆŠ (&d), sendo: ˆ = ! o#o,# .‹ (4e), ˆ = !! oC# .‹‚‹  (4)#q oqC . Nesse modelo, , , … estão associados aos orbitais ocupados no determinante HF; já r,s,... descrevem os orbitais virtuais. Os coeficientes t são valores reais, chamados de 37 amplitudes de cluster. As expressões T1 e T2 são operadores que descrevem as configurações simples e duplamente substituídas. Vale destacar que os valores de t são ajustados de modo que a função de onda Φ seja a solução da equação de Schrӧdinger. Desse modo, tem-se: X)‡Φ0 = B)‡Φ0 (4&). Multiplicando-se o lado esquerdo da igualdade por )‡, obtém-se o resultado: )‡X)‡Φ0 = BΦ0 (44). É possível demonstrar[76] que )‡X)‡é truncada em quarta ordem em T. No cálculo CC, não é possível efetuar a inclusão de todas as ordens de substituições, pois o que ocorre é o truncamento T em algum Tp , com pequenos valores de p. 3.5. Aproximação CCD A aproximação coupled cluster com substituições duplas (CCSD) é a mais simples dentro da abordagem CC. No operador de cluster é realizada a inclusão apenas do termo T2. Esse abordagem foi proposta por Čížek[74,75] e implementada por Bartlett e Purvis[77,78] e Pople et al.[79]. Nessa abordagem, a equação de Schrӧdinger é obtida por: )‡?X)‡?Φ0 = BΦ0 (4@). Ao multiplicarmos essa equação por |Φ0|, tem-se: B = Bst L !! oC# | ‖.‚~#q oqC (4I). A principal dificuldade do método CC está na resolução das equações para as amplitudes do cluster, bastante trabalhosas, pois incluem todos os coeficientes das equações envolvidas. Logo, essas equações têm que ser solucionadas de forma autoconsistente. Dessa 38 forma, o método de Newton-Raphson foi o mais utilizado para a abordagem desse problema.[80] Aplicando-se diversas aproximações nos coeficientes das amplitudes, o método CC recupera os resultados para energia dos métodos MBPT(2), MBPT(3) e MBPT(4). Nesse sentido, o método CCD, apesar de apresentar contribuições duplas e quádruplas para a energia de correlação, possui a limitação de não considerar substituições simples, triplas, etc. Uma maneira de contornar esse problema é incluir no operador do cluster também as substituições simples (CCSD) e as contribuições simples e triplas (CCSDT), com aumento do custo computacional. 3.6. As aproximações CCSD e CCSDT Embora as substituições duplas apresentem maior domínio nos cálculos de energia de correlação em relação a outras substituições, há situações em que as substituições simples são de grande importância, como aquelas que envolvem processos de uma partícula como potencial de ionização e momento dipolo. O método coupled cluster como substituições simples e duplas (CCSD) foi desenvolvido por Purvis e Bartlett.[81] Sua função de onda dá-se por: ΦnnŽ = )‡>‹‡?Φ0 (4M). Vale destacar que, no CCSD, além das substituições simples e duplas, há também a da inclusão de substituições triplas, quádruplas, etc., oriundas de termos desconexos como ˆˆ, ˆ, ˆ, etc. Uma aproximação CC ainda mais robusta, implementada por Noga e Ballet,[82] é a de que inclui substituições simples, duplas e triplas no operador de cluster, denominada CCSDT. Sua função de onda é descrita como: ΦnnŽ(‡) = )‡>‹‡?‹‡Φ0 (4W). Na prática, a inclusão das substituições triplas do CCSD(T) é realizada de forma parcial, porém, mesmo assim, esse é o método de referência para cálculos de camada fechada. 39 O CCSD(T) é o método mais atraente atualmente, apresentando resultados altamente competitivos; entretanto possui elevado custo computacional, o que o torna um método bastante dispendioso.[81] Por essa razão, diversos modelos aproximados mais acessíveis computacionalmente têm sido desenvolvidos nos últimos anos.[65] 3.7. Teoria Quântica de átomos em moléculas (QTAIM) Richard Bader [83] propôs na década de 1960 uma teoria baseada na partição da densidade eletrônica e no estudo de dados experimentais, denominada de Teoria Quântica de Átomos em Moléculas (QTAIM), que pode ser analisada através das derivadas de densidade eletrônica (∇ρ e ∇2ρ). Na QTAIM, a soma das propriedades atômicas de uma determinada molécula descreve as propriedades dessa molécula. Uma limitação dessa teoria é que a caracterização de ligações sigma e pi não pode ser realizada de forma direta, pois a QTAIM não descreve orbitais moleculares. Por meio de dados obtidos por difração de raio-X e da topologia obtida pela QTAIM, Hofmann et al.[84] obtiveram valores similares de densidade de carga nos pontos críticos das moléculas de N(C4H9)B6H7. Do mesmo modo, Matta e Bader[85] relataram diversos exemplos que apresentaram correlação com propriedades medidas experimentalmente. Sendo assim, a QTAIM demonstrou ser um método bastante eficaz na descrição detalhada de propriedades intrínsecas de sistemas moleculares. Um dado primordial na QTAIM é caminho do gradiente,[83,86,87] expresso por: ∇φ = ’S “”“* L ’• “”“– L ’— “”“˜ (4u), onde ’S, ’• e ’— são vetores unitários e ∇φ apontam sempre para a direção de aumento de φ. No que concerne à propriedade do gradiente, a mais interessante é aquela em que os seus caminhos passam pelos pontos ortogonais à superfície, com valores constantes através de curvas de isodensidade. A QTAIM utiliza a matriz de densidade que foi calculada anteriormente por algum método ab initio ou DFT como ponto de partida. A partir dessa informação, calcula novamente a matriz de densidade, porém dividindo o sistema molecular em subsistemas, 40 contendo átomos da molécula de origem.[88] É importante que o método ab initio apresente exatidão considerável, ou seja, que o erro de correlação eletrônica seja muito pequeno. É possível estimar na QTAIM praticamente todas as propriedades quanto-mecânicas médias dos átomos envolvidos, de forma isolada. A soma das propriedades de todos os átomos individualmente resulta na propriedade média molecular.[89] A forma de um determinado sistema molecular é determinada pelas forças atuantes nesse sistema. Como a força que exerce maior domínio é a atração núcleo-núcleo, temos que a densidade ρ(r; X) constitui uma propriedade topológica importante, uma vez que o seu valor representa uma máximo local entre os núcleos. Porém, quando um núcleo e um elétron se unem fortemente, o potencial de Coulomb torna-se altamente negativo, ocasionando a exibição de um cúspide em uma posição do núcleo. 3.8. Superfície interatômica e bacia atômica O subconjunto aberto do subespaço tridimensional constitui a denominada bacia atômica. Em uma molécula, o átomo é representado pela junção entre um atrator e sua bacia associada. A separação de átomos vizinhos ocorre pela superfície interatômica ou superfície de fluxo zero, a qual possui a propriedade específica que a distingue das demais superfícies arbitrárias. Em todos os pontos na superfície de fluxo zero, o vetor n é ortogonal a ∇ρ; logo ∇ρ.n =0. Assim, na superfície de fluxo zero, não ocorre o cruzamento das trajetórias de ∇ρ. 3.9. Pontos críticos A classificação dos pontos críticos ocorre de acordo com o conjunto de autovalores da matriz Hessiana ∇∇ρ. A expressão de autovalor da matriz Hessiana apresenta três soluções (λ1, λ2, λ3) que estão relacionadas ao autovetor i (i=1,2,3), por sua vez associado aos eixos principais da curvatura. Embora os autovalores da matriz Hessiana sejam reais, eles podem assumir valor nulo. Em um ponto crítico, o ranqueamento (ω) é dado pelo número de autovalores ou curvaturas não nulos de ρ, já a assinatura (σ) é dada pela soma dos sinais dos autovalores, conforme Tabela 2, a seguir: 41 Tabela 2: Acrônimos, sinais dos autovalores e representação dos pontos críticos Nome Acrônimo λ1 λ2 λ3 (ω,σ) Atrator nuclear NA - - - (3,-3) Ponto crítico da ligação BCP - - + (3,-1) Ponto crítico do anel RCP - + + (3,+1) Ponto crítico da gaiola CCP + + + (3,+3) O gráfico molecular consiste em um grupo de caminhos de ligação que ligam pares vizinhos de atratores nucleares.[90] Tal gráfico resulta das principais propriedades topológicas de um sistema de distribuição de carga. Os pontos críticos (3,-3) acontecem na posição dos núcleos, ao passo que os pontos críticos (3,-1) ligam determinados pares de núcleos na molécula. O ponto crítico (3,+1) corresponde ao ponto crítico de anel e é encontrado no seu interior somente se os caminhos de ligação formarem um anel. O ponto crítico (3,+3) é encontrado no interior de uma gaiola, desde que os caminhos de ligação reúnam-se a formar superfícies anulares no formato de gaiola. A relação entre o número e a classe de pontos críticos que podem estar presentes em um sistema é obtida pela relação de Poincaré-Hopf (n –b + r – c = 1), em que n representa o número de pontos críticos (3,-3), b o número de pontos críticos (3,-1), r o número de pontos críticos (3,+1) e c o número de pontos críticos (3,+3).[91] Vale destacar que a caracterização experimental dos pontos críticos de uma dada molécula pode ser observada comparando-se o refinamento do multipolo da análise de difração de raios-X com a análise topológica obtida por cálculos QTAIM.[92] 3.10. Elipsidade A elipsidade consiste em uma propriedade bastante útil na QTAIM, que é dada pela medida [(λ1/λ2)-1] na superfície interatômica de um ponto crítico de ligação. Por meio da elipsidade é possível dimensionar o grau da densidade eletrônica presente em um dado plano. 42 Em uma ligação C-C no etano ocorre a equivalência simétrica dos eixos associados λ1 e λ2, de forma que a elipsidade assuma valores próximos de zero. Assim, ela fornece informações quantitativas sobre o caráter pi da ligação C-C. Em seus estudos, Propelier[93] verificou uma elevação considerável da elipsidade nas ligações duplas do etileno e do 1,3-trans-butadieno, além de ter observado que a elipsidade do benzeno assume valores entre uma ligação simples e uma dupla. Foi observado também que, quando o tratamento ocorre em ligações triplas, a elipsidade falha em sua descrição, uma vez no acetileno ela assume valor nulo. Nesse sentido, para descrever o caráter pi em um sistema com triplas ligações, faz-se necessária a utilização de outro parâmetro, que é a polarização quadrupolar da ligação C-C tripla. 3.11. Classificação das interações químicas por QTAIM O sinal do Laplaciano da densidade de carga no ponto crítico da ligação (∇2ρb) é o responsável por determinar a curvatura que exerce domínio na região de ligação descrita pelo ponto crítico de ligação. Uma interação compartilhada entre dois núcleos ocorre quando ∇2ρb<0 e ρb apresenta altos valores (10-1 ua.).[94] Além disso, em uma interação compartilhada |λ1|/λ3>1 e G(rc)/ρ(rc) <1. Já uma interação de camada fechada ocorre quando ∇2ρb>0 e ρb apresenta baixos valores (10-2 ua)[83,93,95]. Em uma interação de camada fechada, |λ1|/λ3<1 e G(rc)/ρ(rc) >1. Vale ainda destacar que, nas interações de van der Waals, as ligações de hidrogênio e as ligações iônicas são de camada fechada, como se pode visualizar na Tabela 3. Tabela 3: Valores de (ρb), (∇2ρb), |λ1|/λ3, G(rc)/ρ(rc) de ligações covalentes C-C no etano, benzeno, etileno, das ligações iônicas em LiCl e NaF e do argônio diatômico. As estruturas foram otimizadas via MP2/aug-cc-pVTZ e a função de onda foi obtida via CCSD(T)/aug-cc-pVTZ. Molécula e interação ρb ∇2ρb ε |λ1|/λ3 G(rc)/ρ(rc) DI Ligação C-C no etano 0,2575 -0,7439 0,0000 1,8101 0,2057 0,9965 Ligação C-C no benzeno 0,3321 -1,2003 0,2144 4,1573 0,3004 1,3892 Ligação C-C no etileno 0,3660 -1,3962 0,3897 8,8927 0,3834 1,8863 LiCl 0,0430 0,2619 0,0000 0,1666 1,4274 0,1663 NaF 0,0436 0,3851 0,0000 0,1336 1,9287 0,1671 Ar2 0,0029 0,0118 0,0000 0,1158 0,7515 0,0290 43 Sistemas moleculares estáveis não apresentam modificações estruturais ocasionadas por vibrações moderadas, uma vez que essas vibrações alteram a configuração nuclear, mas não a conectividade dos átomos. Há, entretanto, sistemas instáveis em que ocorrem mudanças nas suas estruturas, que podem ser analisadas por meio da QTAIM.[96] Um tipo de estrutura instável é a bifurcação, que existe apenas para uma configuração particular dos núcleos ao longo do caminho dissociativo. Tomando como base uma estrutura estável de formaldeído e efetuando a separação dos fragmentos CO e HH de forma progressiva, ocorre a formação de uma estrutura de bifurcação entre os átomos de hidrogênio por meio do aparecimento de um ponto crítico ao longo da ligação química. Esse ponto crítico trata-se de um ponto de inflexão, sendo que a derivada segunda é zero, levando à formação de um ponto crítico degenerado (ω=2). Há também um outro tipo de sistema molecular instável, conhecido como estrutura de conflito,[83,96] no qual aparece um ponto crítico entre um átomo e outro ponto crítico. Estruturas de conflito são características de interações fracas com separações interatômicas maiores que o normal. A elipsidade no ponto crítico degenerado de uma estrutura de bifurcação (e também no ponto crítico da estrutura em conflito) apresenta valor elevado (maior do que 5). Nesse contexto, sistemas moleculares, como o ciclopropano, que apresentam curtas distâncias entre RCP e BCP, são mais instáveis, uma vez que ocorre o aumento da elipsidade do BCP. Na expansão em série da energia de interação eletrostática entre as distribuições de carga surge um novo termo, denominado “momento multipolo”. Diversos sistemas de eixos centrados no núcleo apresentam medidas específicas de distorção da densidade eletrônica que são realizadas através dos momentos multipolos. Quanto à definição de um momento monopolo atômico de distribuição de carga M0(Ω), ela é obtida pela análise da sua população eletrônica N(Ω), conforme equação a seguir: ™0(Ω) = −›œž(Ω)Ω = −K(Ω) (4W). Já a carga atômica q(Ω) é dada pela soma de M0(Ω) e Z(Ω) (momento monopolo nuclear, ou seja, número atômico do núcleo), segundo esta expressão: 44 9(Ω) = ™0(Ω) L ]Ω (4u). A carga atômica fornecida pela QTAIM apresenta a vantagem de não depender tanto do conjunto de bases utilizado, quando comparada com outros métodos que separam no espaço de Hilbert as cargas do conjunto de funções de bases. Em se tratando do momento dipolo atômico M1(Ω), ele é o responsável por indicar a distorção da densidade de carga de um átomo, sendo definido por: ™(Ω) = −›œž(r)Ω .Ω (4d), em que M1(Ω ) é uma vetor de três componentes, onde rΩ é o vetor com centro no núcleo do átomo. 3.12. Os índices de localização (LI) e de deslocalização (DI) Na teoria quântica de átomos em moléculas, tanto o índice de localização (LI) quanto o índice de deslocalização (DI) estão associados de maneira intrínseca ao buraco de Fermi e são obtidos através da integração da densidade do referido buraco. [97,98] O índice de localização fornece o número de elétrons localizados em um átomo. Ao consideramos o conjunto de bacias atômicas (ΩA) oriundas da partição do espaço real, o índice de localização qualitativo λ(A) é obtido pela expressão a seguir: [99,100] λ(A) = − ›¢Tn(., .)œ.œ. (@e)Ω£ , em que as informações referentes à interação quântica entre os elétrons e a medida da correlação eletrônica é fornecida por PXC(r1,r2). O índice de deslocalização fornece o número de elétrons entre dois átomos. Para um par de bacias atômicas, o índice de deslocalização δ(A,B) é calculado pela integração de duas coordenadas eletrônicas de PXC(r1,r2) nas bacias ΩA e ΩB, conforme esta expressão:[97,98] 45 ¤((, U) = = −2 › ›¢Tn(., .)œ.œ.Ω¥ Ω£ (@), . A obtenção do número total de elétrons, Nt, dá-se por:[97,98] K¦ = !λ(A)Ω L 12!δ(A, B)Ω (@&), em que a divisão por dois no somatório de DI ocorre devido à normalização da densidade de troca-correlação.[97,98] 3.13. Topologia do Laplaciano da densidade de carga O Laplaciano da densidade de carga é de fundamental importância na abordagem da QTAIM, uma vez que a integração do Laplaciano na região de uma bacia atômica apresenta valor nulo. Sendo assim, o Laplaciano,∇2ρ(r), em que r é a distância radial a partir do núcleo de um átomo, assume não só valores positivos (dispersão de densidade de carga), como também negativos (concentração de densidade de carga) nas regiões do átomo ao longo da linha radial. O ponto onde ∇2ρ(r)=0 é denominado nó esférico e apresenta relação com formação das camadas do átomo. Na curva de ∇2ρ(r) há regiões de concentração de carga (CC) e regiões de dispersão de carga (CD). Essas regiões ocorrem onde ∇2ρ(r)<0 e ∇2ρ(r)>0, respectivamente. A zona de concentração de carga na camada de caroço (CSCC) e a zona da dispersão de carga na camada de caroço (CSCD) correspondem ao CC e CD internos. Já a zona de concentração de carga na camada de valência (VSCC) e a zona de dispersão de carga na camada de valência (VSCD) correspondem ao CC e CD externos. Destaque-se que as regiões CC e CD descrevem cada camada do número quântico principal de um átomo.[101] A topologia do Laplaciano de uma molécula apresenta pontos críticos que são ligados via caminho do gradiente e que surgem quando ∇(∇2ρ) = 0 e a Hessiana de ∇2ρ apresentam 46 curvaturas fundamentais a partir dos autovalores no ponto crítico. O máximo local é dado pelo ponto crítico (3,-3), em que a derivada segunda da Laplaciana é menor do que zero nos três principais eixos. Já um mínimo ocorre quando a derivada segunda da Laplaciana é maior do que zero, sendo representado pelo ponto crítico (3,+3). O mapeamento dos pares eletrônicos de Lewis pode ser obtido através da análise dos máximos locais do VSCC. Figura 3: Representação do VSCC da molécula de Amônia e dos pontos críticos (3,-3) e (3,-1) da molécula de água. 47 Capítulo 4 – Resumo estendido publicado como Artigo Geral na Revista Processos Químicos (2015) - Qualis C 48 Classic Study on the Growth of Metal Clusters AlxNaN-x (3≤ N ≤ 55) Using Genetic Algorithms Modeled by the Gupta Potential Acassio R. Santos a (PG), Breno R. L. Galvãoc (PQ), Caio L. Firmea (PQ), Fabrício de L. Fariasa (IC) and Jadson C. Belchiorb (PQ) a Chemistry Institute, Universidade Federal do Rio Grande do Norte, Av. Senador Salgado Filho, 3000, Lagoa Nova, (59.072-970), Natal, Rio Grande do Norte, Brazil. E-mail:acassioroch@hotmail.com b Chemistry Department, Universidade Federal de Minas Gerais, Av. Antonio Carlos 6627, Pampulha, (31.270-901), Belo Horizonte, Minas Gerais, Brazil c Chemistry Department, Centro Federal de Educação Tecnológica de Minas Gerais, CEFET- MG, Av. Amazonas, 5253 (30.421-169), Belo Horizonte, Minas Gerais, Brazil Keywords: Genetic Algorithm, Gupta Potential, Nanoalloy, Clusters. INTRODUCTION As a new area of research, nanoscience has been in development since the 1980s. Its main scope is the study of nanoscopic phenomena (1nm = 10-9m), where clusters refer to a new type of materials (nanoparticles), in the order of 2 to 10n (n = 6 or 7) atoms or molecules. These particles may be identical or related to two or more distinct species: monoatomic and monomolecular clusters, polyatomic and/or polymolecular clusters, or molecular clusters. Clusters consisting of a small amount of atoms (2-1000 atoms) are denominated nanoclusters. Such systems are adequate models of researching by the means of theoretical methodologies, as they have a smaller number of particles. On the other hand, they may present difficulties regarding experimental studies, such as structural determination.[1] Theoretical studies concerning clusters are relevant and contribute to the interpretation of experimental measures. Among these studies, we may find those that determine the structures of atomic and/or molecular clusters, possibly using mathematical tools related to spatial geometry applied according to theoretical models. Besides, the quantitative study of 49 potential energy provides us with the basis for the understanding of chemical and physical phenomena of atomic and molecular clusters.[2] Particularly, the determination of isomers of greater cluster stability has ample relevance to the development of new materials, where the more stable structure of cluster corresponds to that of lower energy. In order to determine the cluster structure and energies, one should calculate the conformation of lower energy on the potential energy surface (PES), i. e. the global minimum (GM). Generally, the determination of a minimum is considered an optimization problem. The sheer number of isomers with diverse structures constitute one of the chief difficulties in the determination of the GM in PESs, where it may be verified that there is an almost exponential growth in the number of minimum structures related to the number of particles, i. e., the system size.[3,4] Therefore, optimization methods have been developed in order to contribute to the above-mentioned difficulties, whether they are biased or non-biased. The former utilize chemical or physical information regarding the material and initial conditions such as starting geometry, whilst the latter, in contrast, does not use previous information, and the initial geometry is instead defined at random and thus optimized until a minimum of lower energy is found. This current study proposes the application of non-biased methods, since they may reduce the occurrence of problems such as the determination of local minima, which are recurring in other methods due to the imposition of certain types of geometries or of a PES region. The non-biased methods of optimization have two instruments that may be considered the most effective overall: 1) the exhaustive utilization of local minimization, and 2) the discovered structures of local minima, used as initial conditions to the comparison and selection of new structures to be minimized. Among these methods, the Basin-Hopping [5,6] and the Genetic Algorithm (GA) [7] have been the most promising, [5, 6, 8, 9] incorporating the two above-mentioned instruments. Within this work, we will focus on the applications of GA in the research of nanoalloys of Al and Na. In traditional non-biased methods, such as the Monte Carlo and the Molecular Dynamics Simulated Annealing approaches, the difficulties on finding the global minima (GMs) remain. [10, 11] The interest in the research of metal clusters arises mostly from the possibility of creating new alloys from materials in nanoscale, the so-called “nanoalloys”. The theoretical research of nanoalloys has an important role in materials science, and some of its most relevant objectives are: the prediction of stability in structures, their growth, and assistance in the interpretation of spectroscopic and other experimental measures. In this context, it should 50 be highlighted that a great number of methods were reported in the last few years regarding the effective global optimization of atomic and molecular systems. METHODOLOGY Inspired by Darwin’s theory of evolution by natural selection, the genetic algorithms (GA) concern the field of artificial intelligence, and, more specifically, evolutionary computation. It consists in a technique that solves a problem by using an evolutionary process where the surviving solution is, consequently, the best adapted, i. e., developed. We may compare the GA method with the cells of a living organism – DNA blocks – located in specific positions and components of the same group of chromosomes. These blocks codify certain proteins as well as a given feature of the organism, such as skin color. The parental gene pool crossover is to combine during the reproductive process, resulting in a new chromosome that may mutate due to alterations in the DNA, which are usually caused by some factors such as failure in the copies from the parental gene pool. Biologically formed, the offspring adapts or not to its environment, where survival is a measure to its adaptation [12] . The mathematical modeling of this analogy consists in the application of the GA technique to the resolution of a problem by searching a solution that may come up among a variety of other possible solutions and be considered the best option, according to the marking of possible solutions in conformity with its adequacy to the intended answer. In this context, the GA have often been employed in the determination of maxima to the adequacy values regarding the potential energy surfaces associated to the cluster structure. [13,14] Once the GA standard method is applied, the use of evolutionary genetic operators (OP) becomes quite common, bearing resemblance to the crossover operators among individuals (also known as the mating operator), as well as to the Mutation and Natural Selection operators, used towards a test population with the purpose of optimizing, for instance, the cluster structure. The determination of the best-adapted individuals depends on the “quality” assessment, i. e., on the adequacy of each individual that contributes to the system. In the GA model, such assessment is obtained according to the Fitness Function (FF) [15] that, in the context of the clusters, is used to assess the respective energies. The solutions found in the use of the observed technique (GA) are usually considered adequate, since it is not always possible to find or prove which solution is optimal. The GA 51 used in this research has two other evolutionary operators that distinguish them from the standard GA: the Annihilator (AN) and the History (H). [16,17] Considering the composition and size of the cluster (total amount of atoms), its population is generated at random, given that the spatial coordinates from each atom in each cluster are generated at random inside a search space previously defined. After this stage, the GA method uses a list of atomic Cartesian coordinates from each cluster, with the purpose of generating their respective genes, which evolve to the next generation by the application of genetic evolutionary operators from the standard method of Crossover, Mutation and Natural Selection. In this context, the genetic operators are applied to the current population, immediately after each individual of this population is relocated to its nearest local minimum according to the quasi-Newtonian standard method (BFGS). [18] In the current research, the parental individuals are chosen to the application of the Crossover OP through the roulette-wheel selection method. Therefore, the probability of a choice depends on the adequacy of the individual. Once this process occurs, the use of the Mutation OP may or may not be put into effect, considering that this operator is slow and less frequent, so that its application may occur strictly to a fraction of a population from each generation. This does in no way differ from the GA standard method, and yet, in the Mutation OP, the AN removes both offspring individuals and mutants yielding the same energy when compared to already existent individuals from the same generation, while H “memorizes” information concerning the population of that generation for a next use. Hence, the remaining individuals are classified in a list and arranged according to their FF values, where the last- placed in the list are discarded, so it may persist with a constant number of clusters within the population (Natural Selection). [19] The described procedure is a loop from the standard GA (except for a minor interference of AN and H) and is repeated in a determined number of times until convergence, in which it is not possible to produce any more energy. To the GA, henceforth, the global minimum has been found. However, since it is a stochastic process, it is possible that the true global minimum may not be found. For this reason, the better improvement of the modified GA in this research lies in the realization of new cycles, starting from the point where the AN promotes a mass extinction of the system optimized by the standard method. In this case, the procedure reboots and the initial population is no longer generated at random, but provided by H. This new complete loop may now be repeated until, after a determined number of times, no isomer with energy lower than the last more stable obtained is found. [17] 52 Figure 1: Flow chart indicating the operation of the GA used in this research. The traced box on the left corresponds to the operators of the standard method, while the traced ellipsis to the right indicates the new operators. Thus, considering the basic functioning of the GA, it is pointless to target a problem regarding clusters with overly strict ab initio treatments, such as the DFT method (Density Functional Theory), due to the great quantity of local minimizations that the optimization algorithm is supposed to perform. For this reason, an analytical description, experimentally parameterized to the potential energy among the atoms of the system (demanding much simpler computational requirements), is the most reasonable option. Expressions used in potential energy associated to a many-body system are denominated “empirical many-body potentials”, widely recognized for being capable of modeling in close approximation the bonds between metal atoms, as well as of reproducing accurately the thermodynamic and structural properties, mainly from most of the transition metals. [17,20] Currently, it is recognized that the empirical many-body potential energy expressions are capable of reproducing thermodynamic and structural properties of many transition metals with good accuracy, such as the energy of sublimation, the constant of network, elastic constants and vacancy formation energy. [21,22] This is, to the date, one of the most practical manners of addressing simulations of punctual or extended defects, as much as interface or surface properties in metal alloys. [16] Thus, in the last few years, there has been frequent use of empirical potentials towards the analysis of several problems, in the example of the 53 materials science. The empirical many-body potential energy expressions present a superior efficiency when compared to simpler potential expressions that only add contributions due to pair interaction, from which the latter is usually the most adequate in the reproduction of basic features of metal systems. This occurs because in the former it is included the essential feature of band regarding metal alloys. [20] Such expressions concern the effect of each individual (j) from the system in reference to one of the components (i) of the same system, not inasmuch as it depends solely on the interaction of each pair (ij), but considering also the influence in which the neighbors, so far as to a certain distance from each individual in question (j), exert on the caused effect in the assessed component (i). We must highlight that all these contributions are added to all existent i and j in the system. We have it that the Gupta potential belongs to a class of pair-functional potentials based on the Tight-Binding Second-Moment Model. The Gupta potential may be written in the manner of repulsive pairs Vr and term of many-body coupling Vm. The energy of a given atom i constituted by the interaction with all the neighboring atoms j is, accordingly: (1) where Vr represents pair repulsion among metal ions, and is given by: (2) Subsequently, the many-body term, Vm, represents the cohesion of many bodies given by the valence band, i. e., (3) The parameters A, p, r0, ζ e q are values obtained from adjustments in experimental data of the cohesion energy, parameters of network and elastic constants, independent of the crystalline structure at 0 K, not only to pure metals but also to alloys. The term rij is the distance between atoms i and j; the parameters q and p are dimensionless and determine the 54 extent of the terms of band and repulsion, respectively; whereas A and ζ have energy units and establish the force of these terms. The term r0 is the average distance to the atom’s nearest neighbors in the bulk of the metal. [23] The parameters to homodiatomic interactions (Al-Al, Na-Na) within this research were obtained from the literature. [24,25] Simple arithmetic means were used in order to address the heterodiatomic interactions (Al-Na). [23] In this context, the Gupta potential may be used to model an extensive diversity of solids and clusters, as well as of some alloys. At Table 1 we present the values of the parameters used in this research. Table 1: Empirical parameters of the Gupta potential to the system Al-Na r0 (Ǻ) ζ (eV) p q A (eV) Al-Al 2,8637 1,316 8,612 2,516 0,1221 Na- Na 3,6989 0,2911 10,13 1,30 0,01596 Al-Na 3,2818 0,8035 9,371 1,90 0,06903 The average bond energy is an important quantity that may be associated to the clusters and defined as: [16] (4) The value of Eb must be asymptotically close to the value of the cohesion energy of the bulk related to the atoms that compose the cluster, in the case the clusters are relatively sizable. [16] The second difference of average bond energy among clusters is also a useful quantity to the analyzed structures. [26,27] Its calculation may be performed by the following equation: Through this expression, it becomes possible to assess the energetic behavior of a cluster of N atoms in reference to its “neighbors”, i. e., a cluster of N + 1 and another of N-1. Hence, we may observe stability peaks when ∆2Eb(N) > 0 or instability when ∆2Eb(N) < 0. [16] In spite of being a useful quantity, the efficiency of the second energy difference may be, in some fashion, unreliable, since the stability comparison occurs only to the neighboring clusters. For instance: there might be cases where neighboring clusters may present a very high energetic instability, leading to the demonstration of a stability peak by the cluster of N 55 atoms, not because the structure is stable, but otherwise, because the neighboring clusters were too unstable. In order to analyze the alloy clusters, the excessive energy provides another useful quantity, which may give information regarding the likelihood of the alloy formation in relation to the formation of the corresponding pure cluster [28,29]. In our approach, the excessive energy of an agglomeration of N atoms in a AlxNaN-x alloy is given by: Negative values to Eexc (AxBN-x) favor the formation of the corresponding nanoalloy. In this research, we performed the calculation of the excessive energy to all combinations of AlxNaN-x where 3 ≤ N ≤ 55. RESULTS AND DISCUSSIONS The results obtained regarding the geometries of pure aluminum (Figure 2) and pure sodium (Figure 3) clusters were analogue, and agree with the literature [25,30] concerning systems of this type, where, however, other methods were utilized to deal with the issue. Figure 2: Structures obtained as global minima using the GA regarding pure Aluminum clusters from: a) Al3 to m) Al15. 56 Figura 3: Structures obtained as global minima using the GA regarding pure Sodium clusters from: a) Na3 to m) Na15. The pure aluminum structures present sturdier atom compaction if compared to those of sodium. This may be explained by the differences in the atomic radius between those two elements, which are reflected in the equilibrium of interatomic distances provided by the Gupta potential to distinguished atoms. We have verified that the structures of the pure Al clusters are quite similar to those of pure Na, whereas the most important difference occurs in the clusters of 9 and 11 atoms. We noticed that the dihedral angle in the sodium clusters is sharper when compared to the aluminum clusters. Once we compared the structures obtained by the GA using the many-body Gupta potential to the pure aluminum with recent results from the literature, such as those from Kiohara et al. (2013),[31] we verified that the Gupta potential, despite being a classic method, presented the cluster structure to the systems Al2, Al3, Al6, Al8 and Al9 with striking resemblance to those obtained via ab initio calculations, such as DFT and CCSD(T). Since we did not find aluminum doped with sodium clusters in the literature, we managed to provide a few assumptions and observations concerning the behavior of their structures from the analysis of pure sodium and pure aluminum clusters. The following figures present some results to the second energy difference of the clusters AlxNay regarding different nuclearities. 57 Figure 4– Second energy difference of the clusters aluminum-sodium of 7 atoms. Figure 5– Second energy difference of the clusters aluminum-sodium of 8 atoms. Figure 6– Second energy difference of the clusters aluminum-sodium of 10 atoms. Figure 7– Second energy difference of the clusters aluminum-sodium of 15 atoms. By performing an analysis of the second energy difference to the aluminum-sodium mixtures, we verified that higher stability peaks occurred with a value of x=1, i. e., when we 58 had but a single atom of aluminum in the cluster structure. However, considering a more detailed analysis of the clusters in an alloy, the excessive energy may become a useful quantity, because it provides the likelihood of the bimetallic alloy formation, when compared to the pure corresponding cluster [28,29]. Next, we present the results from the excessive energies (Eexc) concerning varying compositions of clusters of AlxNaN-X (3≤ N ≤ 55). Figure 8– Excessive energy to aluminum-sodium clusters of 7 atoms. Figure 9– Excessive energy to aluminum-sodium clusters of 8 atoms. Figure 10– Excessive energy to aluminum-sodium clusters of 10 atoms. 59 Figure 11– Excessive energy to aluminum-sodium clusters of 15 atoms. Figure 12– Excessive energy to aluminum-sodium clusters of 20 atoms. Figure 13– Excessive energy to aluminum-sodium clusters of 25 atoms. Figure 14– Excessive energy to aluminum-sodium clusters of 30 atoms. 60 Figure 15– Excessive energy to aluminum-sodium clusters of 35 atoms. Figure 16– Excessive energy to aluminum-sodium clusters of 40 atoms. Figure 17– Excessive energy to aluminum-sodium clusters of 45 atoms. Figure 18– Excessive energy to aluminum-sodium clusters of 50 atoms. 61 Figure 19– Excessive energy to aluminum-sodium clusters of 55 atoms. By analyzing the excessive energy graphics, we observe that clusters with 3, 4, 5, 6 and 7 atoms have better stability when the agglomerations are composed of 2 aluminum atoms; in other words, the structures Al2Na1, Al2Na2, Al2Na3, Al2Na4, Al2Na5 are the ones most prone to be formed. When, considering the structures of nuclearities 8, 9, 10, 11, 12 and 13, the composition is altered and then presents stronger trend to form clusters with the composition of 4 aluminum atoms - Al4Na4, Al4Na5, Al4Na6, Al4Na7, Al4Na8 and Al4Na9. Regarding agglomerations of 14 atoms, the most stable structure was Al6Na8, and, to those of 15 atoms, the most stable composition was Al5Na10. To clusters of 20 atoms, the most stable structure presented 5 aluminum atom. Clusters with nuclearity of 25 and 30 presented better stability with 10 atoms of aluminum in the structure. Agglomerations of 35 and 40 atoms presented better stability with 15 aluminum atoms. Finally, structures of 45, 50 and 55 atoms presented themselves more stable with 20 aluminum atoms in the structure of the nanoalloy. In Figure 20 we present alloy structures with nuclearity up to 15 atoms. Figure 20: Structure of the compositions Al2Na1, Al2Na2, Al2Na3, Al2Na4, Al2Na5,Al4Na4, Al4Na5, Al4Na6, Al4Na7, Al4Na8, Al4Na9, Al6Na8 and Al5Na10. 62 We have observed that, regarding sodium structures doped with a single atom of aluminum, the aluminum atoms tend to occupy central positions, whilst sodium atoms tend to occupy peripheral areas within the structure. It is compelling to notice that there is an ordered growth in the structure to the compositions of 4 aluminum atoms (20-f, 20-g, 20-h, 20-i, 20-j and 20-k), in which these atoms form a regular tetrahedron in the center and the sodium atoms agglomerate around those of aluminum, in a very similar fashion to a micelle. We have also verified that the 20-a, 20-b, 20-c, 20-d and 20-e structures presented great similarity to those of pure aluminum and pure sodium with the same nuclearity. With N=8, the structures then started presenting configurations much different from those obtained from pure clusters. CONCLUSIONS In the current research, all clusters to the system AlxNaN-X (3 ≤ N ≤ 55) presented structures that may be considered of lesser energy according to the modified genetic algorithm method, along with the many-body Gupta potential [16,17]. We realized that the pure aluminum structures present sturdier atom compaction when compared to those of sodium, which may be explained by the differences in the atomic radius between sodium and aluminum. Such differences are observed in the equilibrium of interatomic distances provided by the Gupta potential to distinct atoms. Since there is no further research in the literature regarding aluminum-sodium clusters, we started from the obtained structures of pure aluminum and pure sodium in order to predict information regarding to corresponding nanoalloys. Once we considered that the aluminum clusters presented more stability than those of sodium, we focused in the use of pure aluminum structures to serve as a comparative basis for the GA method with the Gupta potential. By confronting the obtained structures using the GA with the many-body Gupta potential to pure aluminum, consonant with recent studies in the literature such as those of Kiohara et al. (2013)[31], we have verified that, despite being a classic method, the Gupta potential presented cluster structures to systems Al2, Al3, Al6, Al8 and Al9 in a very similar manner to those obtained from ab initio calculations, such as the DFT and the CCSD(T). In conformity with what has been explained, it became possible to demonstrate that the genetic algorithms applied to the Gupta potential may be an effective tool in the search for global minima structures. 63 ACKNOWLEDGMENTS Work supported by CNPq, FAPEMIG, FAPERN. [1] Schmid G. Metal Clusters in Chemistry I. Wiley-VCH, Weinheim, 1999. [2] Available at: http://www.tc.bham.ac.uk/~roy/Research/clusters.html. Accessed on July 10, 2015. [3] Berry R. S., Braier P., Hinde R. J. Tumer, G.W.; Johnston, R. L.; WILSON, N. T, J. Chem. Phys. 112, 4773 (2000) P., Israel J. Chem., 30 39, 1990. [4] Hartke B., Chem. Phys. Lett., 258 144, 1996. [5] Wales D. J. e Scheraga H. A., Science, 285 1368, 1999. [6] White R. P. e Mayne H. R., Chem. Phys. Lett., 289 463, 1998. [7] Holland J. H., Sci. Am., 267 66, 1992. [8] Hartke B., Phys. Chem, 214 1251, 2000. [9] White R. P., Niesse J. A. E H. R. Mayne., J. Chem. Phy., 108 2208, 1998. [10] Doye J. P. K., Wales D. J. E Berry R. S., J. Chem. Phys., 103 4234, 1995. [11] Doyce J. P. K. e Wales D. J., J. Chem. Soc. Faraday, Trans., 93 4233, 1997. [12] Disponível em: http://professor.webizu.org/ga/, acessed in May 10, 2014. [13] Hartke, B., Z. Phys. Chem. 214 (2000) 1251. [14] Johnston, R. L.; Roberts, Soft Computing Approaches in Chemistry, edited by H. Cartwright, L. Sztandera (Physica-Verlag, Heidelberg, 2001). [15] Silva, E. S. A.; Duarte, H. A.; Belchior, J. C., Chemical Physics 323 (2006) 553–562. [16] Lordeiro, R. A.; Guimarães, F. F.; Belchior, J. C.; Johnston, R. L., International Journal of Quantum Chemistry 95 (2003) 112. [17] Guimarães, F. F.; Belchior, J. C.; Johnston, R. L.; Roberts, C., Journal of Chemical Physics, 116 nº19 (2002) 8327. [18] Schlegel, H.B., Adv. Chem. Phys. 67 (1987) 249. [19] Rodrigues, D.D.C.; Nascimento, A.M.; Duarte, H.A.; Belchior, J.C., Chemical Physics 349 (2008) 91. [20] Cleri, F.; Rosato, V., Physical Review B, 48 nº1 (1993) 22. [21] Daw, M. S.; Baskes, M. I., Phys. Rev. Lett. 50, 1285 (1983); Phys. Rev. B29, 6443 (1984). [22] Rosato, V.; Guillope, M.; Legrand, B., Philos. Mag. A 59, 321 (1989). 64 [23] Aguado, A.; López, J. M., The Journal of Chemical Physics 133 (2010) 094302-1. [24] Li, Y.; Blaisten-Barojas, E., Physical Review B 57 nº24 (1998) 15 519. [25] Tumer, G.W.; Johnston, R. L.; WILSON, N. T, J. Chem. Phys. 112, 4773 (2000) [26] Wilson, N. T.; Johnston, R. L., Eur. Phys. J D (2000), 12, 161. [27] Rao, B. K.; Jena, P., J. Chem. Phys. (1999) 111 1890. [28] A. Aguado and J. M. L´opez, J. Chem. Phys. 133 (2010) 094302 [29] R. Ferrando, J. Jellinek and R. L. Johnston, Chem. Rev. (Washington, D.C.) 108 (2008) 845. [30] Lai, S. K.; Hsu, P. J.; Wu, K. L.; Liu, W. K.; Iwamatsu, M., Journal of Chemical Physics 117 nº23 (2002) 10715. [31] Valéria O. Kiohara; Edson F.V. Carvalho; Carlos W.A. Paschoal; Francisco B. C. Machado; Orlando Roberto-Neto. Chem. Physics Letters 568–569 (2013) 42–48. 65 Capítulo 5 – Artigo em desenvolvimento para submissão em uma revista com Qualis B2, a ser selecionada 66 Análise dos modos de crescimento de clusters de alumínio-lítio e alumínio- potássio de 7 a 55 átomos pela abordagem do Algoritmo Genético Acassio Rocha Santos a (PG), Caio Lima Firmea (PQ), e Jadson Cláudio Belchiorb (PQ) a Instituto de Química, Universidade Federal do Rio Grande do Norte, Av. Senador Salgado Filho, 3000, Lagoa Nova, (59.072-970), Natal, Rio Grande do Norte, Brazil b Departamento de Química, Universidade Federal de Minas Gerais, Av. Antonio Carlos 6627, Pampulha, (31.270-901), Belo Horizonte, Minas Gerais, Brazil Resumo O presente trabalho consiste em um estudo teórico de clusters metálicos AlxLiy e AlxKy (x + y ≤ 55), otimizados por meio da técnica de inteligência artificial Algoritmo Genético (GA). Com o objetivo de tornar o GA mais eficiente, foram aplicados dois novos operadores em relação ao GA padrão: o OP História e o OP Aniquilação. Constatou-se que as estruturas com 6 e 7 átomos de alumínio em suas composições apresentaram geometria octaédrica e bipirâmide pentagonal, respectivamente. Para os clusters de alumínio-lítio, verificou-se que as menores energias de excesso ao longo do crescimento do sistema ocorreram quando a proporção de átomos de alumínio na estrutura foi 30% a 40% Para os aglomerados de alumínio-potássio, essa proporção foi de 40% a 50%. Palavras-chave: Algoritmo Genético. Potencial Gupta. Nanoligas. Clusters. Introdução Na ciência dos materiais, a gama de propriedades de sistemas metálicos pode ser grandemente ampliada, possuindo misturas de elementos para gerar compostos intermetálicos e ligas. Em muitos casos, os efeitos sinérgicos geram uma melhoria nas propriedades específicas sobre ligas; desse modo, a rica diversidade de composições, estruturas e propriedades das ligas metálicas levaram a significativas aplicações em eletrônica, engenharia e catálise.[1] 67 No contexto de estudos de sistemas metálicos, definem-se como clusters o agrupamento de partículas em escala nanométrica, cujas propriedades variam: tanto as que são apresentadas pela sua composição fundamental (átomos e moléculas); quanto as que são apresentadas no bulk, em que o número de partículas formadas torna-se grande o suficiente, podendo haver mudanças nas propriedades do material. De um modo geral, esperam-se que os clusters apresentem propriedades distintas da matéria condensada, dependendo de seus tamanhos, conformações, estrutura eletrônica e energias de ligação.[2] É importante destacar que o tamanho dos clusters pode alterar significativamente sua forma geométrica e sua estabilidade energética. Sabe-se, por exemplo, que os aglomerados de metais alcalinos, com tamanhos de até milhares de átomos (e clusters também menores de Cu e Ag), em conformidade com o modelo jellium, são relativamente estáveis em certas nuclearidades (os chamados números mágicos).[3] Por outro lado, aglomerados de metais de transição e alguns metais do grupo principal (como Al, Ca e Sr) geralmente possuem números mágicos que correspondem aos grupos constituídos por poliedros concêntricos de átomos. Como os clusters tratam-se de um sitema finito e possuem uma elevada razão entre a superfície e o volume, os efeitos relacionados à enegia de superfície são bastante relevantes, podendo acarretar a formação de estruturas não cristalinas[4] Nesse contexto, verifica-se que a composição química pode alterar tanto as propriedades físicas dos clusters quanto suas conformações estruturais,[5] pois esses agregados podem ser homogêneos ou heterogêneos em relação às partículas que o formam. Os clusters atômicos e moleculares têm dependência direta com os átomos que os constituem, sendo os potenciais empíricos usados de forma eficiente para a sua modelagem. [6,7] Nesse modelo, é preciso determinar o mínimo global da superfície de energia potencial (SEP) para assim realizar de uma abordagem semi-quantitativa.[8] Entretanto, uma das principais dificuldades quanto à determinação do mínimo global na SEP refere-se à grande quantidade de isômeros com estruturas distintas: o crescimento do número de estruturas de mínimo torna-se quase exponencial em relação ao número de partículas que compõem o sistema. [9,10] Uma alternativa eficaz para a resolução desse problema é a utilização de um potencial empírico analítico específico, em conjunto com um método de otimização eficaz. Em seguida, essa metodologia pode ser validada comparando-se os resultados obtidos com aqueles do melhor método de estrutura eletrônica disponível. Dessa forma, o potencial empírico Gupta[11,12] tem sido aplicado de modo satisfatório no modelo de sistemas metálicos; e também construído considerando-se o caráter essencial da banda de 68 ligações metálicas, por meio da aproximação de segundo momento (SMA) para o modelo tight-binding de ligação química.[13-15] Além do potencial empírico adequado, faz-se necessária também a utilização de um método de minimização de energia calculada para a expressão analítica. Com o fim de encontrar solução para os problemas de mínimos locais, diversos métodos têm sido desenvolvidos, tendo o Basin-Hopping[16,17] e o Algoritmo Genético (GA)[18] se mostrado mais eficientes na determinação de mínimos globais.[16,17,19,20] Inclusive, algoritmos eficientes são frequentemente aplicados na busca estruturas de mínimo global para diversos tipos de sistemas, como clusters, cristais e biomoléculas.[21] Há crescente interesse no estudo teórico de clusters metálicos, em razão do potencial de aplicações em áreas como catálise e a nanoeletrônica (por exemplo: em dispositivos de tunelamento de um único elétron).[22] Tal estudo tem como objetivos prever a estabilidade das estruturas, analisar seus modos de crescimento e auxilar na previsão de medidas experimentais, como as espectroscópicas. No presente trabalho, serão utilizados o método do GA acoplado ao potencial Gupta, na análise de clusters de alumínio dopados potássio e lítio. Metodologia Ao se utilizar a técnica GA, a finalidade é encontrar a melhor solução dentro de um conjunto de todas as soluções possíveis, conforme a adequação à resposta procurada. Nesse contexto, o GA é frequentemente aplicado para determinar valores que correspondem aos mínimos da SEP associadas aos clusters. [23,24] Para tanto, ocorre a utilização de operadores genéticos (OP) evolucionários que constituem o GA, como o operador de Cruzamento (ou acasalamento) entre as espécies, de Mutação e de Seleção Natural. Esses operadores são aplicados com o objetivo de otimizar o sistema de uma população teste. O ranqueamento dos indivíduos mais bem adaptados é realizado por meio da adequação deles à Fitness Function (Fi),[25] a seguir:  =  1 − ℎ 2  á  − 1 (). O conjunto de soluções geradas pelo GA são consideradas boas, porém, nem sempre é possível provar ou encontrar uma solução considerada ótima. Por essa razão, a fim de refinar 69 a análise da busca por estruturas de clusters por meio da minimização da energia, bem como aumentar a probabilidade de encontrar estruturas de mínimo global na SEP, utiliza-se neste trabalho o GA modificado, o qual possui dois novos operadores: OP Aniquilador (AN) e o OP História (H), distintos do GA padrão. [8,26] Considerando-se o tamanho e a composição do cluster, a população é obtida de forma aleatória; e as coordenadas espaciais geradas randomicamente dentro de um determinado intervalo. Em seguida, com o propósito de criar novos genes, o GA faz uso de uma relação de coordenadas cartesianas de cada átomo do cluster.[27] Esses genes evoluem para a próxima geração através da aplicação dos operadores evolucionários de Cruzamento, Mutação e Seleção Natural. Encerrada essa etapa, cada indivíduo da população é deslocado para o seu mínimo local por meio do método BFGS e, em seguida, os operadores genéticos agem sobre a população atual.[28] Na abordagem do GA, a depender de sua adequação ao sistema, os indivíduos pais são probabilisticamente escolhidos para dar início à aplicação do OP Cruzamento, por meio do método roleta. Após esse processo, o OP Mutação pode ou não atuar, uma vez que a frequência desse operador é pequena e ele é bastante lento; consequentemente, sua aplicação ocorre apenas em uma fração de cada população gerada. Nesse sentido, para os sistemas estudados, a fração escolhida foi de 15%. Ao que parece, não há nenhuma diferença entre esse método e o método GA padrão, todavia, realizada a ação do OP Mutação, o OP Aniquilador faz a remoção de indivíduos e mutantes com energia igual aos membros já existentes na geração; já o OP História “guarda na memória” informações sobre a população dessa geração para uso em uma etapa futura. Os demais indivíduos que sobraram após a aplicação dos operadores mencionados são ranqueados de acordo com seu valor de Fitness Function. Ao fim desse processo, o OP Seleção Natural faz o descarte dos últimos indivíduos da lista para que a população de clusters permaneça constante. [22] Exceto pela ocorrência de pequena inclusão dos operadores Aniquilador e História, o procedimento do GA modificado é o mesmo do GA padrão, o qual constitui um loop que se repete por uma certa quantidade de vezes até atingir a convergência, quando não se conseguem estruturas com energia menor do que as existentes. Nessas circunstâncias, o GA padrão fornece o resultado do mínimo global encontrado, porém isso pode não ser realmente verídico, razão pela qual se faz necessário o aprimoramento do GA modificado, definido como realização de novos ciclos que ocorrem a partir do ponto em que o OP Aniquilador 70 realiza a extinção de todos os indivíduos do sistema refinado pelo método padrão. O procedimento se reinicia, mas a população inicial não é mais totalmente aleatória, e sim promovida pelo OP História. Em seguida, um novo loop completo é repetido até que, após um certo número de vezes, nenhum isômero com energia menor que o último mais estável seja obtido, momento em que o processo evolucionário é interrompido.[8] Portanto, a abordagem do GA modificado tem apresentado melhores resultados para clusters de água, nanoligas de ouro-cobre e de sódio-potássio.[8,29,30] Vale ressaltar ainda a aplicação de potenciais que reproduzem de forma satisfatória as interações interatômicas. Dessa forma, potenciais empíricos de muitos corpos reproduzem com boa precisão diversas propriedades de boa parte dos materiais de transição, tanto termodinâmicas quanto estruturais. Entre essas propriedades, destacam-se a energia de sublimação, a constante de rede, as constantes elásticas e a energia de formação de vacância.[31,32] Como a reprodução de características básicas de sistemas metálicos envolve a inclusão do caráter essencial da banda das ligações metálicas, o uso do potencial empírico de muitos corpos é mais apropriado. [6] No presente trabalho, adota-se o modelo de potencial de Gupta,[32] da energia de coesão e das propriedades físico-químicas dos materiais. O potencial Gupta está incluído no grupo de potenciais funcionais de pares e baseia-se no modelo do segundo momento de Tight (Tigth-Binding Second-Moment Model). Esse potencial é dado em função dos termos dos pares repulsivos Vr e do termo de acoplamento de muitos-corpos Vm, sendo que a energia de um determinado átomo i, calculada pela interação com todos os átomos vizinhos j, é escrita por:  = ∑ "#( ) − $( )% (&)' , em que Vr descreve a repulsão de pares entre íons metálicos e é definido por: .( ) = !()*+,−+ -.g©.0 − 112 ©≠g (4). O termo Vm descreve a coesão de muitos corpos da banda de valência, sendo dado por: 71 Z( ) = «!¬2)*+,−29 -.g©.0 − 112 ©≠g ­ 12 (@). Os parâmetros A, ρ, r0, ® e q foram obtidos através de dados experimentais de energia de coesão, de parâmetros de rede, de constantes elásticas independentes, da estrutura cristalina a 0 K. Ressalte-se que esses dados foram obtidos para os metais puros e também para as ligas formadas. O termo rij refere-se à distância entre os átomos i e j; já os parâmetros q e p são adimensionais, e descrevem o alcance dos tempos de banda e repulsão, respectivamente. Os parâmetros A e ζ apresentam unidades de energia e determinam a força de Vr e Vm. O termo r0 descreve a distância média dos vizinhos que possuem maior proximidade do átomo no bulk do metal. [33] Os parâmetros para interações homoatômicas Al-Al, Li-Li, e K-K, presentes neste trabalho, foram obtidos da literatura. [34-36] Para as interações heteroatômicas Al-Li, Al-K, foram utilizadas médias aritméticas simples,[33] Desse modo, o potencial Gupta pode ser utilizado para a modelagem de diversos sistemas sólidos e clusters, bem como de ligas metálicas. Na Tabela 1, são explicitados os valores dos parâmetros do potencial Gupta utilizados neste trabalho. Tabela 1: Parâmetros de potencial Gupta para os sistemas Al-Al, Li-Li, K-K, Al-Li e Al-K. Interação r0 (Å) ξ (eV) p Q A (eV) Al-Al 2,8637 1,316 8,612 2,516 0,1221 Li-Li 2,9052 0,325 7,75 0,737 0,03327 K-K 4,3673 0,2626 9,596 1,34 0,02054 Al-Li 2,88445 0,8205 8,181 1,6265 0,07768 Al-K 3,6155 0,7893 9,596 1,928 0,07132 Resultados e Discussões O método do GA com potencial Gupta apresentou resultados semelhantes aos cálculos ab initio para os clusters de lítio puro (Li5,Li6, Li7)[37] e alumínio puro (Al6, Al8 e Al9).[38] Porém, no estudo de ligas metálicas, como no caso de clusters Na-K, [30] o GA apresentou falhas na obtenção de clusters menores do que 7 átomos. 72 Figura 1: Estruturas de menor energia para os sistemas a) Li5, b) Li6, c) Li7, d) Al6, e) Al8 & f) Al9. Com o objetivo de se obter estruturas mais próximas do mínimo global, serão utilizados como ponto de partida as estruturas com 7 átomos, ampliando a abordagem realizada até 55 átomos. A estabilidade dos clusters pode ser avaliada por diversas grandezas, entre as quais uma de grande aplicabilidade para a pesquisa é a energia média de ligação, dada por: B = −¯’‚(K)K (I). Caso os clusters possuam tamanho considerável, o valor de Eb assintoticamente do valor da energia de coesão do bulk referente aos átomos constituintes do cluster[26]. Na Tabela 2, são apresentados os valores de energia média de ligação de clusters de alumínio-lítio e alumínio-potássio com apenas um átomo de alumínio na estrutura. Fixou-se o número de átomos de alumínio, a fim de analisar a variação da energia média de ligação à medida que o número de átomos de lítio e potássio foi aumentado. Além disso, utilizam-se resultados de clusters com nuclearidades menores que 7 átomos, com o intuito de obter uma representatividade maior. 73 Tabela 2: Energia média de ligação de alguns clusters de alumínio-lítio e alumínio-potássio dopados com apenas um átomo de lítio na estrutura. Cluster Eb/eV Cluster Eb/eV Al1Li2 0,8872 Al1K2 0,8572 Al1Li3 0,9320 Al1K3 0,8910 Al1Li4 0,9570 Al1K4 0,9065 Al1Li5 0,9776 Al1K5 0,9192 Al1Li6 0,9932 Al1K6 0,9312 Al1Li7 1,0096 Al1K7 0,9225 Al1Li8 1,0260 Al1K8 0,9188 Al1Li9 1,0372 Al1K9 0,9043 Al1Li10 1,0431 Al1K10 0,8836 Al1Li11 1,0432 Al1K11 0,8612 Al1Li12 1,0556 Al1K12 0,8488 Al1Li13 1,0415 Al1K13 0,8296 Al1Li14 1,0391 Al1K14 0,8160 A análise da Tabela 2 nos mostra que, para o sistema alumínio-lítio, o valor de Eb aumenta até o cluster com 13 átomos. No cluster com 14 átomos ou mais ocorre o inverso: o valor de Eb diminui à medida que a nuclearidade aumenta. Para o sistema alumínio-potássio, o valor de Eb aumenta até o cluster com 7 átomos; ao passo que no aglomerado com 8 ou mais átomos, o valor de Eb diminui com o aumento do número de átomos do sistema. Outra grandeza aplicada no estudo de clusters de ligas metálicas é a segunda diferença de energia, que consiste em comparar a energia de uma determinada estrutura com a de seus vizinhos[39,40]. Seu cálculo é realizado através da seguinte expressão: ∆2B(K) = 2B(K)− B(K− 1)− B(KL 1) (M). A abordagem da segunda diferença de energianos permite avaliar o comportamento energético de um dado cluster com N átomos em relação a seus “vizinhos”, ou seja, um cluster com N + 1 e outro com N – 1 átomos. Se ∆2B(K) > 0, encontramos picos de estabilidade; mas, se ∆2B(K) < 0, aparecem picos de instabilidade. [26] Na Figura 1, a seguir, são apresentados os gráficos de segunda diferença de energia para com nuclearidades 7, 8, 10, 15 e 20. 74 Figura 1: Gráfico da segunda diferença de energia para clusters de alumínio-lítio e alumínio- potássio com 7, 8, 10,15 e 20 átomos. Com base na Figura 1, podemos observar que as estruturas com maior estabilidade apresentaram composições ricas tanto em lítio para o sistema AlxLiy, quanto em sódio para o sistema AlxKy. No sistema alumínio-lítio e no sistema alumínio-potássio, os átomos de alumínio tendem a ocupar o sítio central das estruturas; os átomos de lítio e potássio, ao contrário, ocupam as partes periféricas dos clusters, como verifica-se nas Figuras 2 e 3, a seguir. Isso também é observado nas Figuras 5 e 6 das energias de excesso, expressas adiante. 75 Figura 2: Clusters com maior ∆BC para o sistema alumínio-lítio: a) Al1Li6; b) Al1Li7; c) Al1Li8; d) Al1Li9; e)Al1Li10; f) Al1Li11; g) Al1K12; h) Al1Li13; i) Al1K14; i) Al2K18 Figura 3: Clusters com maior ∆BC para o sistema alumínio-potássio: a) Al1K6; b) Al1K7; c) Al1K8; d) Al1K9; e)Al1K10; f) Al2K10; g) Al2K11; h) Al2K12; i) Al2K13; i) Al4K16 Verifica-se ainda nas Figuras 2 e 3 que, para o sistema alumínio-lítio, as estruturas com maior estabilidade aconteceram quando o sistema possuía apenas um átomo de alumínio na estrutura. A única exceção foi o sistema com nuclearidade igual a 20, cujo pico de estabilidade foi de 2 átomos de lítio no sistema. Nesse contexto, um dado interessante é o de que as estruturas Al1Li6 e Al1Li7, consideradas como as mais estáveis para o sistema de alumínio-lítio, são bastante semelhantes às obtidas por Tai et al. (2010) por meio de cálculos CCSD(T)[41]. Nesse trabalho dos autores, as estruturas Al1Li6 e Al1Li7 mostraram índice de aromaticidade NICS elevado, o que é característico de sistemas aromáticos, sugerindo assim os picos de estabilidade observados para esses clusters. O sistema alumínio-potássio, apresentou, nos clusters com nuclearidades 7, 8, 9, 10 e 11 átomos, picos de estabilidade com apenas 1 átomo de alumínio; nos clusters com nuclearidades 12,13,14 e 15 átomos, picos de 76 estabilidade com 2 átomos de alumínio na estrutura; no cluster com 20 átomos, o pico foi com 4 átomos de alumínio na estrutura. Embora a segunda diferença de energia nos forneça informações quanto à estabilidade energética de determinados clusters, ela apresenta algumas limitações, pois, como tal estabilidade é medida em relação aos clusters vizinhos, falsos picos de estabilidade podem ocorrer caso esses clusters sejam muito instáveis, isto é, possuam alto valor energético. Nesse sentido, têm sido expressiva a quantidade de pesquisas acerca de clusters de ligas metálicas, com foco na energia de excesso, a qual é responsável por informar sobre o quão favorável ou não é a formação da liga quando comparada ao cluster puro correspondente[2,42]. O cálculo da energia de excesso para um cluster de N átomos de uma liga bimetálica AxBN-X é dado por: B)*((VUK−V) = ¯’‚((VUK−V)−V¯’‚((K)K − (K−V)¯’‚(UK)K (W), em que valores negativos para BRS((TUHT) favorecem a formação da nanoliga correspondente. Agrupamentos de gráficos das energias de excesso para clusters de alumínio-lítio e alumínio-potássio são expressos nas Figuras 4 e 5, a seguir. 77 Figura 4: Energia de excesso para clusters de alumínio-lítio e alumínio-potássio com 7, 8, 10, 15 e 20 átomos. 78 Figura 5: Energia de excesso para clusters de alumínio-lítio e alumínio-potássio com 25, 30, 35, 40, 45, 50 e 55 átomos. Constatou-se que os clusters de sistema alumínio-lítio tiveram menor energia de excesso para o sistema com 7 átomos quando esse possuí apenas 2 átomos de alumínio na estrutura, isto é, Al2Li5. Em se tratando dos sistemas com 8, 9, 10, 11, 12 e 13 átomos, a composição dos clusters altera-se e passa a apresentar maior favorabilidade de formações com 4 átomos de alumínio, quais sejam: Al4Li4, Al4Li5, Al4Li6, Al4Li7, Al4Li8 e Al4Li9. Já em aglomerados com 14, 15 e 20 átomos, as estruturas mais estáveis foram Al6Li8, Al5Li10 e Al5Li15, respectivamente. Os clusters com 25, 30, 35, 40, 45, 50 e 55 átomos mostraram-se mais favoráveis para se formarem com as estruturas Al10Li15, Al10Li20, Al15Li20, Al15Li25, Al20Li25, Al20Li30 e Al20Li35, respectivamente. Os clusters de alumínio-potássio apresentaram composições diferentes daquelas do sistema de alumínio-lítio para os clusters com 7, 13, 15, 30, 40, 50 e 55 átomos, ou seja, Al3K4, Al6K7, Al6K9, Al15K15, Al20K20, Al25K25, Al25K30. Nas Figuras 6 e 7 é possível observar estruturas com menores energias de excesso para os sistemas alumínio-lítio e 79 alumínio-potássio, respectivamente. Um dado interessante é o de que o valor da energia de excesso para os clusters de alumínio-potássio é sempre menor do que os valores das energias de excesso dos aglomerados de alumínio-lítio. Isso sugere que os clusters de AlxKy são mais estáveis energeticamente que os de AlxLiy. Figura 6: Clusters com composições que mostram um mínimo sobre o excesso de energia : a)Al2Li5, b) Al4Li4, c) Al4Li5, d) Al4Li6, e) Al4Li7, f) Al4Li8, g) Al4Li9, h) Al6Li8, i) Al5Li10, j) Al5Li15, k) Al10Li15, l) Al10Li20, m) Al15Li20, n) Al15Li25, o) Al20Li25, p) Al20Li30 e q) Al20Na35 80 Figura 7: Clusters com composições que mostram um mínimo sobre o excesso de energia : a)Al3K4, b) Al4K4, c) Al4K5, d) Al4K6, e) Al4K7, f) Al4K8, g) Al6K7, h) Al6K8, i) Al6K9, j) Al5K15, k) Al10K15, l) Al15K15, m) Al15K20, n) Al20K20, o) Al20K25, p) Al25K25 e q) Al25K30. As composições com 4 átomos de alumínio na Figura 6 para os clusters de AlxLiy (6-b, 6-c, 6-d, 6-e, 6-f e 6-g), bem como as composições com 4 átomos de alumínio na Figura 7 para os clusters de AlxKy (7-b, 7-c, 7-d, 7-e, 7-f), parecem seguir um crescimento ordenado dos clusters, em que os átomos de alumínio formam um tetraedro regular no centro da estrutura, ao passo que os átomos de lítio ou potássio aglomeram-se ao redor dos átomos de alumínio de forma semelhante a uma micela. Além disso, para estruturas com 6 e 7 átomos de alumínio em suas composições, esses apresentam geometria octaédrica e bipirâmide pentagonal, respectivamente. Esse fenômeno foi observado também no estudo para os clusters de sódio-potássio[30]. Conclusão Nos clusters de alumínio-lítio, verificou-se que a proporção de 30% a 40% de alumínio na estrutura foi a que apresentou menor energia de excesso ao longo do crescimento do sistema. Nos aglomerados de alumínio-potássio, essa proporção foi de 40% a 50%. Estruturas com Li5, Li6 e Li7 obtiveram resultados semelhantes às estruturas encontradas na 81 literatura, em que foram aplicados cálculos CCSD (T); [37] o mesmo ocorreu com as estruturas de Al1Li5, Al1Li7 e Al1Li8 do GA. [38] O GA com o potencial Gupta também foi capaz de prever a segregação das estruturas em que os átomos de alumínio ocupam sítios centrais, enquanto que os átomos de sódio e potássio ficaram segregados na superfície. Por fim, mesmo sendo o GA com potencial Gupta um método clássico, os dados obtidos por meio da sua aplicação foram similares aos da literatura tratados por métodos ab inito. Referências [1] Ferrando, R.; Jellinek, J.; & Johnston, R. L.; Chem. Rev. (Washington, D.C.) 108 (2008) 845. [2] Bordón, L. O. P.; Computational Studies of Transition Metal Nanoalloys, 2011, 173 f. Tese (Doutorado em Química), School of Chemistry, University of Birmingham. Birmingham. 2011. [3] Cox H.; Johnston, R. L.;Murrell J. N. J.; Solid State Chem. 145 517, (1999) [4].Rodrigues, D.D.C.; Naciemento, A.M.; Duarte, H.A.; Belchior, J.C.,Chemical Physics 349 (2008) 91. [5].Baletto, F.; Ferrando, R., Rev. Mod. Phys. 77 (2005) 371. [6] Cleri, F.; Rosato, V., Phisical Review B, 48 nº1 (1993) 22. [7] Jorgensen, W. L., J. Am. Chem. Soc. 103 335 (1981). [8] Guimarães, F. F.; Belchior, J. C.; Johnston, R. L.; Roberts, C., Journal Of Chemical Physics, 116 nº19 (2002) 8327. [9] Berry, R. S.; Braier, P.; Hinder, R. J. e Cheng H. P., Israel J. Chem., 30 39, 1990. [10] Hartke, B.; Chem. Phys. Lett., 258 144, 1996. [11] Gupta, R. P.; Phys Rev.B 23, 6265 (1981). [12] Cleri, F e Rosato, V.; Phys. Rev. B: Condes. Matter Mater. Phys. 48, 22 (1993). [13] Aguado, A. e López, J. M.; J. Chem. Phys. 133, 094302 (2010). [14] Li, Y.; Blaisten-Barojas, E. e Papaconstantopoulos, D. A. Phys. Rev. B: Condes. Matter Mater. Phys. 57, 15549 (1998). [15] Reyes-Nava, J. A.; Garzón, I. L.; Beltrán, M. R. e Michaelian, K. Rev. Mex. Fís. 48, 450 (2002). [16] Wales, D. J. e Scheraga, H. A.; Science, 285 1368, 1999. 82 [17] White, R. P. e Mayne H. R.; Chem. Phys. Lett., 289 463, 1998. [18] Holland, J. H.; Sci. Am., 267 66, 1992. [19] Hartke, B.;Phys. Chem, 214 1251, 2000. [20] White, R. P.; Niesse J. A. e Mayne, H. R.; J. Chem. Phy., 108 2208, 1998. [21] Wales, D. J.; Scheraga, H. A., Science 285 1368 (1999). [22] Rao, B. K.; Jena, P., J. Chem. Phys. (1999) 111 1890. [23] Wales, D. J.; Hodges, M. P., Chem. Phys. Lett. 286 (1998) 65. [24] Roberts, C.; Johnston, R. L.; Wilson, N. T., Theor. Chem. Acc. 104 123 (2000). [25] Silva, E. S. A.; Duarte, H. A.; Belchior, J. C., Chemical Physics 323 (2006) 553–562. [26] Lordeiro, R. A.; Guimarães, F. F.; Belchior, J. C.; Johnston, R. L., International Journal of Quantum Chemistry 95 (2003) 112. [27] Zeiri, Y., Phys. Rev. E 51 (1995) R2769. [28] Schlegel, H.B., Adv. Chem. Phys. 67 (1987) 249. [29] R. A. Lordeiro, F. F. Guimar ̃aes, J. C. Belchior and R. L.Johnston, Int. J. Quantum Chem. 95 (2003) 112. [30] Silva, M. X.; Galvão, B. R.; Belchior, J. C.; J. Mol. Model. (2014) 20: 2421. [31] Daw, M. S.; Baskes, M. I., Phys. Rev. Lett. 50, 1285 (1983); Phys. Rev. B29, 6443 (1984). [32] Rosato, V.; Guillope, M.; Legrand, B., Philos. Mag. A 59, 321 (1989). [33] Aguado, A.; López, J. M., J. Chem. Phys. 133 (2010) 094302-1. [34] Li, Y.; Blaisten-Barojas, E.; Phisical Review B 57 nº24 (1998) 15 519. [35] Tumer, G.W.; Johnston, R. L.; Wilson, N. T; J. Chem. Phys. 112, 4773 (2000) [36] Li, Y. e Blaisten-Barojas, E.; Phys. Rev. B, 57, 24 (1998). [37] Alexandrova, A. N.; Boldiver, A. I.; J. Chem. Theory Comput. 1, 566 (2005). [38] Kiohara, O. V.; Carvalho, E. F. V; Paschoal, C. W. A.; Machado, F. B. C., Chem. Phys. Leters 568-569 42 (2013). [39] Wilson, N. T.; Johnston, R. L., Eur. Phys. J D (2000), 12, 161. [40] Rao, B. K.; Jena, P.; J. Chem. Phys. (1999) 111 1890. [41] Tai, T. B.; Nhat, P. V.; Nguyen, M. T.; Phys. Chem. Chemical Physics, 12, 11477 (2010). [42] Aguado, A. e López, J. M.; J. Chem. Phys. 133 (2010) 094302. 83 Capítulo 6 – Artigo em desenvolvimento para submissão em uma revista com Qualis A1, a ser selecionada 84 Análise estrutural, potencial de ionização, polarizabilidade e topologia de clusters bimetálicos neutros NaxLiy (4 ≤x+y≤ 10) Frederico Texeira Silva* a (PG), Acassio Rocha Santos (PG) b, Caio L. Firmeb (PQ) e Jadson C. Belchiorc (PQ) a Departamento de Química Fundamental, Universidade Federal de Pernambuco, Av. Prof. Moraes Rego, 1235, Cidade Universitária, (50.670-901, Recife, Pernambuco, Brasil b Instituto de Química, Universidade Federal do Rio Grande do Norte, Av. Senador Salgado Filho, 3000, Lagoa Nova, (59.072-970), Natal, Rio Grande do Norte, Brasil. c Departamento de Química, Universidade Federal de Minas Gerais, Av. Antônio Carlos 6627, Pampulha, (31.270-901), Belo Horizonte, Minas Gerais, Brasil Resumo Neste artigo, desenvolvemos um estudo de clusters NaxLiy, com x+y variando de 4 a 10. Desse modo, exploramos a superfície de energia potencial no nível MP2, com o auxílio de um Algoritmo Genético (doravante, GA). As estruturas do lítio puro, encontradas neste trabalho, estão de acordo com a literatura recente. Para as estruturas resultantes, calculamos o potencial de ionização e a polarizabilidade. O erro médio entre essas duas propriedades e os dados experimentais foi de 6% e 13%, respectivamente. Aplicando-se a QTAIM, a análise topológica forneceu importantes informações sobre as interações dos pares atômicos envolvidos, mostrando que, com aumento do tamanho do cluster em relação ao sistema diatômico, ocorre uma diminuição das interações atômicas. O grau de degenerescência do D3BIA e a análise das cargas atômicas mostram a influência (por transferência de carga) do elemento em menor quantidade no cluster em relação aos demais átomos. Palavras-chave: Clusters metálicos. Algoritmo Genético. MP2. Coupled-Cluster. QTAIM. Introdução A tecnologia atual possibilitou a construção das primeiras nanomáquinas.[1–3] Nesse contexto, sabemos que agregados de átomos ou moléculas na escala nanométrica possuem propriedades únicas que dependem da composição e do tamanho, sendo distintas da matéria 85 condensada. Assim, consideramos que compreender essas propriedades pode contribuir para a construção dos novos tipos de aparelhos mencionados. Ligas de sódio-lítio possuem alta miscibilidade e, separando-se apenas na fase líquida a 573K, podem não só solubilizar hidrogênio, como também são promissoras na construção de anodos.[4–6] Entre estudos teóricos e experimentais sobre clusters NaxLiy destacam-se o de Mohajeri et al., [7], de Benichou et al.[8] e de Pérez et al. [9] Mohajeri et al. [7] avaliaram a estrutura dos isômeros para x+y=6. Já Benichou et al.[8] relataram o potencial de ionização para diversas composições. Por fim, Pérez et. al.9 mediram a polarizabilidade de algumas estruturas selecionadas. Mas, até onde temos conhecimento, não há nenhum estudo teórico sobre clusters de nove e dez átomos para o sistema sódio-lítio. A previsão do potencial de ionização é sensível ao modelo teórico utilizado, podendo servir como um benchmarck do método em questão.[10] Cálculos do referido potencial podem auxiliar a interpretação de dados experimentais[11] e descrever semicondutores.[12] Já a polarizabilidade é importante, por exemplo, no estudo de interações de van der Waals[13] e é calculada com boa exatidão para vários sistemas.[14–16] Clusters apresentam um tipo de ligação mais complexo do que moléculas tradicionais. Sendo assim, as técnicas provenientes da intuição química não conseguem descrever, nem mesmo qualitativamente, a geometria desses compostos.[17] Do ponto de vista teórico, encontrar a geometria mais estável significa descobrir a configuração que represente o mínimo global na superfície de energia potencial do sistema. Para isso, o potencial precisa ser escolhido em termos da qualidade e do custo computacional, pois é significativa a quantidade de cálculos para se encontrar o mínimo global, devido ao grande número de variáveis do sistema. Algumas alternativas para descrever a interação entre os átomos são o uso de potenciais empíricos6, da teoria do funcional da densidade[17–21] e do MP2.[22] Nesses estudos, os métodos de otimização global mais recorrentes são o Algoritmo Genético[23] e o Basin Hoping.[24] Considerando-se que o sistema de sódio puro já foi estudado pelo nosso grupo, [6,22] no presente trabalho tratamos acerca do lítio puro e dos clusters bimetálicos sódio-lítio. Com o objetivo de encontrar o mínimo global desses clusters, aplicamos um Algoritmo Genético acoplado a métodos de estrutura eletrônica. De posse das estruturas finais, calculamos o potencial de ionização e a polarizabilidade. Neste artigo, utilizamos o índice de deslocalização (DI), dados topológicos do ponto crítico de ligação, carga atômica QTAIM, energia atômica de acordo com teorema do virial [25] e topologia do Laplaciano para quantificar e qualificar as interações atômicas para os 86 sistemas NaxLiy com x+y=6. Além disso, foram calculados os índices de aromaticidade D3BIA[26] e NICS27 para os sistemas NaxLiy com x+y=6. Método A superfície de energia potencial foi explorada com um Algoritmo Genético (GA) desenvolvido por nosso grupo [22] sendo a população escolhida composta por 40 indivíduos e os clusters iniciais gerados aleatoriamente em uma esfera de raio 2.5 Å em todos os casos. Duas condições de parada foram determinadas: uma delas se o GA atingisse 300 gerações e a outra se o indivíduo com menor energia permanecesse nessa condição por 50 gerações. Além disso, durante as 50 primeiras gerações, o GA parava as otimizações locais caso o gradiente máximo fosse inferior a 0.5 Hartree/Bohr. Entre 50 e 100 gerações, o critério de parada nas otimizações locais foi de 0.05 Hartree/Bohr. Entre 100 e 200 gerações, o critério foi de 0.005 Hartree/Bohr. Posteriormente, critério de parada passou a ser 0.0005 Hartree/Bohr. Todos os cálculos do GA foram repetidos duas vezes com sementes diferentes. Os três resultados apresentaram-se iguais em muitos casos, possuindo um desvio padrão médio de 1kcal/mol. Condições iniciais diferentes levando a um mesmo resultado indicam que os parâmetros da simulação foram adequados a esse sistema. Quanto ao nível teórico, escolhemos para a descrição das interações entre os átomos o MP2/LANL2DZ. Em alguns casos, refinamos esses cálculos nos níveis teóricos B3LYP/6- 311g+(d), MP2/aug-cc-pVTZ, CCSD/aug-cc-pVTZ, CCSD(T)/aug-cc-pVTZ e CCSD(T)/6- 31++g(d). Os cálculos de estrutura eletrônica deste trabalho, realizamos com os programas GAMESS-US[28] e GAUSSIAN 09.[29] Calculamos, por sua vez, o potencial de ionização através do método da excitação vertical. Avaliamos a polarizabilidade com uma diferenciação numérica dupla da energia sob a influência de um campo finito.30 Em todos os casos, atribuímos a multiplicidade como sendo a menor possível. Cálculos de frequência foram realizados garantindo que as estruturas encontradas eram mínimos na superfície de energia potencial. Como exceção, as estruturas de x+y=9 não convergiram a mínimos no nível teórico MP2, por isso, essas foram refinadas no nível B3LYP/6-311g+(d), o qual garantiu a convergência e permitiu o cálculo das propriedades. As estruturas de x+y=4, bem como Na1Li4, Na3Li2, foram reotimizadas no nível CCSD(T)/aug-cc-pVTZ. Além dessas, refinamos as estruturas Na4Li1 e Na1Li5 nos níveis teóricos CCSD(T)/6-31++g(d) e CCSD/aug-cc-pVTZ, respectivamente. Para os cálculos 87 QTAIM, utilizamos os programas AIMALL[31] e AIM2000[32] a partir de função de onda obtida pelo método CCSD(T)/aug-cc-pVTZ. Resultados e Discussão Estrutura, potencial de ionização e polarizabilidade As estruturas dos clusters Li4, Na3Li, Na2Li2 e Na3Li /se assemelham a losangos e a estrutura do Na4Li se assemelha a um trapézio. Comparadas essas estruturas em três níveis de cálculo - MP2/LANL2DZ, MP2/aug-cc-pVTZ e CCSD(T)/aug-cc-pVTZ – pode-se indicar os parâmetros geométricos avaliados, na Figura 1: Figura 1: em a), temos um esquema da estrutura de mínimo global para os clusters de sódio e lítio com composição x+y=4; em b) temos o mínimo global do Na4Li1, sendo o sódio em marrom e o lítio em amarelo. Os dados geométricos para os diferentes métodos estão expressos na Tabela 1. Na Tabela 1, observa-se que o erro médio da geometria com MP2/aug-cc-pVTZ e MP2/LANL2DZ foi de 2,1% e 3,7%, respectivamente, quando comparados com o CCSD. Os resultados oscilaram em torno da referência, ora superestimando a distância, ora a subestimando. Consideramos a diferença pequena e concluímos que o impacto do aumento de base é benéfico, mas não essencial no estudo aqui proposto. 88 Tabela 1: Comparação de características geométricas para os clusters x+y=4 e o Na4Li. As distâncias estão em Angstrons. R e r referem-se, respectivamente, às distâncias maiores e menores do losango formado por essas estruturas, exceto no caso do Na4Li. Nesse caso, R e r referem-se, respectivamente, a base maior e a base menor do trapézio formado por essa estrutura. A variação na geometria foi definida como: ||       . CCSD(T) MP2/LANL2DZ MP2/pVTZ |∆R|-LANL2DZ |∆R|-pVTZ Li4 R 5.41764 5.64995 5.54116 4.3% 2.3% r 2.66838 2.74175 2.6433 2.7% 0.9% Na1Li3 R 5.71014 5.96171 5.80045 4.4% 1.6% r 2.68675 2.7698 2.66817 3.1% 0.7% Na2Li2 R 6.19648 6.14989 6.00096 0.8% 3.2% r 2.98198 3.06444 2.94235 2.8% 1.3% Na3Li1 R 5.9047 6.43843 6.25498 9.0% 5.9% r 2.96844 3.08736 2.96438 4.0% 0.1% Na4Li1 R 6.45914 6.57189 6.39132 1.7% 1.0% r 3.66908 3.81561 3.80027 4.0% 3.6% Após a comparação entre as características geométricas do lítio puro encontrado neste trabalho e um conjunto apresentado na literatura,[23] verificou-se que os clusters obtidos neste trabalho apresentaram, qualitativamente, a mesma geometria da referência. Figura 2: Mínimo global para estruturas de lítio puro (Lin, onde n varia de 5 a 7). Na Tabela 2 são expressas as distâncias entre os átomos de lítio para as três estruturas mencionadas. Observa-se que as diferenças entre as distâncias obtidas neste artigo obtiveram um erro médio nas distâncias de 10% em relação ao método CCSD apresentado na referência. 89 Tabela 2: Propriedades geométricas do Lin, onde n varia de 5 a 7, para as estruturas de mínimo global. As distâncias estão em Angstrons. A variação na geometria foi definida como: ||  |  | ⁄ Li5 B3LYP* CCSD(T)* MP2-ECP |∆R|(CCSD-MP2) R(Li1-Li2) 2.97 2.98 3.10 4% R(Li1-Li3,4) 2.63 2.72 3.21 18% R(Li2-Li3) 3.05 3.07 2.81 8% R(Li3-Li4) 3.02 3.07 3.26 6% Li6 R(Li1-Li2) 3.05 3.08 2.95 4% R(Li2-Li3) 2.80 2.81 3.70 32% Li7 R(Li1-Li2) 2.91 2.97 3.08 4% R(Li2-Li3) 3.02 3.06 3.16 3% * Valores obtidos da referência. [23] Ainda no que concerne ao cluster de lítio puro, obtivemos algumas estruturas diferentes das relatadas na literatura. É importante destacar que há uma diferença de método e base nesses casos: na literatura, foi DFT-B3LYP; neste trabalho, o método utilizado foi MP2. Na Figura 3, verificam-se as comparações de energias no método MP2/LANL2DZ dos respectivos isômeros, construídos com base na literatura.[17,33] Figura 3: Alguns clusters de lítio puro e seus respectivos isômeros mais estáveis. 90 Destaca-se, na Figura 3, o Li5-I, relatado neste trabalho e na referência[23] como mínimo global, mas que sofreu uma inversão da estabilidade no método usado por Centeno et al.[17] A estrutura Li8-II, não relatada por Centeno et al., apresentou baixa energia no método usado por este trabalho e na referência, possivelmente por ter se revelado um isômero instável.[33] Por outro lado, o Li10 obtido discordou da geometria proposta por ambos os trabalhos.[17,33] No método MP2/LANL2DZ, o mínimo do Li10 é 0.04eV, mais estável do que o apresentado por outros trabalhos. Entretanto, evidentemente isso se inverte no nível de cálculo utilizado pela literatura. As demais estruturas de lítio estão plenamente coerentes com os outros trabalhos. Em resumo, o GA utilizado neste trabalho indicou o mínimo global em todos os casos, de modo que as divergências entre os resultados apresentados e a literatura se devem à escolha do nível de cálculo utilizado. Tal resultado indica que o método utilizado neste trabalho permite a realização de um estudo eficiente dos clusters bimetálicos sódio-lítio para essa mesma quantidade de átomos. Segregação Deshpande et. al. [21] analisaram estruturas dos clusters NanLi e NaLin, para 4≤n≤10 no nível DFT/GGA. Esse estudo apresentou resultados diferentes do que constatamos neste trabalho (Figura 4). O Na4Li e Li4Na investigados se apresentaram planares, diferentes dos C2v propostos por Deshpande et al.. [21] Além disso, o NaLi5 teve o átomo de sódio na borda, e não no eixo principal. Por seu turno, o Na7Li apresentou uma geometria diferente, sendo uma variação de um dos isômeros do lítio puro. No nível utilizado por Deshpande et al. [21], o mínimo global do Li9 foi diferente tanto do que encontramos neste artigo quanto de outros.[17,33] Embora a previsão do mínimo global do Li10 tenha sido a mesma que observamos, os clusters bimetálicos apresentaram estruturas diferentes. Em resumo, o método utilizado por Deshpande et al., [21] apresenta interações Li-Na mais fortes, pois o sódio adquiriu posições ligeiramente mais centrais do que as expressas neste artigo. 91 Figura 4: Comparação entre as estruturas Na1Lix e NaxLi1, para 5 ≤ x ≤ 10. Isômeros dos clusters de sódio-lítio com 6 átomos Na série x+y=6, observou-se que as estruturas adquiriram uma configuração planar entre os clusters Na3Li3 e Na4Li2. Por isso, outros isômeros foram construídos para avaliar essa conversão; aqueles mais estáveis para cada composição são representados na Figura 5. Constatamos que, no caso do Na4Li2, a diferença de energia entre o isômero planar e o octaédrico é bem pequena, inclusive, é a menor diferença entre todos os isômeros que calculamos. Isso indica que poderá haver coexistência dessas duas formas nos sistemas que contenham a referida composição. Já no caso do Na5Li1, a diferença entre o primeiro isômero e a geometria mais estável é substancial, indicando que a geometria plana será predominante nesses sistemas. 92 Figura 5: Isômeros dos clusters NaxLiy, para x+y=6. As variações de energia marcadas com asterisco, ∆E*, indicam que essa geometria possui frequências imaginárias. Propriedades O erro médio do potencial de ionização foi de 6%, pouco maior do que os 4% apresentados por Benichou et al. [8] Atribuímos essa diferença ao tamanho dos sistemas estudados e, assim, focalizamos estruturas com clusters entre 4 e 10 átomos. Benichou et al [8] avaliaram os potenciais dos clusters entre 2 e 6 átomos, permitindo que um nível de cálculo superior fosse utilizado. Mesmo assim, trata-se de uma descrição razoável do sistema, que indica a possibilidade de o método realizar a previsão das demais 23 estruturas, cujos dados experimentais estão indisponíveis. Tabela 3: Potencial de ionização e polarizabilidade de algumas estruturas selecionadas de NaxLiy.O potencial de ionização experimental foi medido por Benichou et. al.8 e a polarizabilidade obtidas nas referências.10,20 IP IPexp |IP-IPexp|/IPexp α/n αexp/n |α-αexp|/αexp Na4 3.75 4.18 10% 21.38 21.00 2% Na5 3.63 3.80 4% 21.33 18.67 14% 93 Na6 4.05 4.02 1% 19.94 18.63 7% Na7 3.67 3.76 2% 15.03 17.11 12% Na8 4.08 4.03 1% 17.56 16.80 5% Na9 3.26 3.73 13% 18.69 17.50 7% Na10 3.31 3.83 14% 18.72 19.40 4% Na3Li1 3.83 4.24 10% 19.31 18.90 2% Na4Li1 3.62 3.95 8% 19.98 Na5Li1 4.05 4.15 2% 18.23 Na6Li1 3.73 3.90 4% 15.76 Na7Li1 3.98 4.00 1% 17.01 Na8Li1 3.27 3.70 12% 14.92 Na9Li1 3.49 3.73 7% 17.70 Na2Li2 3.98 4.37 9% 17.31 15.00 15% Na3Li2 3.71 4.14 10% 18.15 Na4Li2 4.09 4.20 3% 17.36 Na5Li2 3.50 3.94 11% 8.75 Na6Li2 3.99 4.08 2% 16.25 19.20 15% Na7Li2 3.34 3.78 12% 14.16 Na8Li2 3.52 3.81 8% 17.04 Na1Li3 4.08 15.55 13.70 13% Na2Li3 3.84 16.35 Na3Li3 4.01 4.16 3% 16.13 Na4Li3 3.55 3.93 10% 11.70 Na5Li3 4.08 4.04 1% 15.41 17.40 11% Na6Li3 3.46 3.65 5% 13.40 Na7Li3 3.56 3.76 5% 16.19 Li4 4.26 4.31 1% 13.89 12.30 13% Li5 3.83 4.02 5% Li6 4.32 4.20 3% 12.56 8.90 41% Li7 3.84 3.94 3% 10.52 11.40 8% Li8 4.55 4.16 9% 11.77 10.40 13% Li9 3.58 3.29 9% 10.15 9.90 3% Li10 4.01 3.95 1% 11.87 10.40 14% Na4Li4 4.15 14.53 16.50 12% Na3Li5 4.31 13.92 14.30 3% Na2Li6 4.36 13.19 14.60 10% Na1Li7 4.11 12.58 8.20 53% 94 O restante das previsões teóricas, as quais não possuem dados experimentais, seguem na Tabela F1 do material suplementar, após este artigo. Deshpande et al.[19] também realizaram previsões da polarizabilidade. No modelo utilizado por eles, a 0K, o erro médio com o experimento foi de 21%. Desse modo, os autores propuseram uma abordagem de aquecimento teórico. A uma temperatura finita de 600K, usando esse modelo, o erro caiu para 13%. Neste trabalho, obtivemos um erro médio de 15% a 0K. De qualquer modo, os resultados apresentados por Deshpande et al.[19] corroboram qualitativamente o que observamos no presente trabalho. A troca de átomos de sódio por átomos de lítio apresenta uma queda quase linear da polarizabilidade, como se espera intuitivamente, já que átomos de lítio são menos polarizáveis do que átomos de sódio. No entanto, esse comportamento não é observado experimentalmente (Figura 6). Figura 6: Polarizabilidade em função da composição do cluster Na8-xLix. Há duas anomalias na Figura 6. A primeira é um aumento de 15% na polarizabilidade: quando dois átomos de lítio são substituídos por átomos de sódio, a expectativa era a de que trocar átomos de sódio por lítio diminuísse a polarizabilidade, e não aumentasse. A segunda anomalia é a troca de um átomo de sódio (NaLi8) no cluster de lítio puro (Li8). Adicionar átomos de sódio deveria aumentar a polarizabilidade, porém, nesse caso, houve uma diminuição brusca de 21%. Usando o mesmo tipo de base utilizado neste trabalho, Zhao et. al.[34] obtiveram polarizabilidades dentro do erro experimental para clusters de Au. Nossos resultados de lítio puro e sódio puro estão em conformidade com Chandrakumar et. al., [20] o que indica haver uma anomalia nesse sistema, seja nos modelos teóricos usuais na previsão da polarizabilidade, seja na técnica experimental. 95 Estruturas O cluster Li9 possui simetria C4v e é a base da célula unitária do lítio metálico. [33] Com exceção do NaLi8, os clusters ricos em lítio tendem a permanecer em uma simetria C4v distorcida. A partir Na6Li3, os clusters passam a ter um formato de agulha com o lítio no centro, convertendo-se bruscamente para uma estrutura C2v no sódio puro. Observa-se que apenas um átomo de lítio causa um efeito pronunciado nas estruturas de sódio. Esse efeito também é verificado no cluster com 10 átomos. Na Figura 7 todas as estruturas estão apresentadas. Figura 7: Clusters de NaxLiy para x+y=9. Nas estruturas de 10 átomos que constituem a Figura 8, o Na5Li5 apresentou uma geometria tipo D2d. Todas as demais se mostraram como variações do cluster de lítio puro. Novamente, uma mudança drástica na geometria é observada entre o Na9Li e o Na10. Figura 8: Clusters de NaxLiy para x+y=10. 96 Análise Topológica dos sistemas NaxLiy (x+y=6) A Figura 9 mostra o gráfico molecular dos caminhos de ligação para os sistema NaxLiy (x+y=6). Em verde, temos os pontos críticos de ligação (BCP); em vermelho, os pontos críticos de anel (RCP); em azul os pontos críticos de gaiola (CCP); em rosa, os pseudoátomos (atratores não nucleares, NNA). De acordo com a QTAIM, um caminho de ligação une dois átomos ligados entre si.[35] Contrariamente, nos exemplos da Figura 9, verificamos que os caminhos de ligação unem um átomo e um pseudoátomo (atrator não- nuclear), o qual existe apenas na definição topológica quando todos os três autovalores do Laplaciano são negativos em um dado ponto, mas não havendo nenhum núcleo associado a ele.[36-38] Esses caminhos de ligação envolvendo pseudoátomos existem apenas em alguns tipos de clusters metálicos. Compostos de ligação coordenada, covalente, iônica e de interações de van der Waals não possuem pseudoátomos. Diferentemente dos outros clusters NaxLiy(x+y=6), no sistema Na5Li1, há caminhos de ligação unindo os átomos de Li e Na, sem a existência de pseudoátomos (Figura 9), sugerindo que os pseudoátomos se devem à existência de pelo menos dois átomos de lítio nos clusters NaLi. O DI indica a quantidade de elétrons em um par atômico (e.g. numa ligação simples C-C, o DI é igual a um), havendo uma relação direta entre ordem de ligação formal e o DI.[39] Apesar de todo caminho de ligação representar uma interação ligada, [35] nem toda interação ligada é representada pelo caminho de ligação, podendo ser, porém, indicada pelo DI.[40] No sistema Na1Li5, os pares atômicos Na1-Li2 (DI=0,352) e Na1-Li6 têm 35% do DI de uma ligação C-C simples e 40% a mais do DI de uma ligação Ti-C de um bent- metaloceno.[40] O DI desses pares é próximo ao DI do sistema NaLi (DI=0,321). Todos os outros pares Na-Li de todos os clusters NaxLiy (x+y=6) possuem DI com ordem de grandeza 10-2, portanto menores que o DI do NaLi. Em relação às distâncias interatômicas, a distância entre Na e Li no sistema NaLi é de 2,975 Angstroms - que é bem menor do que a interação mais curta Na-Li de qualquer cluster NaxLiy (x+y=6), e.g., a distância Na-Li de 3,279Å no cluster Na3Li3 (ver Figura S1 do material suplementar). Isso mostra que, aumentando o tamanho do cluster, diminui as interações Na-Li (observado pelo DI) porque aumentam as distâncias interatômicas. Contudo, há duas exceções em que as distâncias dos pares atômicos Na1-Li2 e Na1-Li6 (3,149 Å) no Na1Li5 são maiores do que aquela do Na-Li no sistema NaLi, mas todas esses pares têm DI's equivalentes ao do sistema NaLi. 97 Em se tratando das interações Li-Li, todos os pares de todos os sistemas estudados possuem DI's próximos aos valores de DI do sistema Li4 de menor energia (DI entre 0,014 e 0,032 – ver Figura S2 do material suplementar). O sistema Li2 apresenta DI de 0,146, muito superior aos valores de DI dos pares Li-Li dos clusters. Os valores dos pares Li-Li dos clusters são equivalentes àqueles do complexo Ar-Ar, sendo então interações fracas. Para o caso das interações Na-Na, observamos a conhecida relação inversa entre a distância interatômica e o DI. Por exemplo, o par atômico Na2-Na3 do Na3Li3 tem DI 0,0651 e distância interatômica de 4,351 Å; enquanto o par Na1-Na3 do Na4Li2 tem 0,236 e 3,868Å para o DI e distância interatômica, respectivamente; e o par Na4-Na5 do Na5Li1 tem 0,517 e 3,825 Å para o DI e distância interatômica, respectivamente. Além disso, o DI dos pares Na- Na do Na4Li2 se aproxima dos valores de DI do sistema Na4 de menor energia (DI médio de 0,250), ao passo que o valor de DI dos pares Na-Na do Na5Li1 se aproxima do DI do sistema Na2 (DI 0,683). Em suma, comparado ao DI do Ar-Ar, as interações Li-Li e Na-Li (salvo dois pares Na-Li do Na1Li5) são bem mais fracas do que as interações Na-Na nos clusters Na4Li2 e Na5Li1. Figura 9: Gráfico molecular dos caminhos de ligação para os sistemas a) Na1Li5, b) Na2Li4, c) Na3Li3, d) Na4Li2 & e) Na5Li1. Em verde, temos os pontos críticos de ligação (BCP); em vermelho, os pontos críticos de anel (RCP); em azul, os pontos críticos de gaiola (CCP); em rosa, os atratores não nucleares do ponto crítico (NNACP). 98 A Figura 10 apresenta gráfico de isodensidade do Laplaciano para o sistema NaxLiy (x+y=6). A concentração de carga na camada de valência (VSCC) ocorre em interações compartilhadas (ligação covalente). Nesse tipo de interação, a concentração da densidade de carga ocorre na vizinhança da superfície interatômica, estendendo sobre as bacias atômicas dos átomos ligados e formando a sobreposição do VSCC. Através da análise do gráfico da superfície de isodensidade do Laplaciano (Figura 10), observamos nos sistemas Na1Li5, Na2Li4, Na3Li3 e Na4Li2 que a concentração de carga na camada de valência (VSCC) está localizada entre os átomos de lítio para os sistemas. Esses sistemas possuem 3, 2, 1, 3 pontos críticos (3,-3) da topologia do Laplaciano da densidade de carga, respectivamente (Ver Figura S3 do material suplementar). Todos esses pontos críticos estão próximos dos átomos de lítio, exceto no sistema Na4Li2, em que há dois pontos críticos (3,-3), próximos ao Na1, Na2 e Li5. Essa situação nos leva a inferir que a interação Li-Li apresenta topologia ligeiramente mais próxima à de uma ligação covalente em relação às interações Na-Li e Na-Na, apesar da interação Li-Li ser mais fraca do que as demais. Essa afirmação é reforçada pelo fato de que no sistema Na5Li1 não há VSCC (e nem NNA). Nesse sistema, a densidade de carga está distribuída de forma aparentemente homogênea ao longo dos átomos do sistema, como observado por meio do gráfico da Figura 10. Essa homogeneidade do Na5Li1 também é verificada pelos valores semelhantes de ρb na Tabela 4. A Tabela 4 apresenta os valores de (ρb), (∇2ρb), ε, |λ1|/λ3 e G(rc)/ρ(rc) dos pontos críticos de ligação que envolvam os pares átomo-pseudoátomo e os pares atômicos, excluindo aqueles que sejam similares por questão de simetria para os sistemas Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1. O sinal do Laplaciano da densidade de carga no ponto crítico da ligação (∇2ρb) é o responsável por determinar se o ponto critico possui concentração ou dispersão de densidade de carga. Uma interação compartilhada entre dois núcleos ocorre quando ∇2ρb<0 e ρb apresentam altos valores (10-1 ua.).[41] Além disso, em uma interação compartilhada |λ1|/λ3>1 e G(rc)/ρ(rc) <1. Uma interação de camada fechada ocorre quando ∇2ρb>0 e ρb apresentam baixos valores (10-2 ua).42-44 Em uma interação de camada fechada, |λ1|/λ3<1 e G(rc)/ρ(rc) >1. Ressaltamos que as interações de van der Waals, ligações de hidrogênio e ligações iônicas são classificadas como de camada fechada, enquanto ligações covalentes são classificadas como interações compartilhadas. Verificamos ainda na Tabela 4 que, para os sistemas Na1Li5, Na2Li4, Na3Li3 e Na4Li2, os maiores valores de ρb envolvem átomos de lítio na região do VSCC (Figura 10), reforçando 99 a ideia de similaridade topológica das interações Li-Li com a ligação covalente em termos de concentração de densidade de carga na região de ligação. Figura 10: Gráfico de isodensidade do Laplaciano para os sistemas a)Na1Li5, b) Na2Li4, c) Na3Li3, d) Na4Li2 & e) Na5Li1. Analisando os dados da Tabela 4, observamos que as interações dos sistemas Na1Li5, Na2Li4, Na3Li3 e Na5Li1 apresentaram dados topológicos do ponto crítico equivalentes aos de interações intermediárias. Um dado que chama atenção nessa tabela é o de que o sistema Na4Li2 apresentou Laplaciano negativo nos BCP’s (6 e 7). Isso significa que temos uma concentração de carga entre nos átomos Na1 e o NNA II e, por simetria, entre o Na2 e o NNA II. Além disso, apresentaram |λ1| / λ3 > 1, demonstrando que a densidade eletrônica está mais concentrada nas vizinhanças da superfície interatômica. Essa situação ocorre em ligações do tipo compartilhada,[25] sendo que essas “interações Na-NNA” apresentam maior similaridade em termos topológicos às ligações com caráter covalente. As demais interações do sistema Na4Li2 não podem ser classificadas nem como interações de camada fechada, nem compartilhada, pois os dados topológicos dos pontos críticos de ligação não se encaixam em nenhuma das duas classificações, tratando-se de interações intermediárias.[45] 100 Tabela 4: Valores de (ρb), (∇2ρb), ε, |λ1|/λ3 e G(rc)/ρ(rc) dos BCP’s da figura 9 para os sistemas Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1. Na1Li5 BCP Átomos ρb / ua. ∇2ρb / ua. Ε |λ1| / λ3 G(rc)/r(rc) 1 Na1 – NNA I 0,0098 0,0021 0,1338 0,3173 0,1697 2 Li2 – NNAI 0,0088 0,0099 8,9106 0,1659 0,3759 3 Li2 – NNAII 0,0103 0,0157 2,8021 0,1868 0,4520 4 Li4 – NNAIII 0,0118 0,0111 36,2229 0,2443 0,3515 5 Li3 – NNAII 0,0125 0,0103 0,3592 0,2292 0,3284 6 Li5 – NNAIII 0,0125 0,0103 0,3620 0,2292 0,3286 Na2Li4 BCP ρb / ua. ∇2ρb / ua. Ε |λ1| / λ3 G(rc)/r(rc) 1 Na1-NNAI 0,0091 0,0015 0,0824 0,3325 0,1514 2 Li3-NNAI 0,0110 0,0089 0,6279 0,2356 0,3136 3 Li6-NNAI 0,0092 0,0131 2,8156 0,1912 0,4229 4 Li4-NNAI 0,0111 0,0089 0,6255 0,2357 0,3134 Na3Li3 BCP Átomos ρb / ua. ∇2ρb / ua. Ε |λ1| / λ3 G(rc)/r(rc) 1 Na3 - NNAI 0,0088 0,0004 0,0388 0,4382 0,1301 2 Li5 - NNAI 0,0077 0,0077 2,0029 0,1817 0,3418 3 Li5 - NNAII 0,0086 0,0111 3,0393 0,1950 0,3987 4 Na1 - NNAII 0,0089 0,0028 0,3851 0,2804 0,1754 5 Li6 - NNAII 0,0105 0,0068 1,1468 0,2498 0,2876 6 Na2 - NNAII 0,0089 0,0029 0,4152 0,2789 0,1766 Na4Li2 BCP Átomos ρb / ua. ∇2ρb / ua. Ε |λ1| / λ3 G(rc)/r(rc) 1 Li6 – NNAI 0,0092 0,0028 5,7241 0,4858 0,2258 2 Li5 – NNAI 0,0082 0,0067 5,1924 0,3553 0,3084 3 Li5 – NNAII 0,0070 0,0035 16,7431 0,4523 0,2400 4 Na3 – NNAI 0,0080 0,0011 2,0042 0,5546 0,1424 5 Na1 – NNAI 0,0071 -0,0027 1,2384 9,9752 0,0335 6 Na1 – NNAII 0,0075 -0,0015 3,7868 2,1229 0,0842 7 Na2 – NNAII 0,0075 -0,0015 3,7868 2,1229 0,0842 Na5Li1 BCP Átomos ρb / ua. ∇2ρb / ua. Ε |λ1| / λ3 G(rc)/r(rc) 1 Na1-Li6 0,0068 0,0035 128,81 0,4493 0,2443 2 Na5-Li6 0,0068 0,0035 105,68 0,4501 0,2440 3 Na4-Li6 0,0068 0,0035 99,72 0,4499 0,2442 4 Na2-Li6 0,0068 0,0035 131,58 0,4498 0,2440 5 Na3-Li6 0,0068 0,0035 118,83 0,4489 0,2445 Ao analisarmos os dados da Tabela 5, a seguir, observamos que as energias dos átomos nos clusters em cada sistema têm valores significativamente diferentes. Podemos usar a fórmula do grau de degenerescência do D3BIA [26] para analisar o quão uniforme é a energia atômica no sistema estudado. Os sistemas Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1 têm 101 degenerescência igual a 0,500, 0,500, 0,333, 0,333 e 0,667, respectivamente. Nos sistemas em que há preponderância numérica de um elemento em relação ao outro, era de se esperar que o grau de degenerescência fosse determinado pelo número de átomos do elemento preponderante sobre o número de átomos total. Por exemplo, nos clusters Na1Li5 e Na2Li4, as degenerescências esperadas (não levando em conta nenhuma transferência de carga) seriam 0,833 e 0,667, respectivamente. Contudo, no Na1Li5, há três átomos de lítio com energias próximas (Li3 à Li5), de acordo com o critério de degenerescência do D3BIA para sistemas em gaiola,[26] e outros dois de menor energia (Li2 e Li6). Pelo fato de Li2 e Li6 estarem mais próximos do Na1 (Ver Figura S1 do Material Suplementar) que os demais átomos de Li, podemos supor que houve transferência de carga evolvendo Li2, Li6 e Na1. Verificamos na Tabela 5 que Li2 e Li6 são menos positivos que Li3 à Li5, indicando transferência de carga de Li2 e Li6 para Na1. No caso do Na2Li4, há dois grupos de átomos de lítio (cada um com dois Li) de energias próximas, mostrando novamente a influência do sódio na alteração da energia dos átomos de lítio em relação a um sistema de Li puro (e.g., no Li-Li, a energia do lítio é - 7,3878 ua.). A mesma ideia é válida para os sistemas Na4Li2 e Na5Li1. No caso do Na3Li3, em que se esperava degenerescência igual a 0,500, há um valor inferior (0,333) devido ao Na3 e Li6 terem energias muito diferentes dos demais átomos de sódio e lítio, respectivamente, provavelmente porque tanto o Na3 quanto o Li6 estão muito afastados dos demais átomos (Ver Figura S2 do Material Suplementar). Neste caso, em relação ao lítio, o Li4 e o Li5 (que estão mais próximos do Na3 – com distância de 3,282 A) têm as maiores cargas positivas de todo o sistema (0,8486 e 0,8485 ua., respectivamente). Em relação ao sódio, a carga do Na3 é a menos positiva de todos os átomos de Na do sistema (0,2620 ua.). Essas informações sugerem que os átomos Li4 e Li5 estão doando carga para Na3. 102 Tabela 5: Carga atômica e contribuição energética por átomo para os sistemas NaxLiy (x + y=6) Na5Li1 Átomos Na1 Li2 Li3 Li4 Li5 Li6 q (Ω) / e 0,3089 0,5665 0,7097 0,7746 0,7102 0,5634 E(Ω) / H -161,8719 -7,3661 -7,3960 -7,3855 -7,3959 -7,3667 Na2Li4 Átomos Na1 Na2 Li3 Li4 Li5 Li6 q (Ω) ) / e 0,2980 0,2981 0,7320 0,7315 0,8475 0,8476 E(Ω) / H -161,8655 -161,8655 -7,3813 -7,3815 -7,3416 -7,3415 Na3Li3 Átomos Na1 Na2 Na3 Li4 Li5 Li6 q (Ω) ) / e 0,3759 0,3806 0,2620 0,8486 0,8485 0,7048 E(Ω) / H -161,8562 -161,8558 -161,8665 -7,3396 -7,3397 -7,3834 Na4Li2 Átomos Na1 Na2 Na3 Na4 Li5 Li6 q (Ω) ) / e 0,0858 0,0917 0,3057 0,3028 0,8397 0,6491 E(Ω) / H -161,8683 -161,8679 -161,8555 -161,8557 -7,3493 -7,3866 Na5Li1 Átomos Na1 Na2 Na3 Na4 Na5 Li6 q (Ω) / e -0,1593 -0,1940 -0,1563 -0,2787 -0,2553 0,7087 E(Ω) / H -161,8860 -161,8882 -161,8859 -161,8935 -161,8919 -7,3549 Realizamos cálculos dos índices de aromaticidade NICS e D3BIA para os sistemas Na1Li5, Na2Li4 e Na3Li3 porque nos demais não apresentam em suas estruturas de menor energia geometria que possibilite cálculo de NICS e D3BIA. Os valores de NICS foram pouco negativos (-0,60; -0,31 e -1,83, respectivamente), assim como os valores de D3BIA foram baixos (-0,05; 0,21e 0,14 respectivamente). Como valores baixos de NICS e D3BIA estão associados a sistemas não aromáticos, por exemplo, no caso dos tetraedranos de silício e germânio[46], os sistemas Na1Li5, Na2Li4 e Na3Li3 podem ser considerados como não aromáticos ou pouco aromáticos. Conclusões Usamos um algoritmo genético com interações entre átomos descritas pelo método MP2 para estudar as geometrias de clusters NaxLiy com 4 ≤ x+y ≤ 10. Cálculos do potencial 103 de ionização e da polarizabilidade foram realizados, os quais se mostraram de acordo com o experimento de 6% e 13%, respectivamente. Na polarizabilidade, foi observada uma anomalia que não pode ser compreendida com os métodos teóricos deste trabalho nem os empregados pela referência.[19] Estruturas de clusters com 9 e 10 átomos foram reportadas pela primeira vez. Os clusters com 9 átomos ricos em lítio se assemelham a célula unitária deste metal e todos os clusters bimetálicos com 10 átomos encontrados são variações do Li10. Também observamos que a escolha do método possui um papel fundamental na descrição da geometria dos clusters, assim, uma mudança no método ou no funcional pode levar a geometrias finais completamente diferentes. A QTAIM pode ser usada em sistemas de clusters metálicos para analisar: (1) as interações entre pares atômicos através do índice de deslocalização e dados topológicos do ponto crítico de ligação; (2) a concentração de carga em dada região do cluster a partir da topologia do Laplaciano; e (3) transferência de carga entre átomos analisando-se as energias e cargas atômicas. Nos clusters NaxLiy (x+y=6), chama atenção o fato de não haver caminho de ligação envolvendo diretamente os metais, sendo unidos por pseudoátomos, com exceção do Na5Li1. Algumas interações atômicas não foram indicadas pelo caminho de ligação e sua análise foi feita pelo DI. No sistema Na1Li5, os pares atômicos Na1-Li2 e Na1-Li6 têm as interações mais fortes (e equivalentes à do sistema NaLi) de todos os pares Na-Li de todos clusters NaxLiy(x+y=6), enquanto os outros pares Na-Li têm interações dez vezes mais fracas que aquelas do sistema NaLi. As interações Na-Na dos clusters Na4Li2 e Na5Li1 são as mais fortes, quando comparadas com sistemas puros. Portanto, a fórmula do grau de degenerescência do D3BIA e a carga atômica indicaram que os átomos de lítio mais próximos ao átomo de sódio transferem carga para esse último. Referências [1] Ballardini, R.; Balzani, V.; Credi, A.; Gandolfi, M. T.; Venturi, M. Acc. Chem. Res. 2001, 34 (6), 445–455. [2] Fennimore, A. M.; Yuzvinsky, T. D.; Han, W.-Q.; Fuhrer, M. S.; Cumings, J.; Zettl, A. Nature 2003, 424 (6947), 408–410. [3] Hogg, A. C. and B. S. and R. A. F. J. and T. Nanotechnology 2008, 19 (1), 15103. [4] Gireaud, L.; Grugeon, S.; Laruelle, S.; Yrieix, B.; Tarascon, J.-M. Electrochem. commun. 104 2006, 8 (10), 1639–1649. [5] Yagi, J.; Tanaka, T.; Muroga, T.; Sagara, A. Fusion Eng. Des. 2014, 89 (7–8), 1168–1171. [6] Silva, M. X.; Galvão, B. R. L.; Belchior, J. C. J. Mol. Model. 2014, 20 (9), 2421. [7] Mohajeri, A.; Mahmoodinia, M. J. Iran. Chem. Soc. 2013, 10 (6), 1229–1237. [8] Benichou, E.; Allouche, A. R.; Aubert-Frecon, M.; Antoine, R.; Broyer, M.; Dugourd, P.; Rayane, D. Chem. Phys. Lett. 1998, 290 (1–3), 171–179. [9] Pérez, J. F.; Florez, E.; Hadad, C. Z.; Fuentealba, P.; Restrepo, A. J. Phys. Chem. A 2008, 112 (25), 5749–5755. [10] Dzuba, V. A.; Safronova, M. S.; Safronova, U. I.; Kramida, A. Phys. Rev. A 2016, 94 (4), 42503. [11] Mosnier, J.-P.; Kennedy, E. T.; van Kampen, P.; Cubaynes, D.; Guilbaud, S.; Sisourat, N.; Puglisi, A.; Carniato, S.; Bizau, J.-M. Phys. Rev. A 2016, 93 (6), 61401. [12] Kang, Y.; Jeon, S. H.; Cho, Y.; Han, S. Phys. Rev. B 2016, 93 (3), 35131. [13] Tao, J.; Mo, Y.; Tian, G.; Ruzsinszky, A. Phys. Rev. B 2016, 94 (8), 85126. [14] Loboda, O.; Ingrosso, F.; Ruiz-López, M. F.; Reis, H.; Millot, C. J. Comput. Chem. 2016, 37 (23), 2125–2132. [15] Kalugina, Y. N.; Thakkar, A. J. Chem. Phys. Lett. 2016, 644, 20–24. [16] Valiev, R. R.; Minaev, B. F. J. Mol. Model. 2016, 22 (9), 214. [17] Centeno, J.; Fuentealba, P. Int. J. Quantum Chem. 2011, 111 (7-8), 1419–1435. [18] Tai, T. B.; Nhat, P. V.; Nguyen, M. T. Phys. Chem. Chem. Phys. 2010, 12 (37), 11477– 11486. [19] Deshpande, M. D.; Kanhere, D. G.; Panat, P. V; Vasiliev, I.; Martin, R. M. Phys. Rev. A 2002, 65 (5), 53204. [20] Chandrakumar, K. R. S.; Ghanty, T. K.; Ghosh, S. K. J. Phys. Chem. A 2004, 108 (32), 6661–6666. [21] Deshpande, M. D.; Kanhere, D. G.; Vasiliev, I.; Martin, R. M. Phys. Rev. A 2002, 65 (3), 33202. [22] Silva, F. T.; Galvão, B. R. L.; Voga, G. P.; Silva, M. X.; Rodrigues, D. D. C.; Belchior, J. C. Chem. Phys. Lett. 2015, 639, 135–141. [23] Alexandrova, A. N.; Boldyrev, A. I. J. Chem. Theory Comput. 2005, 1 (4), 566–580. [24] Rondina, G. G.; Da Silva, J. L. F. J. Chem. Inf. Model. 2013, 53 (9), 2282–2298. [25] Bader, R. F. W. Atoms in Molecules - A Quantum Theory; Oxford University Press:Oxford, 1990. [26] Araújo, D. M.; Da Costa, T. M.; Firme, C. L. J Mol Model 2015, 21, 248. 105 [27] Chen, Z. F.; Wannere, C. S.; Corminboeuf, C.; Puchta, R.; Schleyer, P. V. Chem. Rev. 2005, 105, 3842-3888. [28] Schimidt, M. W.; Baldridge, K. K.; Boats, J. A.; Elbert, S. T.; Gorgon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery Jr., J. A.; J. Comput. Chem., 14 (1993) 1347. [29] Frisch, M. J.;Trucks, G. W.; Schlegel, H.B.; Scuseria, G. E.;, Robb, M. A.; Cheeseman, J. R.; et al. Gaussian 09, Revision A.02, Gaussian Inc Wallingford CT, 2009, 34, Wallingford CT. [30] Gready, J. E.; Bacskay, G. B.; Hush, N. S. Chem. Phys. 1977, 22 (1), 141–150. [31] AIMAll (Version 12.06.03), Keith, T. A.; TK Gristmill Software, Overland Park KS, USA, 2016 (aim.tkgristmill.com) [32] Biegler-Konig, F.; Schonbohm, J.; Bayles, D.; J. Comput. Chem., 22 (2001) 545−559. [33] Goel, N.; Gautam, S.; Dharamvir, K. Int. J. Quantum Chem. 2012, 112 (2), 575–586. [34] Zhao, Y.-R.; Kuang, X.-Y.; Zheng, B.-B.; Li, Y.-F.; Wang, S.-J. J. Phys. Chem. A 2011, 115 (5), 569–576. [35] Bader, R. F. W.; J. Phys. Chem. A, 2009, 113 (38), 10391–10396. [36] Gatti, C.; Fantucci, P.; Pacchioni, G.; Theor. Chim. Acta, 1987, (72) 433–458. [37] Caio, W. L.; Gatti, C.; MacDougall, P. J.; Bader, R. F. W.; Chem. Phys. Lett., 1987, (141), 380. [38] Pendás, A. M.; Blanco, M. A.; Costales, A.; Sánchez, P. M.; Luaña, V.; Phys. Rev. Lett., 1999, (83) 1930-1933. [39] Firme, C. L.; Antunes, O. A. C.; Esteves, P. M.; Chem. Phys. Let. 2009, 468, 129-133. [40] Dos Santos, H. F. L.; Pontes, D. L.; Firme, C. L.; J. Mol. Model. 2013, (19), 2955-2964 [41] Bader, R. F. W. e Essen, H.; J. Chem. Phys., 1984, (80), 1943- 1984. [42] Bader, R. F. W.; Phys. Rev. B, 1994 , (49), 13348. [43] Popelier, P. L. A.; Atoms in Molecules An Introduction, Manchester: Prentice Hall (2000). [44] Gibbs, G. V.; Spackman, M. A.; Jayatilaka, D.; Rosso, K. M. e Cox, D. F.; J. Phys. Chem. A, 2006, (110), 12259. [45] Firme, C. L.; Ponte, D. L.; Antunes, O. A. C.; Chem. Phys. Let. 2010, 499, 193–198. [46] Monteiro, N. K. V.; De Oliveira, J. F.; Firme, C. L.; NewJ.Chem., 2014, 38, 5892. 106 Material Suplementar Tabela F1: Potencial de ionização e polarizabilidade das demais estruturas de NaxLiy sem dados experimentais. Sistema IP α/n Na1Li4 3,93 15,03 Na1Li5 4,19 13,74 Na1Li6 3,75 11,31 Na1Li8 3,59 10,58 Na1Li9 3,83 12,61 Na2Li4 4,10 14,90 Na2Li5 3,68 12,17 Na2Li7 3,66 11,27 Na2Li8 3,89 13,06 Na3Li4 3,61 12,61 Na3Li6 3,68 11,78 Na3Li7 3,72 13,82 Na4Li5 3,51 12,02 Na4Li6 3,84 14,17 Na5Li4 3,52 12,56 Na5Li5 3,68 15,76 Na6Li4 3,56 15,63 Tabela F2: Valores do RDF, DIU, Degenerescência e D3BIA para os sistemas Na1Li5, Na2Li4 e Na3Li3. Molécula RDF DIU Degenerescência D3BIA Na1Li5 0,0084 -12,187 0,500 -0,0515 Na2Li4 0,0074 55,565 0,500 0,2055 Na3Li3 0,0072 57,971 0,333 0,1393 107 Figura S1: Distâncias interatômicas em Angstrons para os sistemas a) Na1Li5, b) Na2Li4, c) Na3Li3, d) Na4Li2 & e) Na5Li1. Figura S2: Distâncias interatômicas e valores de DI’s para os sistemas a) Li4 & b) Na4. 108 Figura S3: Representação do VSCC (3, -3) para os sistemas a) Na1Li5, b) Na2Li4, c) Na3Li3, d) Na4Li2 & e) Na5Li1. 109 Capítulo 7 – Conclusões gerais 110 Ao se compararem as estruturas obtidas por meio do GA com potencial Gupta para clusters de alumínio puro com resultados recentes da literatura, verificou-se que para os sistemas Al2, Al3, Al6, Al8, Al9, Li5, Li6, Li7, Al1Li5, Al1Li7 e Al1Li8 as geometrias obtidas foram muito semelhantes àquelas obtidas por meio de cálculos de funcional de densidade e ab initio[como CCSD(T)].[23-25] As estruturas AlxLiy, AlxNay e AlxKy apresentaram átomos de alumínio do interior da liga, sendo que os átomos de lítio, sódio e potássio estiveram sempre nas periferias dos clusters. O processo de crescimento das nanoligas com 6 e/ou 7 átomos parecem seguir como base a estrutura octaédrica e bipirâmide pentagonal, respectivamente. Para os clusters de alumínio-sódio e alumínio-lítio foi observado que a proporção de 30% a 40% de alumínio foi a que apresentou menores energias de excesso ao longo do crescimento dos sistemas AlxLiy e AlxNay. Para o sistema de alumínio-potássio, observou-se que essa proporção é de 40% a 50%. De um modo geral, o GA com implementação do potencial Gupta mostrou bons resultados, sendo as estruturas encontradas bem próximas daquelas de mínimo global; porém, tal método apresentou falhas na obtenção das estruturas de Al4, Al5 e Al7, demonstrando que, para sistemas muito pequenos, o potencial Gupta, não descreve tão bem as interações atômicas. Em se tratando do cálculo de pequenas estruturas NaxLiy (x+y≤10), constatou-se que o Algoritmo Genético quântico (Q-GA) teve ótimos resultados, sendo, portanto, mais confiável para o tratamento de estruturas pequenas. Entretanto, por ser um método quântico, o Q-GA apresenta alto custo computacional, não possuindo a mesma abrangência do GA clássico para abordagem de sistemas maiores. Para o sistema de sódio-lítio, foram realizados cálculos topológicos a partir da Teoria Quântica de Átomos em Moléculas (QTAIM) para as estruturas Na1Li5, Na2Li4, Na3Li3, Na4Li2 e Na5Li1 obtidas pelo Q-GA. Observou-se que não há ligação direta entre os átomos, e sim de forma, indireta, através de pseudoátomos que são definidos apenas na QTAIM. Nos clusters NaxLiy (x+y=6) foi observado que não há caminho de ligação envolvendo diretamente os metais, sendo unidos por pseudoátomos, com exceção do Na5Li1. Algumas interações atômicas não foram indicadas pelo caminho de ligação e sua análise foi feita pelo DI. No sistema Na1Li5, os pares atômicos Na1-Li2 e Na1-Li6 têm as interações mais fortes (e equivalentes à do sistema NaLi) de todos os pares Na-Li dos sistemas estudados, enquanto os outros pares Na-Li possuem interações dez vezes mais fracas do que aquelas do sistema NaLi. As interações Na-Na dos clusters Na4Li2 e Na5Li1 são as mais fortes, quando comparadas com sistemas puros. 111 Por fim, a fórmula do grau de degenerescência do D3BIA e a carga atômica indicaram que os átomos de lítio mais próximos ao átomo de sódio transferem carga para esse último. Em trabalhos futuros, pretendemos realizar os cálculos MP2 e CCSD(T) para as estruturas de menores energias de excesso para os sistemas alumínio-lítio, alumínio-sódio e alumínio-potássio, no intuito de analisarmos melhor a eficiência do GA clássico. Também faz-se necessário o estudo de outros sistemas através do Q-GA, de modo a termos maior representatividade sobre a eficiência desse método. Uma outra proposta é implementar o Q- GA acoplado ao pacote ORCA,[102] uma vez que o mesmo trata-se de um programa de estrutura eletrônica mais moderno, sendo mais abrangente do que o GAMESS-US.[103] 112 Referências gerais 113 [1] Martin, T. P.; Näher,U.; Schaber, H.; Chem. Phys. Lett. 199, 470 (1992). [2] Ferrando, R.; Jellinek, J.; Johnston, R. L.; Chem. Rev. (Washington, D.C.) 108 (2008) 845. [3] Bordón, L. O. P.; Computational Studies of Transition Metal Nanoalloys, 2011, 173 f. Tese (Doutorado em Química), School of Chemistry, University of Birmingham. Birmingham. 2011. [4] Cox H.; Johnston, R. L.;Murrell J. N. J.; Solid State Chem. 145 517, (1999) [5] Rao, B. K.; Jena, P., J. Chem. Phys. (1999) 111 1890. [6] Berry, R. S.; Braier, P.; Hinde, R. J.; Cheng, H. P.; Israel J. Chem., 30 39, 1990. [7] Hartke B.., Chem. Phys. Lett., 258 144, 1996. [8] Wales D. J. e Scheraga H. A., Science, 285 1368, 1999. [9] White R. P. e Mayne H. R., Chem. Phys. Lett., 289 463, 1998. [10] Holland J. H., Sci. Am., 267 66, 1992. [11] Hartke B., Phys. Chem, 214 1251, 2000. [12] White, R. P.; Niesse, J. A.; Mayne, H. R.; J. Chem. Phy., 108 2208, 1998. [13] Doye, J. P. K.; Wales, D. J.; Berry, R. S.; J. Chem. Phys., 103 4234, 1995. [14] Doyce, J. P. K.; Wales D. J., J. Chem. Soc. Faraday, Trans., 93 4233, 1997. [15] Rodrigues, D.D.C.; Nascimento, A.M.; Duarte, H.A.; Belchior, J.C., Chemical Physics 349 (2008) 91. [16] Baletto, F.; Ferrando, R., Rev. Mod. Phys. 77 (2005) 371. [17] Lordeiro, R. A.; Guimarães, F. F.; Belchior, J. C.; Johnston, R. L.; International Journal of Quantum Chemistry 95 (2003) 112. [18] Cleri, F.; Rosato, V.; Phisical Review B, 48 nº1 (1993) 22. [19] Solov’yov, I. A.; Solov’yov, A. V.; Greiner, W.; Phys. Rev. A, 65, 053203 (2002). [20] Banerjee, A.; Ghanty, T. K; Chakrabarti, A.; J. Phys. Chem. A, 112, 12303-12311 (2004). [21] Hauser, A. W.; Callegari, C.; Soldán, P.; Ernest, W. E.; J. Chem. Phys., 129, 044307 (2008). [22] H. P. Cheng, R. N. Barnett and U. Landman, Phys. Rev. B: Condens. Matter, 1993,48, 1820. [23] Alexandrova, A. N.; Boldiver, A. I.; J. Chem. Theory Comput. 1, 566 (2005). [24] Kiohara, O. V.; Carvalho, E. F. V; Paschoal, C. W. A.; Machado, F. B. C., Chem. Phys. Leters 568-569 42 (2013). [25] Tai, T. B.; Nhat, P. V.; Nguyen, M. T.; Phys. Chem. Chemical Physics, 12, 11477 (2010). 114 [26] Gireaud, L.; Grugeon, S.; Laruelle, S.; Yrieix, B.; Tarascon, J.-M. Electrochem.commun. 2006, 8 (10), 1639–1649. [27] Yagi, J.; Tanaka, T.; Muroga, T.; Sagara, A. Fusion Eng. Des. 2014, 89 (7–8), 1168– 1171. [28] Silva, M. X.; Galvão, B. R. L.; Belchior, J. C. J. Mol. Model. 2014, 20 (9), 2421. [29] Mohajeri, A.; Mahmoodinia, M. J. Iran. Chem. Soc. 2013, 10 (6), 1229–1237. [30] Benichou, E.; Allouche, A. R.; Aubert-Frecon, M.; Antoine, R.; Broyer, M.; Dugourd, P.; Rayane, D. Chem. Phys. Lett. 1998, 290 (1–3), 171–179. [31] Pérez, J. F.; Florez, E.; Hadad, C. Z.; Fuentealba, P.; Restrepo, A. J. Phys. Chem. A 2008, 112 (25), 5749–5755. [32] Guimarães, F. F.; Belchior, J. C.; Johnston, R. L.; Roberts, C., Journal Of Chemical Physics, 116 nº19 (2002) 8327. [33] Wales, D. J.; Hodges, M. P., Chem. Phys. Lett. 286 (1998) 65. [34] Doye J. P. K., Wales D. J. e Berry R. S., J. Chem. Phys., 103 4234, 1995. [35] Roberts, C.; Johnston, R. L.; Wilson, N. T., Theor. Chem. Acc. 104 123 (2000). [36] Niesse, J. A.; Mayne, H. R., J. Comput. Chem. 18 (1997) 1233. [37] Zhu, C.; Byrd, R. H.; Lu, P.; Nocedal, J., Northwestern University, EECS Tech. Rep. NAM 12, 1995. [38] Zeiri, Y., Phys. Rev. E 51 (1995) R2769. [39] Zeiri, Y.; Fattal, E.; Kosloff, R., J. Chem. Phys. 102 (1995) 1859. [40] Disponível em: http://professor.webizu.org/ga/, acessado em 03 de Dezembro de 2016. [41] Silva, E. S. A.; Duarte, H. A.; Belchior, J. C., Chemical Physics 323 (2006) 553–562. [42] Lordeiro, R. A.; Guimarães, F. F.; Belchior, J. C.; Johnston, R. L., International Journal of Quantum Chemistry 95 (2003) 112. [43] Schlegel, H.B., Adv. Chem. Phys. 67 (1987) 249. [44] Rodrigues, D.D.C.; Nascimento, A.M.; Duarte, H.A.; Belchior, J.C., Chemical Physics 349 (2008) 91. [45] Disponível em: http://iaufes20092.pbworks.com/w/page/8625621/FrontPage. Acesso em 04/12/2016. [46] Johnston, R. L.; Roberts, Soft Computing Approaches in Chemistry, edited by H. Cartwright, L. Sztandera (Physica-Verlag, Heidelberg, 2001). [47] Lai, S. K.; Hsu, P. J.; Wu, K. L.; Liu, W. K.; Iwamatsu, M., J. Chem. Phys. 117 nº23 (2002) 10715. [48] Gupta R. P., Phys. Rev. B, 32, 8278 (1985) 115 [49] Daw, M. S.; Baskes, M. I., Phys. Rev. Lett. 50, 1285 (1983); Phys. Rev. B29, 6443 (1984). [50] Rosato, V.; Guillope, M.; Legrand, B., Philos. Mag. A 59, 321 (1989). [51] Aguado, A.; López, J. M., J. Chem. Phys. 133 (2010) 094302-1. [52] Li, Y.; Blaisten-Barojas, E.; Phisical Review B 57 nº24 (1998) 15 519. [53] Tumer, G.W.; Johnston, R. L.; Wilson, N. T; J. Chem. Phys. 112, 4773 (2000) [54] Li, Y. e Blaisten-Barojas, E.; Phys. Rev. B, 57, 24 (1998). [55] Wilson, N. T.; Johnston, R. L., Eur. Phys. J D (2000), 12, 161. [56] Rao, B. K.; Jena, P.; J. Chem. Phys. (1999) 111 1890. [57] Aguado, A. e López, J. M.; J. Chem. Phys. 133 (2010) 094302 [58] Chen, Z.; Jiang, X.; Li, J. e Li, S.; J. Chem. Phys. 138, 214303 (2013) [59] P.O.Lӧwdin, Adv. Chem. Phys. 2, 207 (1959) [60] I. Shavitt, Methods of Eletronic Structure Theory, e.d, H. F. Schaeter III, Plenum Press, New York (1997) [61] W. Duch, Lecture Notes in Chemistry, vol. 42, Springer, Berlin (1986). [62] A. Szabo e N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Eletronic Theory, Dover Publications Inc. (1996). [63] Presente em: http://www.nobelprize.org/nobel_prizes/chemistry/laureates/1998/ (aceso em 02/12/2016). [64] Schrӧdinger, E.; Ann. Physik 80, 437 (1926) [65] Morgon, N. H; Coutinho, K. Métodos de Química Teórica e Modelagem Molecular, Ed. Livraria da Fisica, São Paulo, (2007) [66] Møller, C.; Plesset, M. S.; Phys. Rev. 46, 618 (1934). [67] Krishman, R.; Pople, J. A.; Int. J. Quantum Chem.14, 91 (1978) [68] Bartlett, R. J.; Ann. Rev. Phys. Chem. 32, 359 (1981). [69] Roothaan, C. C. J.; Rev. Mod. Phys. 26 69 (1951) [70] Pople, J. A.; Nesbet, R. K.; J. Chem. Phys. 22, 571 (1954). [71] Olsen, J.; Christinasen, O.; Koch, H e Jorgensen, P.; J. Chem. Phys. 105, 5082 (1996). [72]Coster, F.; Nuclear Phys. 7, 421 (1958). [73] Kümmel, H.; Lührmann e Zabolitzki, J. G. Phys. Rep. 36, 1 (1978). [74] Čížec, J.; Chem. Phys. 45, 4256 (1966). [75] Čížec, J.; Adv. Chem. Phys. 14, 35 (1969). [76] Jorgensen, P. e Simons, J.; Second-Quantization-Based Methods in Quantum Chemistry, Acad. Press, N. York, 1981. 116 [77] Bartlett, R. J. e Silver, D. M.; Int. J. Quantum Chem. S9, 183 (1975). [78] Bartlett, R. J. e Purvis, G. D.; In. Quantum Chem. 14, 561 (1978). [79] Pople, J. A.; Krishnan, R.; Schlegel, H. B. e Binkley, J. S.; Int. J. Quantum Chem. 14, 545 (1978). [80] Paldus, J.; Čížec, J. e Shavitt, I.; Phys. Rev. A 5, 50 (1972). [81] Purvis, D. G.; Bartlett, R. J.; J. Chem. Phys. 76, 1910 (1982). [82] Noga, J.; Bartlett, R. J.; J. Chem. Phys. 86, 7041 (1987). [83] Bader, R. F. W.; Phys. Rev. B, 49 13348 (1994) [84] Hofmann, K., M. H. Prosenc, et al.; Chem. Comm., 29 3097 (2007). [85] Matta, C. F. e R. F. W. Bader.; J. Phys.l Chem. A, 110 6365 ( 2006). [86] Bader, R. F. W.; Accounts of Chem. Research, 18 9 (1985). [87] Bader, R. F. W.; Chem. Rev., 91 893(1991). [88] Srebrenik, S. e Bader, R. F. W.; Nguyen-Dang, T. T.; J. Chem. Phys. 68, 3667 (1978). [89] Wiberg, K. B.; Bader, R. F. W.; Lau, C. D. H.;. J. American Chem. Soc., 109 1001 (1987). [90] Wiberg, K. B.; Bader, R. F. W.; Lau, C. D. H.;. J. American Chem. Soc., 109 985 (1987). [91] Balanarayan, P. e Gadre, S. R.; J. Chem. Phys., 119, 5037 (2003). [92] Bader, R. F. W., P. L. A. Popelier, Keith, T. A.; Angewandte Chemie, 33 620 (1994). [93] Popelier, P. L. A.; Atoms in Molecules An Introduction, Manchester: Prentice Hall (2000). [94] Bader, R. F. W. e Essen, H.; J. Chem. Phys., 80, 1943 (1984). [95] Gibbs, G. V.; Spackman, M. A.; Jayatilaka, D.; Rosso, K. M. e Cox, D. F.; J. Phys. Chem. A, 110, 12259 (2006). [96] Bader, R. F. W.; Tang, T. H.; Tal, Y.; Biegler-Koenig, F. W.;J. American Chem. Soc., 104, 940 (1982). [97] Fradera, X.; Austen, M. A. e Bader, R. F. W; J. Phys. Chem. A, 103, 304 (1999). [98] Merino, G.; Vela, A. e Heine, T.; Chem. Reviews, 105, 3812 (2005). [99] Bader, R. F. W. e Stephens, M. E.; J. American Chem. Soc., 97 7391 (1975). [100] Bader, R. F. W.; Streitwieser, A.; Neuhaus, A.; Laidig, K. E.; Speers, P.; J. American Chem. Soc., 118 4959 (1996). [101] Malcolm, N. O. J. e Popelier, P. L. A.; J. Phys. Chem. A, 105, 7638 (2001). [102] Neese, Frank; WIREs Comput. Mol. Sci, 2, (2012) 73-78. 117 [103] Schimidt, M. W.; Baldridge, K. K.; Boats, J. A.; Elbert, S. T.; Gorgon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery Jr., J. A.; J. Comput. Chem., 14 (1993) 1347.