Estudio de modelos de ión polarizables para sales fundidas mediante dinámica molecular

ADVERTIMENT . La consulta d’aquesta tesi queda condicionada a l’acceptació de les següents condicions d'ús: La difusió d’aquesta tesi per mitjà del servei TDX (www.tesisenxarxa.net) ha estat autoritzada pels titulars dels drets de propietat intel·lectual únicament per a usos privats emmarcats en activitats d’investigació i docència. No s’autoritza la seva reproducció amb finalitats de lucre ni la seva difusió i posada a disposició des d’un lloc aliè al servei TDX. No s’autoritza la presentació del seu contingut en una finestra o marc aliè a TDX (framing). Aquesta reserva de drets afecta tant al resum de presentació de la tesi com als seus continguts. En la utilització o cita de parts de la tesi és obligat indicar el nom de la persona autora. ADVERTENCIA. La consulta de esta tesis queda condicionada a la aceptación de las siguientes condiciones de uso: La difusión de esta tesis por medio del servicio TDR (www.tesisenred.net) ha sido autorizada por los titulares de los derechos de propiedad intelectual únicamente para usos privados enmarcados en actividades de investigación y docencia. No se autoriza su reproducción con finalidades de lucro ni su difusión y puesta a disposición desde un sitio ajeno al servicio TDR. No se autoriza la presentación de su contenido en una ventana o marco ajeno a TDR (framing). Esta reserva de derechos afecta tanto al resumen de presentación de la tesis como a sus contenidos. En la utilización o cita de partes de la tesis es obligado indicar el nombre de la persona autora. WARNING. On having consulted this thesis you’re accepting the following use conditions: Spreading this thesis by the TDX (www.tesisenxarxa.net) service has been authorized by the titular of the intellectual property rights only for private uses placed in investigation and teaching activities. Reproduction with lucrative aims is not authorized neither its spreading and availability from a site foreign to the TDX service. Introducing its content in a window or frame foreign to the TDX service is not authorized (framing). This rights affect to the presentation summary of the thesis as well as to its contents. In the using or citation of parts of the thesis it’s obliged to indicate the name of the author Estudio de modelos de ión polarizable para sales fundidas mediante dinámica molecular Vicente Bitrián Varea Universitat Politècnica de Catalunya Departament de Física i Enginyeria Nuclear Universitat Politècnica de Catalunya Departament de Física i Enginyeria Nuclear Estudio de modelos de ión polarizable para sales fundidas mediante dinámica molecular Memoria de la tesis doctoral realizada por Vicente Bitrián Varea bajo la dirección del Dr. Joaquim Trullàs i Simó. Barcelona, febrero de 2011 Agradecimientos En primer lugar, quiero expresar mi más sincera gratitud a mi director de tesis, Quim Trullàs, por guiarme en este trabajo con entusiasmo, rigor y paciencia, por la dedicación constante y los consejos que me ha dado, y por la confianza que ha depositado en mí. Quiero agradecer muy especialmente el apoyo de Olga Alcaraz, que me ha brindado su amabilidad, colaboración y consejo siempre que lo he necesitado. Ha sido también un honor compartir trabajo y charlas con el profesor Moisés Silbert, que me ha obsequiado con su inestimable magisterio y su afecto entrañable. Cómo no, me acuerdo en este momento de los que han sido mis compañeros durante todo este tiempo. Ellos ya saben lo importante que es para mí su amistad y el haber compartido este camino con ellos: Michele Catanzaro, Ester Sola, Marco Masia, Andrea Baronchelli, Riccardo Rota, Paolo Moretti, Ramon Forcada, Kostas Sakkos, Claudi Cazorla, Alfredo Soldevilla, Rubén Cabezón, Nuria Serichol, Óscar Lorente, Elena López, Delfi Nieto, Agustí Emperador, Jonàs Sala, Gregory Astrakharchik, Oleg Osychenko, Yarko Lutsyshyn, Clara Prats, Daniel Laria y todos aquellos que, imperdonablemente, pueda olvidar ahora. También deseo agradecer a Ana Calle, Cristina Astier y Montserrat Río su simpatía y disposición, y a los profesores del departamento su valioso asesoramiento en mis tareas docentes, en especial a Jordi Martí, Lluís Ametller, Elvira Guàrdia, Rossend Rey y Manel Canales. Y, finalmente, no puedo dejar de acordarme de mi familia, en especial de mis padres y mi hermano, ya que sin su ayuda, aliento y paciencia este trabajo no habría sido posible. Contenidos Introducción Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido Capítulo 2. Modelos de ión polarizable Capítulo 3. Modelos de ión polarizable para AgBr fundido. Efectos de la polarización inducida Capítulo 4. Estudio de α-AgI y AgI fundido mediante simulaciones de iones polarizables Capítulo 5. Propiedades dieléctricas estáticas en modelos de ión polarizable Capítulo 6. Conclusiones y perspectivas Apéndices. A Acrónimos y símbolos B Unidades CGS y SI C Potenciales interiónicos D Desarrollo de los modelos de ión polarizable E Modelos de Thole F Derivación del teorema de fluctuación-disipación para la función de respuesta dieléctrica Índice Introducción 1 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 9 1.1 Introducción a los haluros alcalinos y haluros de plata y cobre 10 1.1.1 Haluros alcalinos ........................................................................ 10 1.1.2 Haluros de plata y cobre ............................................................. 12 1.2 Modelos de ión rígido para sales 1:1 16 1.2.1 Potenciales de Born-Mayer-Huggins .......................................... 17 1.2.2 Potenciales de Vashishta-Rahman .............................................. 19 1.3 Propiedades estructurales de sales 1:1 fundidas 22 1.3.1 Funciones de distribución radial ................................................. 22 1.3.2 Factores de estructura parciales .................................................. 24 1.3.3 Factores de estructura de Bhatia-Thornton ................................. 27 1.3.4 Factores de estructura experimentales ........................................ 31 1.4 Transporte iónico en sales 1:1 fundidas 35 1.4.1 Autodifusión ............................................................................... 35 1.4.2 Conductividad ............................................................................. 41 Referencias 46 VIII Índice D.5.5 Fuerzas en un PIMe ................................................................... 219 Apéndice E. Modelos de Thole 223 E.1 Interacciones entre cargas y dipolos no puntuales 224 E.2 Tensores multipolares corregidos y funciones amortiguadoras 227 E.3 Energía y fuerzas en un PIM' 228 Referencias 230 Apéndice F. Derivación del teorema de fluctuación-disipación para la función de respuesta dieléctrica 231 F.1 Derivadas funcionales de la función de partición grancanónica 232 F.2 Teoría de la respuesta lineal 236 Referencias 239 Introducción La tesis Estudio de modelos de ión polarizable para sales fundidas mediante dinámica molecular que presentamos ha sido elaborada en el seno del grupo de investigación SIMCON (Simulación por Ordenador en Materia Condensada), integrado en el Departamento de Física e Ingeniería Nuclear de la Universitat Politècnica de Catalunya. Este grupo utiliza la simulación computacional para investigar sistemas como sales fundidas, líquidos moleculares, soluciones iónicas, polímeros o macromoléculas. La complejidad de las simulaciones y de los sistemas objeto de estudio ha ido aumentando progresivamente durante los últimos años en paralelo con el progreso tecnológico de los ordenadores y el desarrollo de nuevos modelos y técnicas de simulación. El método de simulación usado en este trabajo es la dinámica molecular (MD, de Molecular Dynamics). A partir de un modelo microscópico, que esencialmente determina cómo interaccionan las partículas que constituyen el sistema, en las simulaciones MD se generan configuraciones sucesivas integrando simultáneamente las ecuaciones del movimiento clásicas de todas las partículas. El conjunto de configuraciones así generado pertenece a la colectividad microcanónica y, en aplicación de la mecánica estadística, promediando sobre todas ellas se obtienen las propiedades del sistema en equilibrio termodinámico. Otro método ampliamente utilizado es la simulación Monte Carlo (MC), que genera configuraciones explorando el espacio fásico de las posiciones según la distribución de probabilidad canónica. Dichas configuraciones no están ordenadas temporalmente, por lo que en las simulaciones MC, a diferencia de las MD, no se pueden determinar las propiedades dinámicas. Las primeras simulaciones MD, realizadas por Adler y Wainwright en 1959 [1], reducían las moléculas de los gases nobles a esferas duras. Poco después, Rahman [2] y Verlet [3] se basaron en potenciales continuos de 2 Introducción corto alcance de tipo Lennard-Jones para simular el argón líquido. Tras estos trabajos pioneros, el rápido desarrollo de los ordenadores en los años siguientes posibilitó la simulación de sistemas más complejos con una representación menos idealizada de sus partículas. Por ejemplo los sistemas iónicos, cuya simulación entraña una doble dificultad: en ellos existe la interacción coulombiana, que es de largo alcance y cuyo cálculo presenta dificultades técnicas, y al estar compuestos por más de una especie química hay que proponer un potencial de interacción para cada pareja de especies. El año 1971 supuso el inicio del progreso desarrollado durante las últimas décadas en el estudio de sales fundidas. Fue entonces cuando se realizó la primera simulación de una sal fundida, en concreto KCl mediante MC [4]. Durante la década que empezaba y la siguiente se publicaron numerosos trabajos basados en la simulación de haluros alcalinos fundidos cerca de su punto de fusión, tanto con MD como con MC. En un inicio se emplearon modelos de ión rígido (RIM, de Rigid Ion Model), es decir que no tienen en cuenta la deformación de la estructura electrónica de los iones, mayoritariamente con los potenciales propuestos por Fumi y Tosi [5]. Tras estos trabajos iniciales con RIM, se comenzaron a publicar simulaciones MD de haluros alcalinos que incluían la polarización inducida en los iones mediante los modelos conocidos como shell models. Un repaso a la extensa lista de trabajos de esos años se puede consultar en el review de Sangster y Dixon [6] o el libro de March y Tosi [7]. También en 1971, Page y Mika [8] hicieron por vez primera un experimento de difracción de neutrones con la técnica de la substitución isotópica para obtener los factores de estructura parciales de CuCl fundido. El resultado más destacado fue la ausencia de una estructura marcada en la distribución de los cationes. Se trata de una característica relevante porque, como han demostrado trabajos posteriores, es común a los haluros de plata y cobre pero no se observa en los haluros alcalinos. Esta diferencia parece estar relacionada con el hecho de que los haluros de plata y cobre se caracterizan por exhibir antes de fundir una conductividad iónica excepcionalmente alta debida a la difusión de los cationes, lo que se conoce como comportamiento superiónico. Las aplicaciones tecnológicas de las sales superiónicas abarcan ámbitos diversos, desde las baterías recargables de litio ligeras, utilizadas en telefonía móvil, ordenadores portátiles o dispositivos médicos, hasta pilas de combustible o sensores de gas. Introducción 3 En nuestro grupo de investigación, los primeros pasos en la simulación de sales se deben a los Dres. Antoni Giró y Joaquim Trullàs, que implementaron el cálculo de los potenciales de largo alcance mediante la técnica de las sumas de Ewald y confeccionaron la primera versión de nuestro programa de simulación. Su interés se centró desde el principio en los haluros de plata y cobre fundidos, y junto con los profesores Silbert y Stafford (de la Universidad de East Anglia en Norwich), llevaron a término las primeras simulaciones de algunos de estos sistemas (AgI, CuBr, CuCl y CuI) [9], en lo que significó el comienzo de una fructífera colaboración con el profesor Moisés Silbert que sigue hasta la actualidad y ha tenido un papel importante en esta tesis. Desde este primer trabajo se optó por simular la fase líquida de dichas sales con los potenciales propuestos por Vashishta y Rahman [10] para el estado sólido. Ya con la aportación de la Dra. Olga Alcaraz, el estudio de estos sistemas se completó con AgBr y AgCl [11]. A partir de este momento, nuestro grupo comenzó a investigar la importancia de la polarización inducida en las propiedades de las sales fundidas. Para ello se propusieron modelos de ión polarizable (PIM, de Polarizable Ion Model) en los que los iones se representan mediante cargas puntuales y dipolos puntuales inducibles. A diferencia de los haluros alcalinos, que ya se habían simulado con shell models a partir de mediados de la década de los 70, los trabajos acerca del efecto de la polarización en los haluros de plata no aparecieron hasta mucho después. Los dos primeros fueron los que los grupos de Madden [12] y Trullàs [13] dedicaron a AgCl fundido. Aproximadamente en este punto comenzó el trabajo que conforma esta tesis, bajo la dirección del Dr. Joaquim Trullàs. Los objetivos que nos planteamos al iniciarla fueron fundamentalmente dos. El primero consistía en estudiar el origen de la polarización catastrófica, que es una inestabilidad carente de sentido físico observada en las simulaciones de sales binarias fundidas en que se asume la polarizabilidad de ambas especies iónicas. La polarización catastrófica consiste en el crecimiento incontrolado de los momentos dipolares inducidos en los iones, que acaban divergiendo. El segundo propósito era analizar los efectos de la polarización inducida en AgBr y AgI fundidos, y completar así la serie iniciada con AgCl. Para ello simulamos estas sales con diversos PIM, y comparamos los resultados con los obtenidos con RIM. En principio pretendíamos calcular las propiedades microscópicas relativas a la estructura y el transporte iónico, pero ciertos 4 Introducción resultados que obtuvimos nos llevaron a estudiar también las propiedades dieléctricas de AgI y NaI fundidos. Lo que sigue es un breve resumen del contenido de los seis capítulos en que se organiza la tesis. Para poner en contexto nuestro trabajo, en el Capítulo 1 hacemos una introducción a las sales que nos ocupan y a los RIM que hasta ahora se han utilizado para simularlas. Comenzamos presentando las características de los haluros de plata y cobre, que comparamos con las de los haluros alcalinos. A continuación describimos los RIM que se han venido proponiendo para ambos tipos de sales, y recopilamos los resultados más importantes obtenidos con ellos, aprovechando para definir las propiedades estructurales y de transporte iónico que calculamos en nuestras simulaciones y a las que nos referimos durante la tesis. En el Capítulo 2 describimos los modelos en los que se asume que los momentos dipolares iónicos son inducidos únicamente por el campo eléctrico local, a los que llamamos PIMe (la letra e es la inicial de electric field). Resolvemos el problema teórico de la interacción entre dos iones aislados en este tipo de modelos, y el análisis de la amplia casuística que se da en función de si son polarizables o no nos permite justificar el origen de la polarización catastrófica. Para evitarla proponemos un nuevo tipo de modelos que, a la polarización inducida por el campo eléctrico, añaden aquella debida a la deformación de las capas electrónicas causada por su solapamiento con las capas electrónicas de los iones vecinos más próximos. Los denominamos modelos PIMs (donde la letra s, inicial de short range, se refiere a la inclusión de esta contribución de corto alcance a la polarización). En el Capítulo 3 presentamos los resultados obtenidos al simular AgBr fundido con un modelo PIMe y uno PIMs en los que sólo se polarizan los iones de bromo, y nos centramos en el efecto de la polarización en la estructura de la sal, sobre todo en las correlaciones entre iones de plata, y en el transporte iónico. La comparación con los datos experimentales revela que el PIMe es el modelo que permite reproducir mejor el factor de estructura estático coherente y la conductividad iónica. El Capítulo 4 es el dedicado a AgI, que hemos simulado con un RIM, un PIMe en el que sólo los aniones se polarizan y un PIMs con ambas especies polarizables. En primer lugar constatamos que los tres modelos reproducen correctamente el comportamiento superiónico de la fase α-AgI, en la que Introducción 5 únicamente los iones de plata se difunden. También simulamos la fase líquida y determinamos los cambios que se observan como consecuencia de la inducción de polarización. En este caso el acuerdo entre el factor de estructura estático coherente simulado con un PIMe y el experimental es especialmente satisfactorio. En el Capítulo 5 nos ocupamos de las propiedades dieléctricas. En primer lugar presentamos las expresiones teóricas de la función de respuesta dieléctrica y de la función dieléctrica longitudinal para modelos en los que coexisten cargas y dipolos puntuales. Para valorar las diferencias entre la respuesta dieléctrica en los PIM y en los RIM hemos hecho simulaciones de AgI y NaI fundidos, escogidos de entre los haluros de plata y alcalinos por la elevada polarizabilidad del iodo. Para ello hemos debido implementar en el programa de simulación el cálculo de las correlaciones entre las densidades de carga y de dipolo puntual. Finalmente, en el Capítulo 6 resumimos las principales conclusiones del trabajo expuesto en esta tesis, y anticipamos posibles perspectivas de continuación. Además añadimos seis apéndices que contienen datos o desarrollos teóricos de utilidad pero que interrumpirían la exposición si estuvieran incluidos en los capítulos que acabamos de resumir. En el Apéndice A recopilamos la notación básica utilizada a lo largo de la memoria, y el Apéndice B está dedicado al sistema de unidades cegesimal (CGS), que es el utilizado en todas las ecuaciones de la tesis. En el Apéndice C reunimos los valores de los parámetros de los potenciales utilizados o citados. El Apéndice D contiene el desarrollo matemático de las relaciones que definen los modelos PIMe, y el Apéndice E se ocupa brevemente de otro tipo de modelos no utilizados en nuestro trabajo, los basados en las ideas de Thole. Finalmente, en el Apéndice F detallamos la deducción del teorema de fluctuación-disipación de la función de respuesta dieléctrica para modelos con cargas y dipolos puntuales inducibles, utilizado en el Capítulo 5. El contenido de esta tesis ha sido publicado total o parcialmente en forma de artículos en revistas especializadas [14,15,16,17,18,19]. Además, durante el período de investigación publicamos otros trabajos sobre mezclas de sales fundidas [20,21,22], que no han sido incluidos en la tesis porque su contenido es ajeno al tema de la misma. Asimismo, durante la redacción de 6 Introducción esta memoria hemos realizado estudios acerca de los efectos de la polarización inducida en la mezcla eutéctica de AgCl y AgI [23] y en CuI fundido [24]. Las principales conclusiones de ambos trabajos se apuntan en el Capítulo 6. Referencias [1] B. J. Alder and T. E. Wainwright J. Chem. Phys. 31, 459 (1959). [2] A. Rahman Phys. Rev. A 136, 405 (1964). [3] L. Verlet Phys. Rev. 165, 201 (1967). [4] L. V. Woodcock and K. Singer Trans. Faraday Soc. 67, 12 (1971). [5] M. P. Tosi and F. G. Fumi J. Phys. Chem. Solids 25, 45 (1964). [6] M. J. L. Sangster and M. Dixon Adv. Phys. 25, 247 (1976). [7] N. H. March and M. P. Tosi, Coulomb Liquids, Academic, London (1984). [8] D. I. Page and K. Mika J. Phys. C 4, 3034 (1971). [9] A. J. Stafford, M. Silbert, J. Trullàs and A. Giró J. Phys.: Condens. Matter 2, 6631 (1990); J. Trullàs, A. Giró and M. Silbert J. Phys. Condens. Matter 2, 6643 (1990). Las Figuras 2 y 3 en la primera referencia deben ser intercambiadas para coincidir con los pies de figura. Además, en la Tabla 1 de este mismo trabajo el valor de α X para AgI debería ser 6.52 Å 3 y el de H XX para CuI debería ser 399 en unidades de e 2 Å 6 = 14.4 eVÅ 7 . [10] P. Vashishta and A. Rahman Phys. Rev. Lett. 40, 1337 (1978). [11] Ç. Tasseven, J. Trullàs, O. Alcaraz, M. Silbert and A. Giró J. Chem. Phys. 106, 7286 (1997). [12] M. Wilson, P. A. Madden and B. J. Costa-Cabral J. Phys. Chem. 100, 1227 (1996). [13] J. Trullàs, O. Alcaraz, L. E. González and M. Silbert J. Phys. Chem. B 107, 282 (2003). [14] V. Bitrián and J. Trullàs J. Phys. Chem. B 110, 7490 (2006). [15] V. Bitrián, J. Trullàs and M. Silbert J. Chem. Phys. 126, 021105 (2007). [16] O. Alcaraz, V. Bitrián and J. Trullàs J. Chem. Phys. 127, 154508 (2007). [17] V. Bitrián and J. Trullàs J. Phys. Chem. B 112, 1718 (2008). [18] V. Bitrián and J. Trullàs J. Phys.: Conf. Ser. 98, 042006 (2008). [19] V. Bitrián, O. Alcaraz and J. Trullàs J. Chem. Phys. 130, 234504 (2009). Introducción 7 [20] V. Bitrián, J. Trullàs, M. Silbert, T. Enosaki, Y. Kawakita and S. Takeda J. Chem. Phys. 125, 184510 (2006). [21] J. Trullàs, V. Bitrián, Y. Kawakita, S. Takeda, T. Enosaki and M. Silbert J. Non-Cryst. Solids 353, 2993 (2007). [22] V. Bitrián, J. Trullàs and M. Silbert Physica B 403, 4249 (2008). [23] O. Alcaraz, V. Bitrián and J. Trullàs J. Chem. Phys. 134, 014505 (2011). [24] V. Bitrián, O. Alcaraz and J. Trullàs J. Chem. Phys. 134, 044501 (2011). 8 Introducción Capítulo 1 Sales 1:1 fundidas y modelos de ión rígido Las sales 1:1, también llamadas monovalentes, están formadas por dos especies iónicas con cargas +1e (cationes) y -1e (aniones). En este capítulo resumimos algunas de las propiedades de las dos clases de sales 1:1 a las que nos referiremos a lo largo de la tesis: los haluros alcalinos y los haluros de plata y cobre. Gran cantidad de información acerca de las características de las sales 1:1 fundidas ha sido extraída de la simulación por dinámica molecular (MD, de Molecular Dynamics) de modelos de ión rígido (RIM, de Rigid Ion Model), que son aquellos en los que se considera que los iones interactúan mediante un potencial efectivo a parejas. Los dos tipos de potencial que se han utilizado más comúnmente son los de Born-Mayer-Huggins (BMH) para los haluros alcalinos, y los de Vashishta-Rahman (VR) para los haluros de plata y cobre. 16 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido dicen que su fusión, propiciada por la aparición de defectos de Schottky a temperaturas cercanas a T f , frustra la transición a la fase superiónica [7,9]. 1.2 Modelos de ión rígido para sales 1:1 Para poder aplicar la mecánica estadística a un sistema de partículas, ya sea a nivel teórico o mediante la simulación, es necesario disponer de un Hamiltoniano en el que se establece un modelo de interacción entre las partículas. En los Hamiltonianos más simples se supone que la energía potencial puede expresarse en términos de un potencial efectivo aditivo a parejas. En el caso de los sistemas iónicos estos potenciales se denominan de ión rígido y dan lugar a los modelos de ión rígido (RIM), llamados así para diferenciarlos de los modelos de ión polarizable (PIM, de Polarizable Ion Model), en los que se tiene en cuenta la polarización inducida en los iones y, a diferencia de los RIM, la interacción entre dos de ellos depende de los que los rodean. En este capítulo introductorio nos centramos en los RIM. En los RIM se asume que los iones se comportan como esferas blandas cargadas, y que la energía potencial de un sistema de N iones situados en las posiciones {r i } puede escribirse como RIM RIM 1 1 ( ) 2 N N ij ij i j i U r υ = ≠ = ∑∑ con 2 RIM shr ( ) ( ) i j ij ij ij ij ij z z e r r r υ υ = + , (1.1) donde r ij = |r ij | = |r i - r j | es la distancia interiónica y RIM ( ) ij ij r υ es un potencial aditivo a parejas. El primer término de RIM ( ) ij ij r υ , expresado en el sistema cegesimal de unidades (ver Apéndice B), corresponde a la interacción coulombiana, que es de largo alcance, con z i la carga del ión i en unidades de la carga fundamental e. De hecho, al asumir este término se está suponiendo implícitamente que cada ión tiene una distribución finita de carga que siempre es esféricamente simétrica, y que no se superpone con las de otros iones, de modo que el campo que crea en un punto exterior es igual al que crearía una carga puntual z i e situada en su centro, y la energía potencial de interacción electrostática del sistema es idéntica a la energía de un conjunto de cargas puntuales. Por su parte, shr ( ) ij ij r υ es un potencial efectivo de corto alcance que incluye un término característico de esferas blandas para modelizar la repulsión entre dos iones debida al solapamiento Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 17 de sus capas electrónicas externas. Esta repulsión, que a cortas distancias domina sobre todas las demás, es consecuencia del principio de exclusión de Pauli, que exige que los electrones de las capas solapadas ocupen estados cuánticos distintos. Las formas funcionales más habituales de este término son la de decaimiento exponencial ( B exp[- γ r] ) propuesta por Born y Mayer [10], y la de potencia inversa ( H/r n ) basada en ideas de Pauling [11]. Además, a menudo shr ( ) ij ij r υ contiene otras interacciones de corto alcance como las de van der Waals, cuyo origen es la existencia temporal de momentos multipolares de corta duración asociados con transiciones cuánticas de los iones [12]. 1.2.1 Potenciales de Born-Mayer-Huggins La gran mayoría de simulaciones de haluros alcalinos fundidos se han realizado con potenciales de ión rígido cuya contribución de corto alcance incluye el término exponencial de Born-Mayer y términos atractivos correspondientes a las interacciones de van der Waals. Estos potenciales, propuestos por primera vez por Huggins y Mayer [13], y a los que denominaremos potenciales de Born-Mayer-Huggins (BMH), se escriben [ ] 2 BMH 6 8 ( ) exp a b ab ab ab ab ab z z e C D r B r r r r υ γ = + − − − , (1.2) donde a y b indican las especies iónicas (+ o -), y z + = |z - | = 1. En el término exponencial γ ab determina la dureza de la pared repulsiva, y [ ] exp ( ) ab ab ab a b B b β γ σ σ = + , (1.3) donde σ a y σ b indican de algún modo los radios iónicos dado que su suma es una medida de la distancia de acercamiento entre iones a partir de la cual la repulsión comienza a aumentar significativamente, b es una constante común a todos los haluros alcalinos y tipos de interacción, y los factores de Pauling vienen dados por [11] 1 / / ab a a b b z n z n β = + + , (1.4) con n a el número de electrones en la última capa de los iones (igual a 8 para todos los iones haluros y alcalinos excepto el Li + , para el que es 2). Los términos proporcionales a r -6 y r -8 corresponden a las interacciones de van der Waals que dan cuenta de la interacción entre momentos dipolares, y 18 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido entre momentos dipolares y cuadripolares, inducidos por fluctuaciones cuánticas. Los parámetros C ab y D ab se pueden calcular con las aproximaciones de London [14], y fueron estimados por primera vez para los haluros alcalinos por Mayer [15]. Sus valores están tabulados en el Apéndice C. Utilizando los valores de Mayer para C ab y D ab , Fumi y Tosi [16] estimaron el resto de parámetros de los potenciales BMH para los diecisiete haluros alcalinos que cristalizan en la estructura de sal gema. Propusieron el mismo exponente para los tres tipos de interacción, γ = γ ++ = γ +- = γ -- , y permitieron que γ fuera distinto para cada sal. Por tanto, los parámetros a determinar eran los diecisiete γ y los nueve σ a de los iones que componen estas sales (cinco alcalinos y cuatro haluros). Para ello expresaron la energía cohesiva de la red cristalina en función de la energía de interacción entre cada ión y sus primeros y segundos vecinos, y calcularon los valores de γ y σ a que mejor ajustaban la ecuación de estado y su primera derivada respecto al volumen [17]. Los potenciales obtenidos como resultado de este ajuste son conocidos como potenciales de Fumi-Tosi (FT). Los valores de sus parámetros están recogidos en el Apéndice C. A partir de los años 70, los potenciales FT comenzaron a ser utilizados en simulaciones de haluros alcalinos fundidos. Pese a haber sido concebidos para describir la fase sólida, proporcionan buenos resultados de la energía interna en la fase líquida, reproducen bastante bien los factores de estructura extraídos de experimentos de difracción de neutrones, y permiten obtener resultados bastante aceptables de los coeficientes de difusión o la conductividad iónica [18]. Por lo que respecta a los haluros de metales nobles, se intentó también obtener potenciales BMH que reprodujeran sus características. Mayer y Levy utilizaron, como harían años más tarde Fumi y Tosi, las ecuaciones termodinámicas de la energía cohesiva del sólido para parametrizar los potenciales de los haluros de plata y talio [19] y de cobre [20]. Sin embargo, simulaciones posteriores han mostrado que tales potenciales no proporcionan resultados satisfactorios cuando se aplican a la fase líquida, ya que predicen una temperatura de fusión mayor que la experimental [21,22]. Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 19 1.2.2 Potenciales de Vashishta-Rahman Para estudiar el comportamiento superiónico de α-AgI [23] y α-CuI [24], Vashishta y Rahman propusieron un tipo de potencial de ión rígido, que llamaremos VR, con la siguiente forma funcional, 2 VR 4 6 ( ) ab a b ab ab ab ab n z z e H P C r r r r r υ = + − − , (1.5) donde, a diferencia de los potenciales BMH, se consideran cargas efectivas z = z + = |z - | < 1 que pretenden modelizar la incompleta ionización de las partículas e introducen un cierto grado de covalencia en las interacciones [25]. El segundo término, que modeliza la interacción repulsiva de corto alcance, está basado en ideas previas de Pauling [11] con ( ) ab n ab a b H A σ σ = + , (1.6) donde A depende de la sal pero no de las especies a y b, σ a y σ b están relacionados de nuevo con los radios iónicos, y la potencia n ab determina la dureza de la pared repulsiva. El término –P ab /r 4 es, como veremos en el Capítulo 2, una aproximación aditiva a parejas de la energía de interacción debida a la polarización inducida en los iones, con 2 2 2 ( ) 2 ab a b b a e P z z α α = + , (1.7) donde α a es la polarizabilidad electrónica de los iones de la especie a. La aproximación de London [14] para el primer término de van der Waals es 3 2 a b a b ab a b E E C E E α α =+, (1.8) donde E a es una determinada energía media de excitación. Dado que la polarizabilidad de los iones haluros es significativamente mayor que la de los de plata y cobre, tanto en los trabajos originales de Vashishta y Rahman como en los posteriores se ha tomado α + = 0 y α - ≠ 0, de modo que 0 P ++ = y 0 C C ++ +− = = . (1.9) Los parámetros de la ecuación (1.5) fueron determinados por Vashishta y Rahman para AgI y CuI [23,24] utilizando datos de la fase sólida. Por un lado, en base a resultados de la dispersión de fonones y de la 20 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido compresibilidad, estimaron z = 0.6 y n = n ++ = n +- = n -- = 7. Con este valor de z, además de los ya conocidos para la polarizabilidad α - y la energía E - , obtuvieron P +- , P -- y C -- mediante (1.7) y (1.8). Por otra parte, asignaron un valor a σ + y σ - suponiendo que los iones de mayor tamaño están “en contacto” en la fase sólida, y que lo mismo ocurre entre los iones de distinto signo. Así, calcularon σ Ag y σ I a partir de la distancia r 0 entre primeros vecinos en la fase α-AgI, y σ Cu y σ I a partir del valor de r 0 en la fase γ-CuI. Finalmente, determinaron A imponiendo la minimización de la energía de cohesión del sólido para la constante de red experimental. Los valores de los parámetros que Vashishta y Rahman obtuvieron para AgI y CuI están detallados en el Apéndice C. Estos potenciales, que llamaremos AgI-VR7 y CuI-VR7 (donde el número 7 hace referencia al valor de la potencia n), permitieron reproducir mediante simulaciones MD las fases α de AgI y CuI [23,24]. Posteriormente, Parrinello, Rahman y Vashishta [26] recalcularon para AgI los valores de los parámetros del potencial VR. Para ello utilizaron datos de la fase β-AgI en lugar de la α. Las principales novedades de este nuevo potencial, al que llamaremos AgI-PRV, respecto a AgI-VR7 son la mayor diferencia entre radios iónicos y que n ab ya no es igual para los tres tipos de interacción, sino que n ++ = 11, n +- = 9 y n -- = 7. Con él reprodujeron con éxito la transición de fase α↔β. Unos años después, Shimojo y Kobayashi [27] advirtieron que la densidad utilizada en la Referencia [26] era ligeramente menor que la experimental y retocaron los valores de los parámetros del potencial AgI-PRV para reproducir los resultados de Parrinello, Rahman y Vashishta a la densidad experimental. A este nuevo potencial lo llamaremos AgI-SK. Tanto AgI-PRV como AgI-SK están incluidos en el Apéndice C. A partir de experimentos de difracción de neutrones (NS, de Neutron Scattering), Howe et al. [28] observaron que la estructura de los iones Ag + en α-AgI presenta claras similitudes con la de los Cu + en el CuCl fundido [29,30]. Este resultado llevó a Stafford y Silbert [31] a sugerir que las sales superiónicas fundidas pueden presentar características similares a las que tienen en estado sólido, y que los potenciales VR podrían reproducir también su comportamiento líquido. Así, llevaron a cabo cálculos de las funciones de distribución radial y los factores de estructura del AgI fundido mediante la teoría HNC (Hypernetted-Chain) utilizando el potencial AgI- Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 21 PRV. Los resultados de experimentos NS posteriores de Takahashi et al. [32] confirmaron que tales cálculos reproducían cualitativamente algunos de los principales rasgos del factor de estructura de AgI. A partir de entonces Silbert y sus colaboradores han utilizado potenciales VR para estudiar los haluros de plata y cobre fundidos. Stafford y Silbert, siguiendo el mismo procedimiento de Vashishta y Rahman, parametrizaron los potenciales VR para CuCl y CuBr (CuCl-VR7 y CuBr-VR7 en el Apéndice C), y los utilizaron, junto a AgI-PRV y CuI-VR7, para realizar cálculos HNC de las propiedades estructurales de AgI, CuCl, CuBr y CuI fundidos [33]. Sus resultados de las funciones de distribución radial y de los factores de estructura parciales para CuCl coinciden bastante bien con los obtenidos a partir de experimentos NS con substitución isotópica [29,30] (en el momento del estudio no se disponía de resultados experimentales para CuBr y CuI). Con los mismos potenciales, Trullàs y Giró realizaron simulaciones MD de las funciones de distribución radial, que al coincidir salvo pequeños detalles con las g ab (r) teóricas confirmaron la validez de la aproximación HNC [33]. Estas simulaciones MD les permitieron, además, estudiar también las propiedades dinámicas [34]. En particular, obtuvieron una gran diferencia en la difusividad de cationes y aniones, lo cual concordaba con los resultados experimentales de Poignet et al. [35] para los coeficientes de autodifusión del CuCl fundido, que son los únicos datos experimentales disponibles de la difusión de estas sales. Finalmente, Tasseven et al. [22] parametrizaron los potenciales VR para AgCl y AgBr (AgCl-VR6 y AgBrVR6 en el Apéndice C) y llevaron a cabo sus correspondientes simulaciones MD en fase líquida. De este modo completaron el estudio de los haluros de plata y cobre mediante potenciales de ión rígido iniciado por Stafford y Silbert. En este punto es importante remarcar que la diferencia fundamental de los potenciales VR respecto a los BMH no es su forma funcional distinta, sino los valores de sus parámetros. En particular, las cargas efectivas z < 1 y una diferencia de tamaños considerablemente más acentuada. Se ha sugerido que estos rasgos buscan, en cierto sentido, imitar los efectos de la polarización inducida [21]. También hay que resaltar que, mientras que en los potenciales FT el valor de σ a asociado al radio de cada ión es el mismo para todos los haluros alcalinos de los que forma parte, en los VR el valor de σ a para un mismo ión puede cambiar de una sal a otra. 22 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido En los siguientes apartados veremos que los potenciales VR permiten reproducir de forma cualitativa algunas de las principales características estructurales y de transporte iónico de la mayoría de haluros de plata y cobre fundidos. Sin embargo, como veremos en el apartado 1.3.4, incluso en aquellos sistemas para los que estos potenciales funcionan mejor, como es el caso de CuCl o CuBr, sigue habiendo detalles de los factores de estructura experimentales que no se logra reproducir. Para otras sales como CuI, los potenciales VR conducen a resultados insatisfactorios, muy alejados de los experimentales [36]. 1.3 Propiedades estructurales de sales 1:1 fundidas 1.3.1 Funciones de distribución radial En un líquido que ocupa un volumen V y está formado por s especies con N a partículas cada una, cada especie tiene una densidad parcial ρ a = N a /V y una concentración c a = N a /N, donde N = Σ a N a es el número total de partículas, y la densidad total de partículas es ρ N = N/V = Σ a ρ a . Las propiedades estructurales del líquido se describen a partir de la densidad local de partículas de cada especie, ρ a (r), y de sus fluctuaciones alrededor del valor medio, δ ρ a (r), que vienen dadas por 1 ( ) ( ) a N a ia ia ρ δ = = − ∑ r r r y ( ) ( ) a a a ρ ρ ρ δ = − r r , (1.10) donde r ia es la posición de una partícula de la especie a y δ (r) es la delta de Dirac. Las transformadas de Fourier de ρ a (r) y δ ρ a (r) para el vector de onda k son 1 ( ) exp( · ) a N a ia ia i ρ = = − ∑ k k r y ,0 ( ) ( ) a a a N ρ ρ δ δ = − k k k , (1.11) donde δ k,0 es la delta de Kronecker que vale 1 si k = 0 y se anula en cualquier otro caso. Las correlaciones entre las densidades locales de partículas de dos especies a y b se pueden describir mediante la siguiente función de correlación de densidades parciales, Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 23 N N ( ) ( ) 1 ( ) ( ) ( ) a b ab ab a b ab a b G c c g c c ρ ρ δ δ ρ ρ = = + r 0 r r r , (1.12) donde los corchetes 〈…〉 denotan el promedio estadístico. En sistemas homogéneos 〈 ρ a (R+r) ρ b (R)〉 es igual para todos los puntos R, en particular R = 0, y por tanto 〈 ρ a (r) ρ b (0)〉 = (1/V)∫〈 ρ a (R+r) ρ b (R)〉 dR. El término δ ab δ (r) incluye la autocorrelación de cada partícula consigo misma (contribución self), y la función de distribución de parejas g ab (r) da cuenta de las correlaciones entre partículas distintas (contribución distinct). En sistemas isotrópicos como las sales fundidas, g ab (r) sólo depende de r = |r|, y g ab (r) = g ab (r). Entonces g ab (r) se denomina función de distribución radial. La densidad media de partículas de la especie a situadas a una distancia r respecto a una partícula b distinta, ρ ab (r), está relacionada con g ab (r) por ρ ab (r) = ρ a g ab (r). El límite de g ab (r) para r→∞ es igual a 1. Figura 1.4: Funciones de distribución radial de NaCl fundido a T = 1170 K (izquierda) y CuCl fundido a T = 773 K (derecha), obtenidas mediante simulaciones MD con los potenciales NaCl-FT y CuCl-VR7: g ++ (r) (curva continua), g +- (r) (curva de puntos), g -- (r) (curva de trazos). 0 2 4 6 8 10 12 r (Å) 0 1 2 3 4 g ab ( r ) 2 4 6 8 10 12 14 r (Å) NaCl NaNa NaCl ClCl CuCl CuCu CuCl ClCl 24 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido Excepto en algunos pocos casos en los que han sido extraídas de forma indirecta de experimentos NS, las g ab (r) han sido obtenidas mediante cálculos teóricos o, más a menudo, simulaciones. En la Figura 1.4 se muestran los resultados de las simulaciones MD de NaCl y CuCl fundidos con sendos RIM. El principal rasgo que comparten todas las sales 1:1 fundidas simuladas utilizando RIM es la estructura de capas alternadas de diferente signo alrededor de un ión, que se manifiesta en el hecho de que g +- (r) está en oposición de fase respecto a g ++ (r) y g -- (r). Cada ión está rodeado por una primera capa de iones de signo contrario situados aproximadamente a la distancia del primer pico de g +- (r). La segunda capa, a distancias similares a las de los primeros picos de g ++ (r) y g -- (r), la forman iones del mismo signo que el central. Sin embargo, mientras que en los haluros alcalinos g ++ (r) ≈ g -- (r), en los haluros de plata y cobre los aniones muestran una estructura mucho más marcada que los cationes, siendo los picos de g -- (r) más altos y estrechos que los de g ++ (r). Además, la penetración de cationes en la primera esfera de coordinación de otro catión es más significativa en los haluros de plata y cobre que en los alcalinos. Como se aprecia en la Figura 1.4, la subida inicial de g CuCu (r) se adentra en el primer pico de g CuCl (r), cosa que no ocurre con g NaNa (r) respecto a g NaCl (r). Estas diferencias son debidas fundamentalmente a que los iones de plata y cobre son considerablemente más pequeños que los iones haluros. 1.3.2 Factores de estructura parciales Las transformadas de las funciones de correlación entre fluctuaciones de densidades locales, (c a c b ) -1/2 〈 δρ a (r) δρ b (0)〉/ ρ N = G ab (r) - ρ N (c a c b ) 1/2 , son los factores de estructura parciales de Ashcroft-Langreth [37], ( ) ( ) 1 ( ) a b ab a b SN c c ρ ρ δ δ − = k k k , (1.13) que haciendo uso de las ecuaciones (1.11) también se pueden escribir Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 25 ,0 1 1 1 1 ( ) exp( · ) exp( · ) a b N N ab ia jb a b ia jb a b S i i N c c N c c δ = = = − − ∑ ∑ k k k r k r . (1.14) La expresión (1.14) se reduce a S ab (k) = (c a c b ) -1/2 〈 ρ a (k) ρ b (-k)〉/N para k ≠ 0, es decir a la transformada de Fourier de G ab (r). Los factores de estructura parciales trasladan al espacio recíproco la información proporcionada por las g ab (r), con las que están relacionadas a través de N ( ) F ( ) ( ( ) 1) ab ab a b ab S c c g δ δ ρ   = + −   k r r , (1.15) donde F indica la transformada espacial de Fourier. En sistemas isotrópicos los factores de estructura son funciones de k = |k|, S ab (k) = S ab (k), y la relación (1.15) se convierte en 2 N0 sin( ) ( ) ( ( ) 1) 4 πd ab ab a b ab kr S k c c g r r r kr δ ρ ∞ = + − ∫ , (1.16) donde δ ab es la contribución self, es decir la de los términos ia = jb en los sumatorios de (1.14), y es la única que sobrevive a k → ∞. Conviene indicar que algunos autores [38] definen S ab ( k ) sin el factor (c a c b ) -1/2 . La inclusión de este factor hace que S ab (k) → δ ab cuando k → ∞. En las simulaciones MD con condiciones periódicas de contorno, los vectores del espacio recíproco están restringidos por la condición k = 2π n /L, con n un vector de componentes enteras y L el lado de la celda cúbica de simulación [39]. En las simulaciones realizadas en esta tesis, para k > 4 Å -1 hemos calculado los S ab (k) mediante (1.16). Sin embargo, para valores de k pequeños, los efectos espurios del tamaño finito de la celda de simulación en la cola de g ab (r) se vuelven relevantes, y los S ab (k) calculados presentan una oscilación artificial. Por ello, para k < 4 Å -1 hemos obtenido los S ab (k) a partir de la expresión (1.14), adecuadamente promediada sobre los n k vectores de onda que tienen el mismo módulo | k | = k , 1 1 1 1 1 ( ) exp( · ) exp( · ) a b N N ab ia jb k ia jb ka b S k i i n N c c = = =     = −       ∑ ∑ ∑ k k r k r . (1.17) 32 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido Figura 1.8: Factores de estructura S c (k) de NaCl a 1148 K, KCl a 1073 K, RbCl a 1023 K y CsCl a 968 K, obtenidos mediante experimentos NS (curvas negras continuas) y simulaciones MD con los potenciales FT (curvas rojas de puntos). Para obtener S c (k) a partir de un experimento NS es necesario sustraer de la intensidad medida aquella que no es debida a la colisión de los neutrones con la muestra, como la de fondo o la dispersada por el contenedor de la muestra [43, 49]. También es necesario sustraer la contribución de los neutrones que han sido dispersados más de una vez antes de ser detectados (corrección de multiple scattering [50]), y tener en cuenta que una fracción de los neutrones incidentes han sido absorbidos tanto por el contenedor como por la muestra [51]. Además, no se pueden ignorar los efectos de la inelasticidad. A pesar de que la mayoría de colisiones pueden considerarse aproximadamente elásticas, existe un intercambio de energía entre la 0.4 0.8 1.2 1.6 S c ( k ) 012345 k (Å -1 ) NaCl 123456 k (Å -1 ) KCl 0.0 0.4 0.8 1.2 S c ( k ) RbCl CsCl NS [45] MD NS [46] MD NS [47] MD NS [48] MD Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 33 muestra y algunos neutrones, lo que complica la relación teórica entre la intensidad medida y el factor de estructura (1.30) y exige la aplicación de ciertos términos correctivos que dependen en gran medida del montaje experimental [52,53]. En la Figura 1.8 se comparan los S c (k) experimentales de diversos cloruros alcalinos fundidos con los obtenidos a partir de simulaciones MD de modelos RIM. La forma de los S c (k) simulados es similar a la de los experimentales, y la posición de los picos y valles, así como la amplitud de las oscilaciones posteriores al pico principal, es aproximadamente la misma en ambos. Se puede afirmar, pues, que los modelos de ión rígido reproducen bastante bien la estructura de los haluros alcalinos fundidos. Las discordancias más destacables surgen al comparar las alturas de algunos picos y, excepto para CsCl, el comportamiento de S c (k→0). De hecho, los valores de los S c (k) experimentales a k→0 no son fiables ya que, a juzgar por la κ T de los haluros alcalinos fundidos [54] y el límite (1.28), S c (k→0) debería valer ∼ 0.1, aproximadamente como en las simulaciones. En la Figura 1.9 se representan los S c (k) de los haluros de plata y cobre fundidos. Las simulaciones de CuCl, y en menor medida las de CuBr, reproducen bastante bien los resultados experimentales. En ambos casos la posición de los dos primeros máximos coincide con la determinada experimentalmente, aunque hay diferencias en la altura. Para CuCl se reproducen también las oscilaciones posteriores a los dos primeros picos, pero las de CuBr se amortiguan demasiado rápidamente. Sin embargo, el comportamiento a pequeñas k del S c (k) simulado no concuerda con el observado experimentalmente. En particular, no se logra reproducir el shoulder que estas dos sales exhiben en k ∼ 1 Å -1 . Por su parte, las simulaciones de CuI proporcionan resultados alejados de los experimentales. A pesar de que la posición de los picos en k ∼ 1.7 Å -1 y k ∼ 3.0 Å -1 concuerda, sus alturas son muy diferentes, y algo parecido ocurre con las oscilaciones a k grandes. Además, tampoco se logra reproducir el pico inicial en k ∼ 1 Å -1 (al que llamaremos prepico en capítulos posteriores). 34 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido Figura 1.9: Factores de estructura S c (k) de los haluros de plata y cobre fundidos, obtenidos mediante experimentos NS (curvas negras continuas y de puntos) y simulaciones MD con modelos RIM (rojas de trazos). En la leyenda de cada gráfica se indica la referencia de los datos experimentales y la temperatura correspondiente. Los potenciales utilizados en las simulaciones son AgCl-VR6, AgBr-VR7, AgI-SK, CuCl-VR7, CuBr-VR7 y CuI-VR7. 123456 k (Å -1 ) CuCl CuBr 0.4 0.8 1.2 1.6 S c ( k ) CuI 012345 k (Å -1 ) AgCl 0.4 0.8 1.2 S c ( k ) 0.0 0.4 0.8 1.2 S c ( k ) AgBr AgI NS [55] (773 K) NS [56] (773 K) MD (773 K) NS [55] (753 K) NS [56] (773 K) MD (753 K) NS [55] (933K) NS [56] (923K) MD (923K) NS [57] (778K) NS [30] (773K) NS [57] (853K) NS [58] (788K) MD (853K) MD (773 K) NS [57] (938 K) NS [36] (933 K) MD (938 K) Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 35 Los S c (k) de AgCl y AgBr obtenidos mediante experimentos NS tienen un pico principal con forma de altiplano, poco puntiagudo y muy ancho, en el que se aprecia un ligero rizado de tres picos. Las simulaciones no logran reproducir este rasgo particular. En su lugar muestran dos picos con un pronunciado mínimo entre ellos, el segundo de los cuales (en k ∼ 3.5 Å -1 ) está desplazado a k mayores respecto al experimental. Por lo que respecta a los resultados MD para AgI, la posición del primer pico en k ∼ 1.75 Å -1 coincide con el shoulder experimental que marca el inicio del pico principal, pero el máximo absoluto en las simulaciones está ligeramente desplazado hacia la derecha respecto al experimental. El prepico en k ∼ 1 Å -1 tampoco está presente en los resultados MD. Vale la pena remarcar finalmente que, como queda patente en la Figura 1.9, existen discrepancias entre los resultados de diferentes experimentos para una misma sal, ya sea debido a los diversos montajes experimentales o a las correcciones realizadas durante el análisis de los datos. Los resultados obtenidos para k pequeñas son los que presentan una mayor incertidumbre, ya que el ruido de fondo es especialmente intenso en la región de ángulos θ pequeños (es decir, k pequeñas), lo que afecta a la precisión de S c (k→0). Además, en ocasiones los valores de S c (k→0) se renormalizan imponiendo el cumplimiento de la relación (1.28). Por todo ello, conviene cierta prudencia a la hora de extraer conclusiones de resultados experimentales. 1.4 Transporte iónico en sales 1:1 fundidas 1.4.1 Autodifusión A diferencia de lo que ocurre en un sólido, los iones de un líquido se difunden. El coeficiente de autodifusión D a cuantifica esta propiedad, y se define a través de la conocida relación de Einstein 2 ( ) 1 lim [ ( ) (0)] lim 6 2 a a ia ia t t Q t D t t t →∞ →∞ = − =r r , (1.32) donde Q a (t) es el desplazamiento cuadrático medio 2 1 1 1 ( ) [ ( ) (0)] 3 a N a ia ia ia a Q t t N= = − ∑ r r , (1.33) 36 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido y r ia (t) es la posición de un ión cualquiera de la especie a en el instante t. Q a (t) está relacionado con la función de autocorrelación de velocidades, 1 1 1 ( ) ( )· (0) 3 a N a ia ia ia a t t N = Λ = ∑ v v , (1.34) donde v ia (t) es la velocidad de un ión a en t. Λ a (t) tiene un valor inicial Λ a (0) = k B T/m a , con m a la masa iónica de la especie a, y suele presentarse normalizada por este valor, B ( ) ( ) ( ) (0) a a a a a t m C t t k T Λ = = Λ Λ. (1.35) Q a (t) y Λ a (t) están relacionadas por Λ a (t) = (1/2) (d 2 Q a (t)/dt 2 ) [59], lo que permite obtener la siguiente relación de Green-Kubo para D a , B 0 0 ( )d ( )d a a a a k T D t t C t t m ∞ ∞ = Λ = ∫ ∫ . (1.36) En las Figuras 1.10 y 1.11 se representan las funciones Ca(t) y Qa(t) obtenidas para NaI y AgI a partir de simulaciones MD de modelos RIM. Cuando se describe el comportamiento dinámico de los iones se distinguen tres escalas de tiempo distintas. Para valores muy pequeños de t, o régimen de partícula libre, los iones se mueven a velocidad prácticamente constante e igual a la inicial, de modo que Qa(t) ∼ (kBT/ma) t 2 y Ca(t) ∼ 1. A medida que las interacciones del ión con sus vecinos propician la desviación de su trayectoria, Qa(t) se aparta del comportamiento cuadrático y toma valores menores. Esta escala temporal coincide con la región en la que Ca(t) presenta el decaimiento inicial y las oscilaciones posteriores, y se conoce como régimen cinético o de dinámica molecular. Finalmente, a t suficientemente grandes, o régimen hidrodinámico, cuando se han producido una serie de colisiones sucesivas, la velocidad del ión está descorrelacionada con la inicial y Ca(t) = 0, mientras que el desplazamiento cuadrático medio es el típico de un movimiento browniano a tiempos grandes, es decir lineal con t. De acuerdo con la relación de Einstein (1.32), la pendiente de Qa(t) para valores de t suficientemente grandes es 2Da, y por tanto Qa(t) ∼ 2Da t. Así pues, Da es una propiedad macroscópica característica del régimen hidrodinámico. Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 37 Figura 1.10: Desplazamientos cuadráticos medios de NaI a T = 1011 K y AgI a T = 923 K, obtenidos a partir de simulaciones MD con los potenciales NaI-FT y AgI-SK: Q + (t) (curvas continuas) y Q - (t) (curvas a trazos). En los ejes insertados en la esquina superior izquierda se representan las mismas funciones en la región de t pequeños. Las curvas finas corresponden al comportamiento límite de Q a (t) en el régimen de partícula libre (ver texto). Figura 1.11: Funciones de autocorrelación de velocidades normalizadas de NaI y AgI fundidos, calculadas en las simulaciones MD referidas en el pie de la Figura 1.10: C + (t) (curvas continuas) y C - (t) (curvas de trazos). 0 1 2 t ( p s ) 0 2 4 6 Q a ( t ) (Å 2 ) 123 t ( p s ) 0.0 0.1 0.2 0.0 0.1 0.2 0.0 0.1 0.2 0.0 0.1 0.2 NaI AgI 0.0 0.2 0.4 t (ps) -0.4 0.0 0.4 0.8 C a ( t ) 0.2 0.4 0.6 t (ps) NaI C Na ( t ) C I ( t ) AgI C Ag ( t ) C I ( t ) 38 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido Tabla 1.3: Coeficientes de autodifusión D + y D - (en 10 -5 cm 2 /s), conductividad iónica σ (en Ω -1 cm -1 ) y parámetro ∆ de diversas sales 1:1 fundidas, calculados en simulaciones MD de modelos RIM con los potenciales detallados, a las densidades ρ N (en Å -3 ) y temperaturas T (en K) indicadas. Las cantidades entre paréntesis corresponden a resultados experimentales de D a de haluros alcalinos [60] y CuCl [35], y de σ y ∆ [61]. Potencial ρ N T D + D - σ ∆ NaCl NaCl-FT 0.0312 1148 8.5 (10.7) 7.8 (7.8) 3.5 (3.8) 0.14 (0.19) KCl KCl-FT 0.0244 1073 6.4 (7.2) 6.5 (6.4) 2.7 (2.2) 0.01 (0.23) RbCl RbCl-FT 0.0221 1023 3.6 (5.1) 4.0 (4.5) 1.3 (1.6) 0.13 (0.17) CsCl CsCl-FT 0.0196 973 3.1 (4.0) 3.6 (4.4) 1.2 (1.3) 0.04 (0.17) NaI NaI-FT 0.0216 1011 7.71 (8.5) 5.56 (4.8) 2.7 (2.5) -0.14 (0.05) AgCl AgCl-VR6 0.0406 773 4.20 2.21 4.4 (4.0) -0.40 AgBr AgBr-VR7 0.0354 753 3.86 2.00 3.3 (3.0) -0.30 AgI AgI-SK 0.0281 923 6.51 2.21 3.4 (2.5) -0.36 CuCl CuCl-VR7 0.0443 773 9.38 (6.4) 2.66 (2.1) 9.8 (3.7) -0.53 (0.18) CuBr CuBr-VR7 0.0357 853 12.08 3.69 8.4 -0.37 CuI CuI-VR7 0.0281 938 9.50 2.51 5.0 -0.30 En la Tabla 1.3 se detallan los coeficientes de autodifusión en diversas sales 1:1 fundidas. Los efectos del tamaño de los iones en su difusividad prevalecen sobre los de la masa [62]. Así, los iones de menor tamaño poseen un coeficiente de autodifusión mayor que el de los iones más grandes. Por ejemplo, en la Figura 1.10 se aprecia claramente, a partir de la relación de Einstein, que D Na > D I y D Ag > D I , respectivamente. Al igual que en NaI y AgI, en la mayoría de haluros alcalinos los iones de menor radio son también los más ligeros, pero en sales como AgCl, AgBr o CuCl ocurre al revés, y la especie más pesada es la más pequeña y la que presenta una mayor difusividad. La diferencia más notable entre los haluros alcalinos y los haluros de plata y cobre es que, mientras que en los primeros D + y D - son similares, en los haluros de plata y cobre D + es entre dos y cuatro veces Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 39 mayor que D - . Este hecho es consecuencia de la mayor diferencia de tamaño entre aniones y cationes en los haluros de plata y cobre. Parece como si en estas sales fundidas se conservara, en cierto modo, la diferencia de comportamiento que cationes y aniones muestran en las fases superiónicas (en las que D + ≠ 0 y D - = 0). La diferencia entre los valores de los coeficientes de autodifusión calculados por las relaciones de Einstein y Green-Kubo es del orden del error estadístico, que suele ser de alrededor del ~ 1%. Los D a obtenidos en las simulaciones de los haluros alcalinos son en todos los casos ligeramente inferiores a los valores experimentales. Por otra parte, en el caso de CuCl, que es el único haluro de plata o cobre del que se dispone de datos experimentales ocurre lo contrario, y los D a calculados exceden a los experimentales. Los casos de KCl y RbCl son interesantes porque en ellos la especie con el coeficiente de autodifusión mayor no es la misma en los resultados de las simulaciones (Cl) que en los experimentales (K y Rb respectivamente). La razón de este comportamiento opuesto es que, mientras en los potenciales FT utilizados en las simulaciones el tamaño de los iones de cloro es menor que el de los de potasio y rubidio, el radio iónico estimado empíricamente por Pauling [63] es menor para K y Rb que para Cl. Hay que aclarar que cuando decimos que en un cierto potencial el tamaño de los iones de una especie es menor que el de los de otra nos referimos a que el potencial de interacción a pares entre los primeros es menos repulsivo que aquel entre los segundos, y les permite acercarse a distancias menores. Este tamaño no viene determinado sólo por el valor de σ a , sino también por el de los demás parámetros del potencial. Por eso, a pesar de que en los potenciales FT σ K < σ Cl , el tamaño de los iones de Cl es menor que el de los de K porque RIM ClCl υ es menos repulsivo que RIM KK υ . La masa iónica, por su parte, juega un papel principal en el régimen de partícula libre. En este régimen los desplazamientos cuadráticos medios de los iones más ligeros toman valores mayores debido a su menor inercia, como se observa en la Figura 1.10. En aquellas sales en las que los iones de menor masa son también los más grandes (AgCl, AgBr y CuCl), los dos Q a (t) se intersecan, ya que el de los iones ligeros comienza siendo mayor 40 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido que el de la otra especie, y acaba siendo menor debido a su menor coeficiente de autodifusión, como se aprecia en la Figura 1.12 para AgCl. Figura 1.12: Desplazamientos cuadráticos medios Q Ag (t) (curvas continuas) y Q Cl (t) (curvas de trazos) de AgCl a T = 773 K, extraídos de simulaciones MD con el potencial AgCl-VR6. Por lo que respecta al régimen cinético, las C a (t) son las funciones más adecuadas para describir la dinámica individual. Uno de sus rasgos más distintivos es el primer mínimo que presentan tras el decaimiento inicial. El hecho de que sea negativo indica que cuando ha transcurrido un intervalo de tiempo igual al valor de t en que se encuentra el mínimo, los iones han colisionado, en promedio, una vez con un ión vecino y su velocidad ha cambiado de sentido. A este cambio de sentido se le llama backscattering. El ejemplo de las C a (t) de NaI fundido en la Figura 1.11 sirve para ilustrar una característica común a todos los haluros alcalinos: la C a (t) de la especie más ligera (Na en este caso) decae más rápido y muestra un backscattering más pronunciado que la de la especie más pesada [64]. Este comportamiento se atribuye al movimiento vibratorio que los iones ligeros llevan a cabo en el interior de la capa formada por sus primeros vecinos más pesados [38,61]. Por el contrario la C a (t) de los iones pesados (como I en NaI) presenta un backscattering menos profundo y a tiempos mayores. Debido a su mayor 0 1 2 t (ps) 0 1 2 3 Q a ( t ) (Å 2 ) 0.0 0.2 0.4 0.0 0.2 0.4 AgCl Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido 41 inercia, la velocidad de estos iones cambia de dirección menos bruscamente, y su movimiento entre dos colisiones sucesivas se ve menos afectado por las interacciones con sus vecinos. En cambio, la C I (t) en AgI, pese a que los iones de I son más pesados que los de Ag, presenta un backscattering más pronunciado que C Ag (t). Este comportamiento se explica de la siguiente forma. En los haluros de plata y cobre fundidos los iones grandes (aniones) llevan a cabo un movimiento vibratorio en el interior de una doble capa formada por los primeros vecinos de signo opuesto, más pequeños, y los segundos vecinos de su misma especie. Por ello, C - (t) muestra un comportamiento marcadamente oscilatorio, con un mínimo de backscattering pronunciado. En cambio, los iones pequeños (cationes) se difunden a través de los huecos de la estructura compacta aniónica, de modo que C + (t) tiene un decaimiento más lento y un backscattering muy poco marcado. Por tanto, la diferencia de tamaños es el principal factor que explica la diferencia en la dinámica individual de estos sistemas en el régimen cinético, y las citadas características de las C a (t) son independientes de que los aniones sean los iones más ligeros, como ocurre en AgCl, AgBr y CuCl, o los más pesados como en AgI, CuBr y CuI, si bien es cierto que cuanto más ligeros son los aniones más pronunciado es su backscattering. Un detalle interesante a tener en cuenta es que una C a (t) marcadamente oscilatoria, a pesar de que recuerde al comportamiento de los iones en un sólido, no implica necesariamente una menor difusividad. De hecho, en el caso de NaI fundido, y en general en la mayoría de haluros alcalinos, ocurre al contrario. Si bien el área bajo la curva de C Na (t) es menor que bajo C I (t), la menor masa de los iones de Na hace que el área bajo Λ Na (t) = (k B T/m Na )C Na (t) sea mayor que bajo Λ I (t), y por tanto D Na > D I . 1.4.2 Conductividad La conductividad iónica σ se puede calcular mediante expresiones análogas a las relaciones de Einstein y de Green-Kubo para D a [38], 2 2 N NZZ ZZ B B 0 ( ) lim ( )d 2 t e e Q t t t k T t k T ρ ρ σ ∞ →∞ = = Λ ∫ . (1.37) En la relación de Einstein para σ , en el segundo miembro de (1.37), 48 Capítulo 1. Sales 1:1 fundidas y modelos de ión rígido [49] M. A. Howe, R. L. McGreevy and W. S. Howells J. Phys.: Condens. Matter 1, 3433 (1989). [50] I. A. Blech and B. L. Averbach Phys. Rev. A137, 1113 (1965). [51] H. H. Paalman and C. J. Pings J. Appl. Phys. 33, 2635 (1962). [52] G. Placzek Phys. Rev. 86, 377 (1952). [53] J. G. Powles Molec. Phys. 26, 1325 (1973). [54] G. J. Janz, Molten Salts Handbook, Academic, New York (1967). [55] M. Inui, S. Takeda, Y. Shirakawa, S. Tamaki, Y. Waseda and Y. Yamaguchi J. Phys. Soc. Japan 60, 3025 (1991). [56] Y. Kawakita, T. Enosaki, S. Takeda and K. Maruyama J. Non-Cryst. Solids 353, 3035 (2007). [57] Y. Shirakawa, M. Saito, S. Tamaki, M. Inui and S. Takeda J. Phys. Soc. Japan 60, 2678 (1991). [58] D. A. Allen and R. A. Howe J. Phys.: Condens. Matter 4, 6029 (1992). [59] U. Balucani and M. Zoppi, Dynamics of the Liquid State, Clarendon, Oxford (1994). [60] R. E. Young and J. P. O’Connell Ind. Eng. Chem. Fundam. 10, 418 (1971). [61] G. J. Janz, F. W. Dampier, G. R. Lakshminarayanan, P. K. Lorenz and R. P. T. Tomkins, Molten Salts, National Bureau of Standards Reference Data Series Nº 15, Washington DC (1968). [62] O. Alcaraz and J. Trullàs J. Chem. Phys. 113, 10635 (2000). [63] L. Pauling, The Nature of the Chemical Bond, Ithaca, Cornell (1945). [64] G. Ciccotti, G. Jacucci and I. R. McDonald Phys. Rev. A 13, 426 (1976). [65] J. Trullàs and J. A. Padró Phys. Rev. B 55, 12210 (1997). [66] J. P. Hansen and I. R. McDonald Phys. Rev. A 11, 2111 (1975). En la Figura 4 de este trabajo, así como en la Figura 10.4 de [38], se representa – δ ZZ (t). Capítulo 2 Modelos de ión polarizable En el capítulo anterior hemos visto que los modelos de ión rígido (RIM) proporcionan una descripción razonablemente buena de muchas de las propiedades de los haluros alcalinos fundidos. No ocurre lo mismo, sin embargo, en el caso de los haluros de plata y cobre, cuyas simulaciones con RIM no reproducen algunos de los rasgos más distintivos de los factores de estructura experimentales, como el triple rizado del pico ancho de AgCl y AgBr, o los prepicos de AgI y los haluros de cobre. Por ello es necesario modelizar las interacciones iónicas de una forma más realista que en los RIM. Un primer paso en esta dirección consiste en tener en cuenta que la distribución de carga de un ión no tiene por qué ser siempre esférica como se supone en los RIM, sino que puede deformarse en función de su entorno, es decir polarizarse. Los modelos de sistemas iónicos en los que se tiene en cuenta este efecto son los modelos de ión polarizable (PIM, de Polarizable Ion Model). La mayoría de los modelos en los que se tiene en cuenta la polarización inducida pueden clasificarse en tres tipos [1]: los Shell Models, los modelos de carga fluctuante, y los de dipolo puntual inducido. En este capítulo, después de una breve introducción a los tres, describimos en detalle los modelos de dipolo puntual inducido que hemos utilizado en este trabajo. Comenzaremos por los modelos en los que se asume que los momentos dipolares son inducidos únicamente por el campo eléctrico local, a los que denominamos PIMe (la letra e, inicial de electric field, hace referencia al hecho de que la polarización es inducida sólo por el campo eléctrico). En ellos puede producirse una polarización excesiva de los iones, fenómeno 50 Capítulo 2. Modelos de ión polarizable que se conoce como polarización catastrófica. Una manera de evitarla es añadir, a la polarización inducida por el campo eléctrico, aquélla debida a la deformación de las capas electrónicas causada por su solapamiento con las capas electrónicas de los iones vecinos más próximos, lo que da lugar a los modelos que denominamos PIMs (la letra s, inicial de short range, denota la inclusión de esta contribución de corto alcance a la polarización). 2.1 Tipos de modelos de ión polarizable 2.1.1 Shell Models Los primeros PIM que se propusieron fueron los Shell Models (SM) [2], que, entre muchos otros sistemas sólidos y líquidos [1], han sido utilizados para estudiar los haluros alcalinos fundidos [3-10]. En los SM cada ión polarizable i se describe mediante dos cargas puntuales de signo opuesto y valor constante conectadas por una fuerza recuperadora armónica de constante k i , de modo que su distancia de separación d i puede variar. El valor de la carga positiva es z i e+Q i y el de la negativa –Q i , de modo que la carga iónica neta es z i e. En cierto modo, un modelo de este tipo se justifica físicamente asumiendo que la carga positiva representa el núcleo atómico y las capas electrónicas internas, y que la negativa representa el centro de la capa electrónica externa. El ión se polariza cuando el centro de la capa externa se desplaza respecto al núcleo. La energía potencial electrostática de un SM es la de interacción entre los pares de cargas de todos los iones. A esta hay que añadir la energía interna de polarización de cada ión, o sea la necesaria para polarizarlo, 2 /2 i i k d . También hay que sumar la energía de corto alcance, es decir la repulsión y la dispersión de van der Waals, que al ser debidas a interacciones entre capas externas se modelizan en función de las distancias entre cargas negativas. En las simulaciones MD de los SM se asigna una masa nula a las capas externas, de modo que estas se adaptan instantáneamente a cualquier cambio en la configuración de los núcleos y ocupan las posiciones que minimizan la energía del sistema. Por otro lado las posiciones de los núcleos, en los cuales se localiza la masa iónica y sobre los cuales actúa la fuerza resultante de cada ión, se hacen evolucionar integrando numéricamente sus ecuaciones del movimiento. Capítulo 2. Modelos de ión polarizable 51 2.1.2 Modelos de carga fluctuante Los modelos de carga fluctuante (FCM, de Fluctuating Charge Model) están muy extendidos en el estudio de sistemas moleculares. En un FCM cada molécula se representa mediante cargas situadas en las posiciones atómicas (y en ocasiones también en posiciones no ocupadas por ningún átomo), cuyo valor varía en función de su entorno. Tales variaciones están sujetas a la condición de conservación de la carga molecular neta y se calculan imponiendo que la energía del sistema sea mínima en todo momento. Los FCM han sido ampliamente utilizados en el estudio del agua líquida. Los dos modelos más extendidos para el agua [11] consisten en tres cargas puntuales situadas en las posiciones del oxígeno y los dos hidrógenos (Simple Point Charge), o en dos cargas situadas en los dos hidrógenos y una tercera en un punto de la bisectriz del ángulo HOH que no coincide con la posición del átomo de oxígeno (Four-point Transferable Intermolecular Potential). Los cuatro puntos a los que hace referencia el nombre de este último modelo son las tres cargas y la posición del oxígeno, ya que se asume que las interacciones de corto alcance intermoleculares tienen lugar entre los átomos de oxígeno. En lo concerniente a las sales 1:1, Wilson y Madden [12] utilizaron una variante de FCM para estudiar las propiedades estructurales y dinámicas de NaCl y LiF, tanto en fase sólida como líquida. En ella modelizaron cada ión, además de con una carga constante z i e, con dos cargas opuestas variables equidistantes del centro del ión, situadas en los extremos de una varilla rígida rotatoria. En este caso, a la energía potencial electrostática debida a la interacción coulombiana entre todas las cargas hay que sumar la energía interna de polarización, así como la debida a las interacciones efectivas de corto alcance entre pares de iones. 2.1.3 Modelos de dipolo puntual inducido En el caso de sistemas iónicos, los modelos de dipolo puntual inducido (PPDM, de Polarizable Point Dipole Model) se basan en un RIM. La novedad respecto a los RIM es que cada ión se representa, además de con una carga puntual, con un dipolo puntual inducido. Entonces, a las interacciones de corto alcance y las coulombianas entre cargas puntuales 52 Capítulo 2. Modelos de ión polarizable propias del RIM, hay que añadir aquellas entre cargas y dipolos y entre dipolos puntuales, así como la energía interna de polarización, es decir las relacionadas con la inducción de polarización. El primer trabajo de simulación de sales 1:1 con un PPDM fue el de Wilson et al. sobre las fases sólida y líquida de AgCl (modelos PIM1 y PIM2 en [13], basados en el potencial BMH de Mayer [14]). Más tarde Trullàs et al. [15] llevaron a cabo un estudio de esta misma sal fundida con un PPDM basado en el potencial de ión rígido AgCl-VR6 (ver Apéndice C). Los PPDM, al igual que los SM y el FCM de Wilson y Madden [12], asumen que las interacciones electrostáticas pueden ser descritas de forma adecuada a nivel dipolar, obviando momentos de orden superior. Algunos autores han incluido también cuadripolos puntuales [13], pero sus efectos en las propiedades de las sales 1:1 fundidas no parecen significativos, aunque no ocurre lo mismo en los sólidos. De entre los diferentes tipos de PIM citados en este apartado, en esta tesis nos hemos limitado a estudiar los PPDM. A pesar de ser un tipo particular de PIM, nos referiremos a ellos con esta denominación genérica. En los apartados que siguen definiremos las diferentes versiones que hemos utilizado. 2.2 Modelo de iones con dipolos puntuales inducidos por el campo eléctrico (PIMe) En los modelos de iones con dipolos puntuales inducidos, a cada ión i de los N que forman el sistema se le asocia una carga puntual z i e y un dipolo puntual de momento dipolar µ µµ µ i inducido por todos los demás iones. La elección más simple es considerar el campo eléctrico como única causa de la polarización, y suponer la siguiente relación constitutiva q ( ) i i i i i i α α µ = = + µE E E , (2.1) donde α i es la polarizabilidad electrónica del ión i y E i el campo eléctrico local en la posición r i creado por el resto de iones j (j ≠ i), que es la suma del campo creado por el resto de cargas puntuales z j e, Capítulo 2. Modelos de ión polarizable 53 q q N i ij j i≠ = ∑ E E con q 3 j ij ij ij z e r = E r , (2.2) más el campo creado por el resto de dipolos µ µµ µ j , N i ij j i µ µ ≠ = ∑ E E con 5 3 · 3 j ij j ij ij ij ij r r µ = − µrµ E r . (2.3) Como se aprecia al substituir (2.3) en (2.1), el momento dipolar de un ión i depende de los momentos µ µµ µ j del resto de iones j ≠ i, y µ µµ µ j depende a su vez de los momentos de los demás iones, entre ellos µ µµ µ i . Las N relaciones vectoriales (2.1) forman pues un sistema de 3N ecuaciones lineales, cuya resolución implica un aumento de la complejidad de las simulaciones respecto a las de los RIM. Nos referimos a los modelos en los que se aplica la relación (2.1) con el acrónimo PIMe (la letra e hace referencia al hecho de que la polarización es inducida sólo por el campo eléctrico). Como se demuestra en el Apéndice D, la energía potencial total de un PIMe es PIMe RIM ind(e) U U U= + , (2.4) donde el primer miembro es la energía potencial de un RIM con potencial efectivo de ión rígido RIM ij υ , RIM RIM 1 1 ( ) 2 N N ij ij i j i U r υ = ≠ = ∑∑ con 2 RIM shr ( ) ( ) i j ij ij ij ij ij z z e r r r υ υ = + , (2.5) y U ind(e) es la energía debida a la inducción de polarización, que viene dada por 2 ind(e) q µ 1 1 1 1 1 · · 2 2 N N N i i i i i i i i i U µ α = = = = − − + ∑ ∑ ∑ µEµE. (2.6) El primer término en (2.6) es la energía de interacción entre cargas y dipolos puntuales (U qµ en el Apéndice D), el segundo la de interacción entre dipolos puntuales (U µµ ), y el último es la energía interna que los iones adquieren al polarizarse (U self ). 54 Capítulo 2. Modelos de ión polarizable El potencial efectivo a parejas RIM ( ) ij ij r υ que aparece en (2.5) suele ser, en el estudio de sales 1:1 fundidas, de los tipos BMH o VR citados en el Capítulo 1. En el caso de los PIM basados en potenciales VR se omite el término –P ij /r 4 , que no es más que una aproximación a parejas de la expresión (2.6) de U ind(e) , de modo el potencial de ión rígido que se asume en estos modelos es de la forma 2 VR0 6 ( ) ij i j ij ij ij ij n ij ij ij z z e H C r r r r υ = + − . (2.7) A los potenciales con esta forma funcional los llamamos potenciales VR0. A diferencia de lo que ocurre con U RIM , U ind(e) no se puede expresar en términos de un potencial aditivo a parejas. Introduciendo las expresiones de los campos (2.2) y (2.3) en (2.6) es posible escribir U ind(e) como la suma de las interacciones entre todas las parejas de iones i y j, como se ve en la ecuación (D.55) del Apéndice D. Sin embargo, la interacción de una cierta pareja no depende sólo de ambos iones, ya que los momentos µ µµ µ i y µ µµ µ j dependen a su vez de todas las demás partículas. Para deducir U PIMe es necesario asumir la hipótesis de adiabaticidad, es decir que en cada instante y para cada configuración los dipolos puntuales inducidos están en equilibrio con su entorno. En el Apéndice D, partiendo de la relación constitutiva (2.1) e imponiendo la adiabaticidad, hemos deducido que la energía interna de polarización (U self ) es la del último término en (2.6). La adiabaticidad implica que los valores de los momentos dipolares dados por (2.1) son los que minimizan U PIMe , y que por tanto debe cumplirse PIMe ind(e) 0 i i U U = = µ µ µ µµ µ µ µ ∇ ∇ ∇ ∇∇ ∇ ∇ ∇ (∀ i). (2.8) Esta relación es fundamental en el cálculo de las expresiones de las fuerzas, detalladas en el Apéndice D. Podríamos haber procedido como hacen otros autores [12] que, en lugar de definir un PIMe a partir de la relación constitutiva, lo definen directamente estableciendo U PIMe y, a continuación, deducen la relación constitutiva imponiendo la adiabaticidad mediante (2.8). En las simulaciones de sistemas macroscópicos se utilizan condiciones periódicas de contorno para evitar los efectos de superficie. Las Capítulo 2. Modelos de ión polarizable 55 interacciones de corto alcance modelizadas por shr ( ) ij ij r υ se calculan mediante el convenio de la mínima imagen [16], que sólo tiene en cuenta las interacciones entre aquellas partículas cuya separación es menor de una cierta distancia de truncamiento inferior o igual a la mitad del lado de la celda de simulación. En cambio, en el cálculo de los campos eléctricos, la energía potencial U ind(e) y las fuerzas electrostáticas, no se puede aplicar dicho truncamiento porque las interacciones electrostáticas son de largo alcance y en el resultado influyen también los iones de las celdas muy alejadas. Estas contribuciones se calculan mediante el método de las sumas de Ewald [6,16]. Las expresiones de las sumas de Ewald necesarias para simular un RIM o un PIMe están recogidas por ejemplo en el trabajo de Trullàs et al. sobre AgCl fundido [15]. 2.3 Cálculo de los momentos dipolares inducidos en un PIMe La simulación de los PIMe requiere el cálculo de los momentos dipolares que los iones adquieren en cada configuración del sistema. Para ello se pueden utilizar tres métodos, que difieren entre sí en su exactitud y en el consumo de tiempo de cálculo computacional que requieren. 2.3.1 Método exacto de inversión de la matriz Introduciendo (2.2) y (2.3) en (2.1), los N momentos µ µµ µ i vienen dados por las N relaciones vectoriales siguientes 3 5 3 3( · ) N j ij j ij ij j i i j i ij ij ij z e r r r α ≠   = + −       ∑ r µr r µ µ ( ) i ∀ . (2.9) Este sistema de 3N ecuaciones puede rescribirse en forma matricial como [1] q =· A µE , (2.10) donde A es una matriz 3N×3N, y µ = (µ µµ µ 1 ,…,µ µµ µ N ) y q q q 1 ( ,..., ) N = E E E son vectores cuyas 3N componentes son las de los momentos dipolares y los campos q i E respectivamente. La matriz A se compone de N 2 submatrices 3×3, 2 [1 ] ( ) ij ij ij ij i δ δ α = − − r A I T , (2.11) 56 Capítulo 2. Modelos de ión polarizable siendo I la matriz identidad, δ ij la delta de Kronecker y T 2 el tensor dipolar definido en el Apéndice D, cuyas componentes son 2 5 3 ( ) 3r r T r r β γ βγ βγ δ = −r, (2.12) donde los subíndices β y γ indican las componentes del vector r = r ij . Haciendo explícitos todos los elementos de las matrices y las componentes de los vectores, la ecuación (2.10) se escribe 3 q 1 1 ( ) N ij j i j A E βγ γ β γ µ = = = ∑∑ con 2 ( ) [1 ] ( ) ij ij ij ij i A T βγ βγ βγ δ δ δ α = − − r . (2.13) Entonces, dado que A y E q dependen sólo de las posiciones iónicas, una vez conocidas estas en una determinada configuración los momentos se obtienen a partir de la inversión de A en (2.10), -1 q · = µA E . (2.14) En simulaciones con N ∼ 1000 iones, el rango de A es ∼ 10 3 y su inversión requiere una gran cantidad de tiempo computacional, por lo que es necesario aplicar otros métodos menos precisos pero más rápidos para calcular los momentos, como el procedimiento iterativo de predicción y corrección. 2.3.2 Método iterativo de predicción y corrección En la gran mayoría de simulaciones llevadas a cabo en este trabajo hemos utilizado un procedimiento iterativo de predicción y corrección [17]. En este procedimiento primero se calcula el conjunto de los campos { q i E }, que sólo dependen de las posiciones iónicas, y se hace una predicción inicial de los N momentos, {µ µµ µ i(P) }. Lo más sencillo es escoger los valores del paso de tiempo anterior, o si se trata del primer paso en una simulación {µ µµ µ i(P) = 0}. A continuación se utilizan estos momentos dipolares predichos para calcular mediante la ecuación (2.3) los correspondientes { µ(P) i E }, y se comienza el siguiente proceso iterativo: (i) a partir de { µ(P) i E } se calculan los momentos dipolares corregidos, (C) q µ(P) ( ) i i i i α = +µE E . (2.15) Capítulo 2. Modelos de ión polarizable 57 (ii) Los momentos {µ µµ µ i(C) }, que sólo coinciden con los {µ µµ µ i(P) } si son la solución que buscamos, se introducen en (2.3) para calcular los correspondientes { µ(C) i E }. (iii) Los campos corregidos { µ(C) i E } pasan a ser los nuevos { µ(P) i E } en el paso (i), y el ciclo iterativo se repite hasta que se cumple 2 µ(C) µ(P) 2 µ(C) ∆ i i i − ≤ E E E ( ) i ∀ , (2.16) donde el parámetro de convergencia ∆ es usualmente del orden de 10 -4 o 10 -6 . Finalmente, los valores de los momentos dipolares que se toman son los calculados mediante (2.15) con los { µ(C) i E } finales que satisfacen (2.16). La figura 2.1 muestra un ejemplo típico del proceso de convergencia de las tres componentes del momento dipolar de un anión cualquiera en una simulación MD de una sal 1:1 fundida utilizando un PIMe en el que sólo los aniones son polarizables. El número de iteraciones necesarias en este tipo de simulaciones oscila entre 10 y 20 si la precisión es ∆ ∼ 10 -6 , o entre 5 y 10 si ∆ ∼ 10 -4 . Figura 2.1: Componentes del momento dipolar de un anión polarizable en unidades de e Å (puntos unidos con rectas) calculadas mediante el método de predicción y corrección en una simulación de una sal 1:1 fundida utilizando un PIMe en el que sólo los aniones son polarizables. Los valores obtenidos por el método exacto de inversión de la matriz se indican con rectas horizontales discontinuas. 0 4 8 12 16 i t e r a c i o n e s 0.609 0.610 0.611 0.612 µ y µ z + 0.34 µ x + 0.69 64 Capítulo 2. Modelos de ión polarizable En el caso de un catión y un anión aislados con sólo uno de ellos polarizable, por ejemplo α 1 = α + = 0 y α 2 = α - ≠ 0 como se supone en los potenciales VR para los haluros de plata y cobre, r C = 0 y υ ind(e) (r) 2 2 4 /2 z e r α − + = − coincide con el termino –P +- /r 4 de los potenciales VR, cuya inclusión en el potencial de ión rígido viene a ser por lo tanto una aproximación a parejas de U ind(e) en (2.6). Como se observa en el ejemplo de la Figura 2.4, a pesar de que υ ind(e) (r→0) → -∞ la mayor intensidad de la repulsión de corto alcance hace que υ PIMe (r) tome valores positivos a cortas distancias. Figura 2.4: (a) υ RIM (r), υ ind(e) (r) y υ PIMe (r), y (b) módulo del momento dipolar µ µµ µ 2 , para un PIMe en el caso de dos iones aislados, un catión de Ag con α 1 = 0 y un anión de Br con α 2 = 4.16 Å 3 , basado en el potencial de ión rígido AgBr-VR7 (ver Apéndice C). En el caso de dos iones polarizables del mismo tipo, de la especie a por ejemplo, y por tanto z 1 = z 2 = z a y α 1 = α 2 = α a , estamos ante un caso particular en el que se cumple 2 2 1 2 2 1 z z α α = . Teniendo en cuenta que en el denominador común de las ecuaciones (2.19)-(2.21) aparece r 6 -4 2 a α = (r 3 +2 α a )(r 3 -2 α a ), y que el numerador pasa a ser proporcional al factor (r 3 -2 α a ), estas ecuaciones pueden escribirse como 012345 r (Å) -4 -2 0 2 υ (eV) 012345 r (Å) 0 2 4 µ (e Å) (a) (b) µ 2 υ PIMe υ ind(e) υ RIM Capítulo 2. Modelos de ión polarizable 65 1 2 3 ˆ ( ) ( ) ( 2 ) a a a z e r r r r α α = = − + µrµ y 2 2 ind(e) 3 ( ) ( 2 ) a a a z e rr r α υ α = − +, (2.24) donde la singularidad en r C = (2 α a ) 1/3 ha desaparecido, como se visualiza en el ejemplo de la Figura 2.5. Vale la pena subrayar que, a diferencia de lo que ocurre en el caso de un anión polarizable y un catión que no lo es, υ ind(e) (r) en (2.24) no coincide, salvo para valores grandes de r, con el término –P aa /r 4 = - α a z a2 e 2 /r 4 presente en la interacción anión-anión de los potenciales VR. Figura 2.5: (a) υ RIM (r), υ ind(e) (r) y υ PIMe (r), y (b) módulos de los momentos dipolares µ µµ µ 1 y µ µµ µ 2 , para un PIMe en el caso de dos aniones aislados de Br, con α 1 = α 2 = 4.16 Å 3 , basado en el potencial de ión rígido AgBr-VR7 (ver Apéndice C). A partir de los resultados (2.24) se podría concluir que la polarización catastrófica no puede ocurrir en los PIMe para sales binarias con sólo una especie polarizable. Tal conclusión es precipitada porque el resultado anterior únicamente es válido para iones aislados, mientras que en una sal fundida hay que tener en cuenta el papel fundamental que juega la influencia de los demás iones, especialmente en el caso de la interacción entre dos iones del mismo signo ya que entre ellos habrá otros de signo contrario que también contribuirán a su polarización. Para dos iones aislados, la condición 2 2 1 2 2 1 z z α α = es una condición de equilibrio inestable 012345 r (Å) -8 0 8 16 υ (eV) 1 2 3 4 5 r (Å) 0 2 µ (e Å) (a) (b) υ PIMe υ ind(e) υ RIM µ 1 µ 2 66 Capítulo 2. Modelos de ión polarizable ya que cualquier mínima variación en las polarizabilidades hace que la singularidad en r C reaparezca. Del mismo modo, cualquier diferencia entre los entornos respectivos de dos iones iguales inmersos en un sistema de N iones puede romper el equilibrio que tendrían si estuvieran aislados, provocando la polarización catastrófica. En efecto, la polarización catastrófica tiene lugar también en modelos PIMe en los que únicamente se inducen momentos dipolares en los aniones, como ocurre para NaCl o NaI fundidos [25]. Las simulaciones MD de ambas sales muestran que en ciertas configuraciones, y como consecuencia de interacciones con terceras partículas, dos aniones se aproximan a distancias cercanas a r C . A partir de entonces sus momentos dipolares adquieren sentidos paralelos y comienzan a crecer de forma anómala. En la Figura 2.6 hemos representado los momentos dipolares de los dos Cl - que originan la polarización catastrófica del NaCl en un caso particular cualquiera. Presentan una dependencia con la distancia muy similar a la de los iones polarizables del mismo signo pero de especies distintas de la Figura 2.3. Vale la pena subrayar que la distancia a la que se produce la divergencia es aproximadamente r C = (2 α Cl ) 1/3 , de modo que parece que la relación (2.23) sigue siendo válida en estos casos. Figura 2.6: Proyección en la dirección de ˆ r de los momentos inducidos en los dos iones Cl - que originan la polarización catastrófica en las simulaciones de NaCl fundido con un PIMe basado en el potencial NaCl-FT (ver Apéndice C), durante los pasos previos a la polarización catastrófica. Ambos momentos están prácticamente sobre la dirección de ˆ r , y los valores negativos de una de las curvas indican que a esas distancias son antiparalelos. La línea vertical a puntos indica r C = (2 α Cl ) 1/3 . 1 2 3 4 r (Å) 0 5 10 µ (e Å) NaCl Cl - Cl - r Capítulo 2. Modelos de ión polarizable 67 2.5 Modelos de iones con polarización de corto alcance (PIMs) Los PIM en los que los iones sólo se polarizan por el campo eléctrico dejan de ser realistas cuando los iones se encuentran a distancias típicas del solapamiento entre sus capas electrónicas puesto que, como hemos visto, puede producirse la polarización catastrófica. La corrección de esta anomalía se puede abordar desde dos puntos de vista; el propuesto originalmente por Thole [28], o el desarrollado paralelamente por Wilson y Madden [12,13,29] y por Tosi y sus colaboradores [30-33]. En ambos casos, pese a partir de supuestos distintos, se llega a una corrección del momento dipolar inducido debida a las interacciones a cortas distancias entre iones vecinos. La corrección es muy parecida en los dos casos y ha sido corroborada mediante cálculos ab initio [29,34,35]. En los modelos de Thole [36,35] a cada ión se le asocia, en lugar de una carga puntual, una distribución extensa de carga con simetría esférica, y en lugar de un dipolo puntual, otras dos distribuciones también esféricamente simétricas de signo opuesto, que de alguna manera pretenden dar cuenta de la amplitud espacial de la densidad electrónica. Cuando la separación entre dos iones es grande, su interacción se describe igual que en un PIMe. Sin embargo, cuando se encuentran a cortas distancias se produce una intersección entre sus distribuciones de carga que provoca un apantallamiento de los campos eléctricos y las energías de interacción que evita la polarización catastrófica. Los fundamentos de los modelos de Thole están resumidos en el Apéndice E. En el enfoque alternativo propuesto por Wilson y Madden se considera que las distribuciones de carga electrónica de los iones no se pueden superponer, pero sí que se pueden deformar como consecuencia de las interacciones de corto alcance con los vecinos más cercanos. Esta deformación, que no es de origen electrostático sino mecánico (al menos desde un punto de vista de física clásica), da lugar a la que llamamos polarización de corto alcance. La misma idea está presente en el modelo propuesto por Tosi y sus colaboradores, que resulta de llevar un Shell Model al límite de ión puntual. En ambos casos la ecuación constitutiva para el momento dipolar de un ión se puede escribir como s i i i i α = + µEµ , (2.25) 68 Capítulo 2. Modelos de ión polarizable de modo que al término inducido por el campo eléctrico se le añade la contribución s i µ debida a las interacciones de corto alcance que deforman el ión. 2.5.1 Modelos PIMsq Tanto Wilson y Madden como Tosi y colaboradores suponen que los efectos de las interacciones polarizadoras de corto alcance son aditivos a parejas y radiales, de modo que la relación constitutiva de s i µ puede escribirse de forma general como s ˆ ( ) N i i ij ij ij j i s r α ≠ = − ∑ µr , (2.26) donde ˆ ij r = r ij /r ij y s ij (r) es una función de corto alcance de la distancia r entre iones. Mientras Tosi et al. [30] usan una s ij (r) proporcional a la fuerza de corto alcance que j ejerce sobre i, Wilson y Madden [12] sugieren q 2 ( ) ( ) j ij ij z e s r f r r =, (2.27) donde q 4 q q 0 ( / ) ( ) exp( / ) ! k ij ij ij k r f r r k ρ ρ = = − ∑ (2.28) es la función propuesta originalmente por Tang y Toennies en su trabajo sobre el amortiguamiento de las interacciones de van der Waals [37]. La función q ( ) ij f r es decreciente, vale 1 cuando r = 0 y es significativamente diversa de 0 sólo a valores de r pequeños. De este modo s i µ amortigua las contribuciones a µ µµ µ i debidas a los campos eléctricos q 2 ˆ ( / ) ij j ij íj z e r= E r creados por las cargas puntuales de los vecinos más próximos al ión i. A los modelos definidos por las ecuaciones (2.25)-(2.28) los denominamos PIMsq. La letra s, inicial de short range, indica la inclusión de la polarización de corto alcance, y q hace referencia a que la contribución del ión j a s i µ es proporcional a q ij E y de sentido opuesto. Los PIMsq han sido utilizados por Wilson, Madden y sus colaboradores para estudiar diferentes sales 1:1 [13], 1:2 [18,19,29,38] y 1:3 [39] fundidas. Capítulo 2. Modelos de ión polarizable 69 A pesar del amortiguamiento de los momentos dipolares introducido en los PIMsq, este no siempre garantiza que se evite la polarización catastrófica, como ocurre, por ejemplo, en nuestros trabajos sobre AgBr [26] y AgI [45] fundidos. Por este motivo en algunos trabajos más recientes Madden y sus colaboradores [21,40-44] han optado por añadir un parámetro constante c ij > 1 que multiplica a la función q ( ) ij f r en (2.27), y que hace que el amortiguamiento sea más intenso. Alternativamente, en esta tesis hemos desarrollado otro tipo de modelos, que describimos a continuación. 2.5.2 Modelos PIMsqµ µµ µ En los modelos que llamamos PIMsqµ la ecuación (2.25) sigue siendo válida, pero ahora s i µ se modeliza mediante la siguiente fórmula general s q q µ µ [ ( ) ( ) ] N i i ij ij ij ij ij ij j i f r f r α ≠ = − + ∑ µE E , (2.29) con q ij E y µ ij E definidos en (2.2) y (2.3). La función de corto alcance q ( ) ij f r es la misma que en la expresión (2.28), y µ 4 µ µ 0 ( / ) ( ) exp( / ) ! k ij ij ij k r f r r k ρ ρ = = − ∑ . (2.30) Es decir, la polarización de corto alcance también amortigua las contribuciones inducidas por el campo eléctrico creado por los momentos dipolares de los iones cercanos. La relación (2.29) , a diferencia de la (2.26), no es aditiva a parejas porque en ella aparecen los momentos dipolares. El valor de X ij ρ , con X = {q,µ}, en (2.28) y (2.30) es una medida de la escala de longitud característica del amortiguamiento de los momentos dipolares. La distancia interiónica a la cual el amortiguamiento se reduce a la mitad, X X ( ) 1/2 ij ij f λ = , es X X 4.671 ij ij λ ρ = . (2.31) Como X ij λ tiene una interpretación más intuitiva que X ij ρ , la utilizaremos usualmente para definir los modelos con que trabajemos. En el límite 70 Capítulo 2. Modelos de ión polarizable qµ 0 ij ij λ λ = = se recupera la ecuación constitutiva de un PIMe, y el caso qµ ij ij λ λ = = ∞ corresponde a un RIM, ya que los momentos se anulan. Por su parte, los PIMsq se recuperan tomando q 0 ij λ ≠ y µ 0 ij λ = . En la Tabla 2.1 se introduce la nomenclatura que utilizamos a lo largo de la tesis para los nombres de los PIM de sistemas binarios. En los modelos en cuyo nombre aparece un 1 sólo los aniones son polarizables, mientras que en los que hay un 2 también lo son los cationes; las letras tras el número indican si se trata de un PIMe, PIMsq o PIMsqµ. En todos los casos, la definición del modelo requiere proporcionar las polarizabilidades α + y α - y las longitudes características de amortiguamiento X ij λ . Por simplicidad, en esta tesis hemos asumido la misma longitud característica para las tres interacciones en sales binarias, X X X X λ λ λ λ ++ +− −− = = = , y en los PIMsqµ hemos tomado q µ λ λ = . Tabla 2.1: Nomenclatura de los PIM utilizados en la tesis. α + α - q λ µ λ PIM1e 0 ≠ 0 0 0 PIM2e ≠ 0 ≠ 0 0 0 PIM1sq 0 ≠ 0 ≠ 0 0 PIM2sq ≠ 0 ≠ 0 ≠ 0 0 PIM1sqµ µµ µ 0 ≠ 0 ≠ 0 = λ q PIM2sqµ µµ µ ≠ 0 ≠ 0 ≠ 0 = λ q 2.5.3 Energía y fuerzas en los PIMs Puesto que en un PIMs las densidades de carga no se intersecan entre sí, se puede aplicar la aproximación dipolar para distribuciones disjuntas detallada en el Apéndice D, y la energía de interacción electrostática entre iones es la misma que en un PIMe. También sigue siendo válida la misma expresión para la energía interna de polarización (U self ) que aparece en (2.6). De Capítulo 2. Modelos de ión polarizable 71 hecho, esta expresión es válida para todos aquellos modelos en los que los momentos dipolares son inducidos linealmente por el campo eléctrico, independientemente de que haya otras posibles causas de polarización, como es el caso de la relación constitutiva (2.25). Esto es así porque el trabajo necesario para polarizar un ión debe ser el mismo con independencia de la causa de la polarización, y por tanto igual al obtenido en la ecuación (D.48) del Apéndice D suponiendo que sólo actúa el campo eléctrico. Obviamente la energía potencial de un PIMs, U PIMs , también incorpora los términos shr ( ) ij ij r υ propios del RIM en el que se fundamenta. Por tanto, U PIMs contiene todos los términos de U PIMe , pero además incluye un término, que denotaremos por U s , que da cuenta de las interacciones responsables de la polarización de corto alcance, es decir PIMs PIMe s U U U = + , (2.32) que teniendo en cuenta (2.4) también se puede escribir PIMs RIM ind(s) U U U= + , (2.33) donde ind(s) ind(e) s U U U = + (2.34) es la energía debida a la inducción de polarización en un PIMs. La forma de U s se deduce imponiendo la condición de adiabaticidad, es decir que los momentos dipolares dados por la relación constitutiva (2.25) sean los de equilibrio y minimicen por tanto U PIMs . Si calculamos el gradiente de la energía (2.32) y exigimos que sea cero, llegamos a s PIMe s / / i i i i i i i U U α α = − = − + = − µ µ µEµ∇ ∇ ∇ ∇∇ ∇ ∇ ∇ , (2.35) donde para obtener el tercer miembro se ha utilizado la expresión (2.4) de U PIMe con la (2.6) de U ind(e) , y para pasar al último miembro se ha hecho uso de (2.25). Entonces, a partir de (2.29) se concluye que la expresión de U s que satisface (2.35) para todo i es s sq s µ U U U = + , (2.36) con sq q q 1 · ( ) N N i ij ij ij i j i U f r = ≠ = ∑ ∑ µE y s µ µ µ 1 1 · ( ) 2 N N i ij ij ij i j i U f r = ≠ = ∑ ∑ µE . (2.37) Esta expresión de U s asegura, por tanto, que en cada instante de tiempo y en cualquier configuración se cumple 72 Capítulo 2. Modelos de ión polarizable PIMs ind(s) 0 ( ) i i U U i = = ∀ µ µ ∇ ∇ ∇ ∇∇ ∇ ∇ ∇ . (2.38) La fuerza que un ión k experimenta en un sistema de N iones se obtiene calculando el gradiente negativo de la energía potencial (2.32), y se puede escribir como PIMs PIMe s k k k = + F F F , (2.39) donde PIMe k F es la fuerza de un ión en un PIMe, cuya expresión es la de la ecuación (D.61) del Apéndice D, y s sq s µ k k k = + F F F (2.40) es la fuerza debida a la polarización de corto alcance, con sq q qµq 4 4 · · ( ) ( ) N N j kj k kj k kj kj kj kj kj k kj j kj j k j k kj kj f r f r z e z e r r ≠ ≠   = − + −       ∑ ∑ µrµr F F r r ɺ (2.41) y µ sµ µ µµ µ · ( ) ( ) N N k kj k kj kj kj kj kj kj j k j k kj f r f r r ≠ ≠ = − − ∑ ∑ µE F F r ɺ , (2.42) donde el punto sobre las funciones f kj indica la derivada respecto a r kj , y q µ µ µ kj k kj j jk z e z e = − F E E (2.43) y µµ 7 5 5 5 ( · )( · ) 3( · ) 3( · ) 3( · ) 15 k kj j kj k kj j j kj k j k kj kj kj kj kj kj kj r r r r = − + + + µrµrµrµ µ rµ µ µ r F r (2.44) son las contribuciones del ión j a las fuerzas q µ k F y µµ k F sobre el k definidas en el Apéndice D. Como en los PIMe, la condición de adiabaticidad (2.38) implica que la energía interna de polarización no da lugar a un término correspondiente de la fuerza, ya que en su expresión no aparecen explícitamente las posiciones iónicas. La implementación computacional de un PIMs a partir de un PIMe es muy sencilla. Los nuevos términos en las expresiones de los momentos dipolares, la energía potencial y la fuerza, es decir s i µ , U s y s k F , son simplemente correcciones de corto alcance ya que contienen las funciones fij(r). Por lo tanto, en las simulaciones MD de los PIMs todas estas nuevas Capítulo 2. Modelos de ión polarizable 73 contribuciones se calculan mediante el convenio de la mínima imagen. Las sumas de Ewald utilizadas para las interacciones de largo alcance son las mismas que en un PIMe. 2.6 Polarización catastrófica en un PIMs La resolución analítica del problema de dos iones aislados en un PIMs permite establecer si estos modelos son efectivos a la hora de evitar la polarización catastrófica. Usando la misma notación que en el apartado 2.4, los momentos dipolares inducidos en un par de iones asilados son 6 q µ 12 2 2 1 1 1 12 6µ2 2 5 1 2 12 ( ) 2 ˆ ( ) ( ) ( 4 [ ( )] ) r u r z e z eu r r u r r r α αα α   = −   −   µr r , (2.45) 6 q µ 12 1 1 2 2 2 12 6µ2 2 5 1 2 12 ( ) 2 ˆ ( ) ( ) ( 4 [ ( )] ) r u r z e z eu r r u r r r α αα α   = − −   −   µr r , (2.46) donde q q 12 12 ( ) 1 ( ) u r f r = − y µ µ 12 12 ( ) 1 ( ) u r f r = − . (2.47) Por su parte, la energía potencial total del par de iones es PIMs RIM ind(s) 12 ( ) ( ) ( ) r r r υ υ υ = + , (2.48) donde 6 q 2 2 2 2 2 ind(s) µ 12 1 2 2 1 1 2 1 2 12 6µ2 4 7 1 2 12 [ ( )] ( ) 21 ( ) ( ) ( 4 [ ( )] ) 2 r u r z z e z z e r u r r u r r r α α α α υα α   + = − +   −   (2.49) es la energía debida a la inducción de polarización en un PIMs. En el caso de un PIM2sq ( µ 12 ( ) 0 f r = y µ 12 ( ) 1 u r = ) los momentos dipolares de los dos iones y su energía de interacción son los mismos que en un PIM2e pero multiplicados por un factor q 12 ( ) u r o q 2 12 [ ( )] u r respectivamente. Sigue existiendo por tanto una singularidad en r C = (4 α 1 α 2 ) 1/6 , aunque los valores de los momentos y de la energía son menores que en un PIM2e. En cambio, esta singularidad puede evitarse utilizando un PIM2sqµ ( qµ 12 12 ( ) ( ) 0 f r f r = ≠ ). Como se desprende de los resultados (2.45)-(2.49), la 80 Capítulo 2. Modelos de ión polarizable [27] V. Bitrián, J. Trullàs, M. Silbert, T. Enosaki, Y. Kawakita and S. Takeda J. Chem. Phys. 125, 184510 (2006). [28] B. T. Thole Chem. Phys. 59, 341 (1981). [29] P. A. Madden and M. Wilson Chem. Soc. Rev. 25, 339 (1996). [30] Z. Akdeniz, G. Pastore and M. P. Tosi Nuovo Cimento 20, 595 (1998). [31] Z. Akdeniz and M. P. Tosi Z. Naturforsch. 54, 180 (1999). [32] R. Ruberto, G. Pastore and M. P. Tosi Phys. Lett. A 373, 1083 (2009). [33] R. Ruberto, G. Pastore and M. P. Tosi Physica B 405, 974 (2010). [34] P. W. Fowler and P. A. Madden Phys. Rev. B 31, 5443 (1985). [35] M. Masia, M. Probst and R. Rey J. Chem. Phys. 123, 164505 (2005). [36] C. J. Burnham, L. Jichen, S. S. Xantheas and M. Leslie J. Chem. Phys. 110, 4566 (1999). [37] K. T. Tang and J. P. Toennies J. Chem. Phys. 80, 3726 (1984). [38] P. A. Madden and M. Wilson J. Phys.: Condens. Matter 12, A95 (2000). [39] F. Hutchinson, M. Wilson and P. A. Madden Mol. Phys. 99, 811 (2001). [40] A. Aguado, M. Wilson and P. A. Madden J. Chem. Phys. 115, 8603 (2001). [41] A. Aguado, L. Bernasconi, S. Jahn and P. A. Madden Faraday Discuss. 124, 171 (2003). [42] B. Morgan and P. A. Madden J. Chem. Phys. 120, 1402 (2004). [43] N. Ohtori, M. Salanne and P. A. Madden J. Chem. Phys. 130, 104507 (2009). [44] V. Sarou-Kanian, A. L. Rollet, M. Salanne, C. Simon, C. Bessada and P. A. Madden Phys. Chem. Chem. Phys. 11, 11501 (2009). [45] V. Bitrián and J. Trullàs J. Phys. Chem. B 112, 1718 (2008). Capítulo 3 Modelos de ión polarizable para AgBr fundido. Efectos de la polarización inducida En este capítulo presentamos los resultados de simulaciones MD de AgBr fundido utilizando dos modelos en los que los aniones son polarizables. Ambos están basados en un modelo de ión rígido con potencial de tipo VR (al que denominaremos simplemente RIM). En uno de ellos (PIM1e) sólo el campo eléctrico induce la polarización en los aniones, mientras que en el otro (PIM1sq) también se tiene en cuenta la polarización de corto alcance. Este estudio ha sido publicado en Journal of Physical Chemistry B [1]. Su objetivo es determinar la influencia de la polarización inducida en las propiedades simuladas de AgBr fundido, y ver si el acuerdo con los resultados experimentales disponibles mejora al considerar modelos de ión polarizable. 82 Capítulo 3. Modelos de ión polarizable para AgBr fundido 3.1 Introducción Hace unos años Tasseven et al. [2], siguiendo el procedimiento establecido por Vashishta y Rahman en sus trabajos sobre las fases α de AgI [3] y CuI [4], parametrizaron potenciales VR para AgCl y AgBr, y calcularon los correspondientes factores de estructura estáticos coherentes S c (k) mediante la aplicación de la teoría HNC (Hypernetted Chain). Sin embargo los S c (k) que obtuvieron no reproducen el pico principal muy ancho y con un ligero rizado de tres picos que muestran los obtenidos por difracción de neutrones [5,6] (ver la Figura 1.9). Por la misma época, Wilson et al. [7] mostraron a partir de simulaciones MD que el triple rizado en el pico principal de S c (k) se pude reproducir para AgCl, aunque de forma exagerada, si se utilizan PIM en los que los efectos de los momentos dipolares y cuadripolares inducidos se añaden al potencial BMH propuesto por Mayer [8]. Este potencial de ión rígido predice una temperatura de fusión muy superior a la real, pero la adición de la polarización hace que esta temperatura disminuya. Sin embargo, los coeficientes de autodifusión que se obtienen son demasiado pequeños, y los valores de la conductividad iónica estimados a partir de la aproximación de Nernst-Einstein son significativamente menores que los experimentales. Posteriormente, Trullàs et al. [9] mostraron que el triple rizado en el pico del factor de estructura de AgCl también se reproduce en simulaciones MD si los efectos de los momentos dipolares inducidos en los aniones son añadidos al potencial VR parametrizado en [2], aunque las simulaciones exageran este rasgo al igual que ocurre en [7]. Además, obtuvieron unos resultados para la conductividad bastante acordes con los valores experimentales. El programa de simulación utilizado por el grupo de Trullàs en aquellas fechas sólo tenía en cuenta la polarización inducida por el campo eléctrico, es decir sólo permitía simular PIMe (ver Capítulo 2, en particular la Tabla 2.1), y por tanto no se podían estudiar sales en las que los cationes también fueran polarizables porque, como hemos visto en el apartado 2.4, en los PIM2e se produce la polarización catastrófica. En ese momento fue cuando se inició esta tesis, con el primer objetivo de ampliar el programa para poder simular también PIMs. De hecho, siguiendo las ideas de Wilson et al., la primera versión ampliada únicamente contemplaba la posibilidad de simular Capítulo 3. Modelos de ión polarizable para AgBr fundido 83 PIMsq, es decir aquellos modelos en los que la polarización de corto alcance amortigua la inducida por el campo eléctrico de las cargas puntuales, pero no la inducida por el campo de los dipolos puntuales asociados a cada ión. Esta primera versión es la que se aplicó para realizar el estudio de AgBr fundido que presentamos en este capítulo. De igual manera que Wilson et al. [7] hicieron con AgCl, primero utilizamos el potencial BMH propuesto para AgBr por Mayer [8], pero comprobamos que, incluso tras añadir la polarización inducida de los iones de bromo, el sistema simulado no fundía en el estado termodinámico para el cual hay datos experimentales disponibles. Por esta razón los modelos de ión polarizable que estudiamos están basados en un potencial VR con el que la sal sí funde a la temperatura y densidad adecuadas. En dos de ellos se incluye sólo la polarización de los iones de bromo; en uno (PIM1e) es inducida únicamente por el campo eléctrico, mientras que en el otro (PIM1sq) también se tiene en cuenta la polarización de corto alcance. En un tercero (PIM2sq) los cationes también son polarizables tanto por el campo eléctrico como por las deformaciones de corto alcance. Sin embargo nos encontramos con que, tal como hemos comentado en el apartado 2.6, tras tandas de simulaciones bastante largas en PIM2sq se producía de forma inesperada la polarización catastrófica. Por esta razón los resultados que presentamos se limitan a PIM1e y PIM1sq, que comparamos con los obtenidos con el modelo de ión rígido del tipo VR (al que denominamos RIM en este capítulo). 3.2 Modelos de AgBr simulados 3.2.1 Modelo de ión rígido Los modelos de AgBr que estudiamos se basan en el potencial efectivo a parejas que denominamos AgBr-VR7. Este potencial es del tipo VR y su forma funcional es VR VR0 4 ( ) ( ) ab ab ab P r r r υ υ = − (3.1) con 2 VR0 6 ( ) ( ) n a b a b ab ab n z z e A C r r r r σ σ υ + = + − , (3.2) donde n es el mismo para las tres interacciones, y 84 Capítulo 3. Modelos de ión polarizable para AgBr fundido 2 2 1 2 ( ) ab a b b a P z z α α = + y 3 2 /( ) ab a a a b a b C E E E E α α = + . (3.3) El significado físico de cada término está descrito en el apartado 1.2.2. Los valores de los parámetros del potencial AgBr-VR7 son en parte diferentes a los del propuesto por Tasseven et al [2], al que llamamos AgBrVR6 (ver Apendice C). Mientras que en AgBr-VR6 n = 6, en AgBr-VR7 impusimos n = 7. La razón de este cambio es que, en un trabajo sobre la mezcla de AgBr y AgI [10], para obtener de manera sencilla potenciales de interacción a partir de los de las sales puras, nos convino que n = 7 para AgBr como en la formulación original de los potenciales VR para AgI [3,4]. Para parametrizar el potencial AgBr-VR7 seguimos el procedimiento de Vashishta y Rahman [3,4], utilizado también por el grupo de Silbert en [11] y [2]. Igual que en [2] asumimos z Ag = -z Br = 0.66 [12], α Ag = 0 (y por tanto P AgAg = 0 y C AgAg = C AgBr = 0), α Br = 4.16 Å 3 [13], y E Br = 8.67 eV [8]. Los radios iónicos σ Ag = 0.846 Å y σ Br = 2.042 Å fueron estimados a partir de las distancias entre iones vecinos en la red cristalina de sal gema, suponiendo que dos iones de bromo contiguos están en contacto, y que lo mismo ocurre con uno de bromo y otro de plata contiguos. Es decir, 2 σ Br = a 0 /(2 1/2 ) y σ Ag + σ Br = a 0 /2 donde a 0 = 5.7745 Å es la constante de red a temperatura ambiente [14]. Con todos los valores anteriores fijados sólo faltaba determinar el valor de A. Primero estimamos su valor imponiendo la minimización de la energía de la red cristalina cuando la constante de red toma su valor experimental a 0 . Las aproximaciones termodinámicas asumidas y las sumas sobre las redes de Bravais involucradas en este cálculo se encuentran en [15]. Sin embargo las simulaciones MD revelaron que el sistema no fundía a la temperatura y densidad experimentales. Por ello fuimos reduciendo el valor de A hasta A = 0.1289 eV, para el cual el sistema funde en el estado termodinámico adecuado. En este capítulo al modelo de ión rígido con el potencial AgBr-VR7 lo denominamos simplemente RIM. A pesar de las diferencias introducidas, los resultados obtenidos con él no difieren de manera significativa de los obtenidos con AgBr-VR6. Capítulo 3. Modelos de ión polarizable para AgBr fundido 85 3.2.2 Modelos de ión polarizable Como se argumenta en el apartado 2.2, los modelos de ión polarizable basados en potenciales VR se construyen añadiendo las interacciones debidas a la polarización inducida a las del potencial de ión rígido VR0 de la ecuación (3.2). Los valores de los parámetros de este potencial son los mismos que los del AgBr-VR7 (excepto P ab = 0). En la Tabla 3.1 se resumen las características de los tres modelos de ión polarizable con los que trabajamos. En PIM1e y PIM1sq, α Ag = 0 y α Br = 4.16 Å 3 [13], que son los mismos valores utilizados en la parametrización de AgBr-VR7. Pero mientras el primero es un PIMe, el segundo es un PIMsq en el que la longitud característica de las deformaciones de corto alcance que amortiguan la polarización inducida por el campo eléctrico de las cargas iónicas es igual a λ q = 2.33 Å. Este valor es ligeramente mayor que el propuesto para AgCl por Wilson et al. [7], ya que σ Br > σ Cl . En el caso del tercer modelo, PIM2sq, ambas especies son polarizables con α Ag = 1.67 Å 3 [12] y α Br = 4.16 Å 3 , y el valor de λ q es el mismo que en PIM1sq. En ninguno de los tres modelos se tiene en cuenta la polarización de corto alcance que amortigua la inducida por el campo de los momentos dipolares, es decir λ µ = 0. Tabla 3.1: Polarizabilidades y longitud característica de la polarizacion de corto alcance en los PIM para AgBr. α Ag α Br λ q λ µ PIM1e 0 4.16 Å 3 0 0 PIM1sq 0 4.16 Å 3 2.33 Å 0 PIM2sq 1.67 Å 3 4.16 Å 3 2.33 Å 0 3.3 Detalles de las simulaciones Los distintos modelos de AgBr fueron simulados a 753 K y 883 K, con densidades iónicas ρ N = 0.0354 Å -3 y ρ N = 0.0346 Å -3 respectivamente [16]. 86 Capítulo 3. Modelos de ión polarizable para AgBr fundido Estas son las temperaturas a la cuales Inui et al. [6] realizaron sus experimentos de difracción de neutrones. La dependencia de los resultados con la temperatura es la esperada. A temperaturas más altas los picos de las funciones de distribución radial y de los factores de estructura parciales son más bajos y ligeramente más anchos, y los coeficientes de autodifusión y la conductividad iónica son mayores. Dado que estas diferencias no aportan demasiado a la discusión de las propiedades de la sal, y puesto que disponemos también de los datos experimentales del factor de estructura de Keen et al. a 703 K [17] y de Kawakita et al. a 773 K [18], sólo mostramos los resultados obtenidos a 753 K. Las simulaciones MD las llevamos a cabo considerando N = 1000 iones (500 de cada especie) en una celda cúbica de lado L con condiciones periódicas de contorno. Utilizamos el algoritmo de Beeman para integrar las ecuaciones del movimiento con un paso de tiempo de 5×10 -3 ps, y cada modelo se simuló en el equilibrio durante 10 6 pasos (5 ns). Los campos eléctricos y las energías y fuerzas electrostáticas se calcularon mediante el método de las sumas de Ewald. Todos los detalles acerca de las sumas de Ewald pueden encontrarse en [9]. Para obtener los momentos dipolares de los iones en cada configuración usamos el método iterativo de predicción y corrección, descrito en el apartado 2.3.2, con una precisión ∆ ∼ 10 -4 . La convergencia requiere unas 7 u 8 iteraciones en las simulaciones con PIM1e y PIM1sq, número que puede ascender a más de 40 en el caso de PIM2sq. Para simular un líquido cerca del punto de fusión es importante disponer de una configuración en la que los iones exhiban un comportamiento líquido a partir de la cual iniciar el proceso de equilibrio a la temperatura y densidad deseadas. Para conseguirla partimos de una estructura de sal gema y la simulamos a una temperatura mayor y una densidad menor para asegurarnos de que el sistema funde. Después, simulamos a una temperatura cada vez menor y una densidad cada vez mayor, hasta llegar a los valores deseados. Para asegurar que el sistema es realmente líquido y no se solidifica, como observamos en el caso de modelos basados en el potencial BMH, a cada temperatura es necesario realizar tandas de simulación sucesivas del orden de 10 4 pasos, y comprobar que los coeficientes de autodifusión no disminuyen de una a otra. Capítulo 3. Modelos de ión polarizable para AgBr fundido 87 Para comprobar que el sistema funde en todos los modelos, realizamos simulaciones preliminares con N = 216 iones, mucho más rápidas que con N = 1000. A pesar de que las funciones de distribución radial obtenidas con N = 216 son prácticamente idénticas a las obtenidas con N = 1000, es necesario utilizar un millar de partículas para calcular mejor los factores de estructura y los coeficientes de transporte. Polarización catastrófica de PIM2sq En el caso de PIM2sq, a pesar de que las simulaciones preliminares de 300×10 3 pasos (1.5 ns) con N = 216 iones se mantenían estables, en aquellas con N = 1000 se producía la polarización catastrófica al cabo de unos 20×10 3 pasos (0.1 ns). Más tarde comprobamos que esta también se observaba en simulaciones más largas con N = 216. Un análisis detallado de las configuraciones previas a la polarización catastrófica muestra que inicialmente hay un par de iones de signo opuesto, cuya separación es del orden de la típica entre primeros vecinos, con momentos dipolares bastante elevados y prácticamente paralelos. A partir de un determinado instante comienza un proceso de retroalimentación entre estos dos momentos dipolares y los de los iones de su entorno. Los primeros tienden a orientar a los segundos en su mismo sentido, y los segundos hacen crecer a su vez a los primeros. Este proceso provoca que, en cierto modo, la barrera que se observa en el caso de un catión y un anión aislados (ver la Figura 2.7) se reduzca, y que acabe produciéndose la polarización catastrófica. En este proceso llega un punto en que el método iterativo no converge y hay que aplicar el método exacto de inversión de la matriz descrito en 2.3.1. Las simulaciones MD de AgCl fundido llevadas a término por Wilson et al. [7] usando un modelo equivalente a nuestro PIM2sq, pero basado en el potencial BMH, no detectaron la polarización catastrófica. Este hecho puede ser debido a una de las siguientes tres razones: (a) sólo usaron 216 iones y sus simulaciones no fueron lo suficientemente largas para detectar una configuración crítica; (b) el método del Lagrangiano extendido que aplican, descrito en 2.3.3, evita las configuraciones críticas; o (c) no hay configuraciones críticas en el PIM2sq que utilizaron para AgCl fundido. En nuestra opinión, la opción (a) parece la más plausible, seguida por la (b). 88 Capítulo 3. Modelos de ión polarizable para AgBr fundido 3.4 Estructura de AgBr fundido A continuación analizamos las propiedades estructurales calculadas para RIM, PIM1e y PIM1sq. Como era previsible, PIM1sq exhibe rasgos estructurales intermedios entre los de RIM y PIM1e, pues la polarización de corto alcance cancela en buena medida los efectos de la polarización inducida por el campo eléctrico. 3.4.1 Funciones de distribución radial En la Figura 3.1 se representan las funciones de distribución radial g ab (r) calculadas a 753 K. En los tres modelos se aprecia la estructura de capas alternadas de signo opuesto alrededor de un ión. Sin embargo, la ordenación de carga alrededor de un catión es menos marcada en PIM1e. Su g AgAg (r), cuyos picos son más bajos y sus valles menos profundos que en RIM y PIM1sq, no está en fase con g BrBr (r) ni en oposición de fase con g AgBr (r). Su primer pico, que se encuentra en 3.55 Å, está desplazado a menores distancias respecto a los de RIM y PIM1sq, ambos en 3.95 Å. Además, la g AgAg (r) de PIM1e tiene un segundo máximo en 5.50 Å, entre el primero y el segundo de g AgBr (r), es decir que alrededor de cada catión hay una doble capa de cationes entre la primera y la segunda capa de aniones. Y el tercer máximo está en 8.25 Å, más a la derecha que el segundo máximo de g BrBr (r) y el segundo mínimo de g AgBr (r). Por su parte, la estructura anión-anión mostrada por g BrBr (r) es prácticamente la misma en los tres modelos. Esta similitud se puede explicar como resultado del empaquetamiento de los aniones. Estos, debido a su tamaño, mayor que el de los cationes, disponen de poco espacio para situarse y sólo pueden distribuirse de una determinada manera. Así pues, pese a que en PIM1e sólo tenemos en cuenta la polarizabilidad de los aniones, es la estructura catión-catión la que se ve más afectada por las interacciones debidas a la polarización inducida. Los rasgos enumerados de la g AgAg (r) para PIM1e indican un apantallamiento de la repulsión electrostática entre cationes, cuyo origen analizamos con más detalle al final de este subapartado. El hecho de que un ión de plata en PIM1e pueda acercarse más a un catión vecino que en RIM implica que, a su vez, se está alejando de otros cationes de su primera capa. Es por esto que el primer pico de g AgAg (r) en RIM se descompone en dos máximos en el caso de PIM1e. El Capítulo 3. Modelos de ión polarizable para AgBr fundido 89 de la izquierda (en 3.55 Å) corresponde a las parejas de cationes que se acercan, y el de la derecha (en 5.50 Å) a las parejas que se alejan. Figura 3.1: Funciones de distribución radial de AgBr fundido a 753 K. Izquierda: g AgAg (r) (curvas continuas), g AgBr (r) (curvas de puntos) y g BrBr (r) (curvas de trazos), para los tres modelos. Derecha: comparación de estas mismas funciones para RIM (curvas de puntos), PIM1sq (curvas de trazos) y PIM1e (curvas continuas). 1 2 3 4 g ab ( r ) RIM 1 2 3 g ab ( r ) PIM1e 0 1 2 3 g ab ( r ) 0 2 4 6 8 10 r (Å) PIM1sq 0 2 4 6 8 10 r (Å) Ag-Ag Ag-Br Br-Br RIM PIM1sq PIM1e 1 2 3 4 g AgAg ( r ) 1 2 3 g AgBr ( r ) 0 1 2 3 g BrBr ( r ) Ag-Ag Ag-Br Br-Br 96 Capítulo 3. Modelos de ión polarizable para AgBr fundido En la Figura 3.4 se representa S c (k) junto con los parciales, lo que permite relacionar los rasgos de uno con los de los otros. Por ejemplo, se ve con claridad que el máximo central de S c (k) para PIM1e está asociado al segundo máximo de S AgAg (k). Se advierte también que en la escala de los factores de estructura parciales los rasgos de S c (k) apenas se aprecian. Es necesario aumentar la escala como se hace en la Figura 3.5 para realzar sus características. 3.4.4 Comparación con propiedades estructurales experimentales De los tres modelos simulados, el S c (k) que concuerda mejor con los datos experimentales es el de PIM1e, a pesar de que el rizado del pico principal, muy suave en los datos experimentales, está exagerado en los resultados de la simulación. En la Figura 3.6 se compara el S c (k) para PIM1e con los obtenidos a partir de experimentos de difracción de neutrones por Inui et al. a 753 K [6], Keen et al. a 703 K [17] y Kawakita et al. [18] a 773 K. El ancho pico principal de los resultados de Inui es más alto que el de los otros dos, siendo el de Kawakita el más bajo. Estas diferencias son probablemente debidas a las diferentes correcciones aplicadas en el análisis, principalmente en lo que respecta al multiple scattering. Mientras Inui et al. lo corrigieron haciendo uso de datos termodinámicos, Kawakita et al. utilizaron el método de Blech y Averbach [21]. (De hecho, poco después de ser publicados los resultados de haluros de plata y cobre del grupo de Inui [6,22], Shirakawa et al. [23] repitieron el análisis de los mismos datos modificando la corrección de multiple scattering y obtuvieron valores más bajos del pico principal). También se observan diferencias a k pequeñas, que podrían ser debidas a que la S c (k) de Inui, que en esta región es la más baja, fue renormalizada para que su valor de S c (0) coincidiera con el calculado a partir de datos termodinámicos. A pesar de las discrepancias, los tres S c (k) experimentales muestran un shoulder en k ~ k M , tras el cual aumentan lentamente, exhibiendo un shoulder o un máximo muy suave aproximadamente en k ~ 2.5 Å -1 , hasta llegar al máximo absoluto en k ~ 3 Å -1 . Estos tres rasgos son los que dan lugar al ligero rizado citado con anterioridad, menos evidente en este caso que en el de AgCl fundido (ver Figura 1.9). En los datos experimentales de Keen se observa un shoulder a k’s pequeñas, parecido al que exhibe en k ~ 1 Å -1 el S c (k) simulado con PIM1e y que tiene su origen en el shoulder de S AgAg (k). A falta de datos a k’s menores, no sabemos si los altos valores que toma el S c (k) de Kawakita en Capítulo 3. Modelos de ión polarizable para AgBr fundido 97 k → 0 son parte de un shoulder, si se trata de un error asociado a la técnica utilizada (difractometría de tiempo de vuelo, que proporciona resultados poco fiables en esta región del dominio de las k), o si es consecuencia de no haber corregido el factor de estructura a k pequeñas [10]. Figura 3.6: Factores de estructura S c (k) de AgBr fundido obtenidos a partir de experimentos de difracción de neutrones por Inui et al. [6] a 753 K (círculos blancos), Keen et al. [17] a 703 K (círculos negros) y Kawakita et al. [18] a 773 K (círculos rojos). La curva continua es el S c (k) obtenido en las simulaciones MD con PIM1e a 753 K. Keen et al. [17], a partir del S c (k) extraído de la difracción de neutrones, realizaron simulaciones RMC (Reverse Monte Carlo) para obtener las tres funciones de distribución radial g ab (r). Por otro lado, Di Cicco et al. [24] determinaron el primer pico de g AgBr (r) para AgBr fundido a 725 K tras analizar los datos de sus experimentos de espectroscopía de absorción de rayos X (EXAFS). En la Figura 3.7 comparamos los resultados del primer pico de g AgBr (r) obtenido con estas dos técnicas con el de nuestras simulaciones MD con PIM1e. El pico más alto y estrecho corresponde a los datos EXAFS, que además muestran un mayor volumen excluido. La posición del primer pico en PIM1e es 2.60 Å, cercana a la de los datos RMC, 2.63 Å, y a la de los EXAFS, 2.67 Å, pero su altura es menor. El 0 1 2 3 4 5 6 k (Å -1 ) 0.0 0.4 0.8 1.2 S c ( k ) NS Inui [6] NS Keen [17] NS Kawakita [18] MD (PIM1e) 98 Capítulo 3. Modelos de ión polarizable para AgBr fundido hecho de que las g AgBr (r) experimentales presenten un mayor volumen excluido y sean más altas parece sugerir la conveniencia de utilizar una pared más repulsiva en la interacción Ag-Br del modelo. Figura 3.7: Funciones g AgBr (r) de AgBr fundido obtenidas mediante simulaciones RMC a 703 K a partir de datos de difracción de neutrones [17] (círculos negros) y mediante el análisis de datos EXAFS a 725 K [21] (círculos blancos). La curva continua corresponde a nuestras simulaciones MD con PIM1e a 753 K. 3.5 Transporte iónico en AgBr fundido 3.5.1 Autodifusión En la Tabla 3.2 se incluyen los coeficientes de autodifusión D a obtenidos para los distintos modelos de AgBr simulados a 753 K. En todos los casos los cationes, más pequeños que los aniones, se difunden más, con D Ag unas dos (o más) veces mayor que D Br , lo que recuerda en cierto modo el comportamiento en una fase superiónica. El coeficiente de autodifusión de los cationes D Ag es mayor en PIM1e que en RIM y PIM1sq, siendo el de este último el menor de los tres. La mayor difusividad de los iones de plata en PIM1e es consistente con el 1.6 2.0 2.4 2.8 3.2 r (Å) 0 1 2 3 4 g AgBr ( r ) RMC [17] EXAFS [21] MD (PIM1e) Capítulo 3. Modelos de ión polarizable para AgBr fundido 99 hecho de que su g AgAg (r) es la menos estructurada, y se explica por la fuerza que el momento dipolar de un anión ejerce sobre un catión. Como se ve en la Figura 3.2, esta fuerza no es radial sino que se dirige hacia el extremo negativo del dipolo, y permite que los cationes eviten en cierta medida las colisiones con los aniones vecinos y se difundan con mayor facilidad entre sus intersticios. Esta fuerza apantalla además, como hemos visto, la repulsión entre cationes. Es decir, de alguna manera los dipolos inducidos en los aniones aumentan el espacio libre disponible para los cationes, que se difunden a través de canales más rápidos que en RIM. Este efecto no se manifiesta en PIM1sq debido al amortiguamiento de la polarización. D Ag es menor en PIM1sq que en RIM porque la pared repulsiva de VR0 AgBr ( ) r υ , al no incluir el término –P AgBr /r 4 , está desplazada a r mayores que la de VR AgBr ( ) r υ , de forma que el espacio libre para los cationes es menor en PIM1sq que en RIM. Por lo que respecta a los aniones, su coeficiente D Br es menor en los PIM que en RIM, lo que también está relacionado con el hecho de que el potencial VR0 de los PIM comporta que el tamaño efectivo de los aniones sea mayor que en RIM. La polarización inducida contribuye al aumento del D Br en PIM1e respecto al de PIM1sq, pero no lo suficiente como para superar al de RIM. Tabla 3.2: Coeficientes de autodifusión D a y conductividades iónicas σ calculadas en las simulaciones MD de AgBr a 753 K. RIM PIM1sq PIM1e D Ag (10 -5 cm 2 /s) 3.85 2.34 5.00 D Br (10 -5 cm 2 /s) 2.00 1.09 1.42 σ (Ω cm) -1 3.3 a 1.9 a 2.8 a σ NE (Ω cm) -1 2.56 1.50 2.81 ∆ -0.30 -0.28 0 > ɶ a Valor experimental de σ = 2.97 (Ω cm) -1 [16] 100 Capítulo 3. Modelos de ión polarizable para AgBr fundido Figura 3.8: Funciones de autocorrelación de velocidades C Ag (t) (panel superior) y C Br (t) (panel central), y función de autocorrelación de la corriente de carga C ZZ (t) (panel inferior), de AgBr a 753 K para RIM (curvas de puntos), PIM1sq (curvas de trazos) y PIM1e (curvas continuas). En la Figura 3.8 se comparan las funciones de autocorrelación de velocidades normalizadas C a (t) = Λ a (t)/Λ a (0). Su comportamiento es el característico de los haluros de plata y cobre expuesto en el apartado 1.4.1. En los tres modelos C Br (t) oscila tras el decaimiento inicial, ya que los aniones, mayores que los cationes, experimentan un movimiento de vibración en la doble caja formada temporalmente por los primeros vecinos 0.0 0.2 0.4 0.6 t ( p s ) -0.4 0.0 0.4 0.8 C ZZ ( t ) 0.0 0.4 0.8 C Br ( t ) 0.0 0.4 0.8 C Ag ( t ) RIM PIM1sq PIM1e Capítulo 3. Modelos de ión polarizable para AgBr fundido 101 de signo opuesto y los primeros del mismo signo [25]. En cambio los cationes se difunden a través de la estructura altamente empaquetada de los aniones, de modo que C Ag (t) presenta un decaimiento más lento y un mínimo de backscattering bastante menos profundo. Este contraste entre C Ag (t) y C Br (t) se ve reforzado en este caso porque los iones de bromo son, además de más grandes, más ligeros que los de plata, lo que favorece el comportamiento vibratorio de los primeros. El decaimiento inicial de C Ag (t) es prácticamente igual para RIM y PIM1sq, y sus mínimos de backscattering se encuentran aproximadamente en el mismo tiempo, aunque el segundo es más profundo. En cambio, la C Ag (t) para PIM1e exhibe la caída inicial más lenta y el backscattering menos profundo y a tiempos mayores. Estas características concuerdan con el aumento de la difusividad que comporta la inclusión de la polarización inducida y que ya hemos discutido antes. Por otro lado, las tres C Br (t) oscilan con prácticamente la misma frecuencia, y con un backscattering más pronunciado cuanto mayor es el valor de D Br . 3.5.2 Conductividad iónica Los valores obtenidos para la conductividad iónica σ están recogidos en la Tabla 3.2. Como en trabajos anteriores [9,2], y como hemos comentado en 1.4.2, en las ecuaciones (1.37)-(1.41) para el cálculo de σ hemos tomado z = 1 en las funciones C ZZ (t) y Q ZZ (t). El valor de la conductividad obtenido con PIM1e, 2.8 (Ω cm) -1 , está en buen acuerdo con el determinado experimentalmente [16], 2.97 (Ω cm) -1 . Por el contrario, el obtenido con RIM, 3.3 (Ω cm) -1 , es mayor, y el de PIM1sq, 1.9 (Ω cm) -1 , menor. Estas relaciones pueden explicarse en términos de la forma de las funciones C ZZ (t) = Λ ZZ (t)/Λ ZZ (0) en la Figura 3.8. Cuanto más profundo es su mínimo, menor es la conductividad iónica. En este punto hay que recalcar que conviene interpretar con prudencia la concordancia entre la conductividad iónica experimental y la obtenida en las simulaciones, ya que esta última depende fuertemente del potencial de ión rígido utilizado. Así, pequeños cambios en el potencial pueden propiciar, independientemente de que la polarización sea tenida en cuenta o no, cambios significativos en la conductividad calculada. 102 Capítulo 3. Modelos de ión polarizable para AgBr fundido Atendiendo a la Tabla 3.2 vemos que σ es mayor que la aproximación de Nernst-Einstein, σ NE = ρ N e 2 (D + + D - )/2k B T , para RIM y para PIM1sq, y que son prácticamente iguales para PIM1e. Por tanto el parámetro ∆ que indica el grado de discrepancia entre σ NE y σ = σ NE (1-∆), relacionado con las correlaciones entre velocidades de iones distintos, es ∆ < 0 para RIM y PIM1sq y ∆0 > ɶ para PIM1e. Estos resultados están en la misma línea que los de simulaciones anteriores de haluros de plata y cobre fundidos usando RIM [2,26], en las que se obtuvieron sistemáticamente valores negativos de ∆, y de AgCl usando un PIMe, en las que ∆ > 0 [9]. Desafortunadamente no se dispone de datos experimentales sobre los coeficientes de autodifusión de estas sales, con la excepción del CuCl fundido (para el cual ∆ > 0) [27], con lo que los valores MD de ∆ no pueden ser comparados con los experimentales. Figura 3.9: Función δ ZZ (t) de AgBr a 753 K para RIM (curva de puntos), PIM1sq (curva de trazos) y PIM1e (curva continua). En la Figura 3.9 se representa la función de correlación temporal δ ZZ (t) = {Λ ZZ (t)-1/2[Λ + (t) + Λ - (t)]}/Λ ZZ (0) definida en el apartado 1.4.2, que da cuenta de las correlaciones entre velocidades de iones distintos y cuya integral cambiada de signo es proporcional al parámetro ∆. Se advierte que en todos los casos δ ZZ (t) comienza tomando valores positivos. Este comportamiento ha sido observado también en simulaciones MD de haluros -0.2 0.0 0.2 δ ZZ ( t ) 0.0 0.2 0.4 0.6 0.8 t (ps) RIM PIM1sq PIM1e Capítulo 3. Modelos de ión polarizable para AgBr fundido 103 de cobre fundidos usando RIM [28], así como para AgCl fundido utilizando tanto RIM como PIM [9]. A pesar de que la forma de δ ZZ (t) es similar en los tres modelos, las diferencias de amplitud de sus oscilaciones positivas y negativas hacen que sus integrales, y por tanto los valores de ∆, cambien. 3.6. Conclusiones Hemos visto que en un modelo PIM2sq, en el que tanto los aniones como los cationes son polarizables, el amortiguamiento de la polarización inducida por el campo eléctrico de las cargas no garantiza evitar la polarización catastrófica. Esta se produce en configuraciones críticas poco probables pero que acaban dándose tras tandas de simulación largas. Como veremos en el próximo capítulo, esta anomalía puede resolverse si se contemplan modelos PIM2sqµ, en los que también se amortigua la polarización inducida por el campo eléctrico creado por los dipolos. De los tres modelos de AgBr estudiados, PIM1e es el que mejor reproduce el factor de estructura S c (k) y la conductividad iónica experimentales. Este hecho indica que el papel jugado por la polarización inducida es importante en AgBr fundido y que no puede omitirse en su descripción. El efecto más destacado de la polarización inducida en los aniones es el apantallamiento de la repulsión electrostática entre dos cationes vecinos, que permite que estos últimos se acerquen más de lo que lo harían en RIM. Entonces, la primera capa de cationes que rodea a un catión, en PIM1e se desdobla en dos, de modo que g AgAg (r) deja de estar en fase con g BrBr (r) y en oposición de fase con g AgBr (r) como ocurre en RIM. En el espacio recíproco, la ordenación de carga menos marcada se pone de manifiesto en la disminución y ensanchamiento del pico principal del S AgAg (k) de PIM1e, que se traduce en un pico de S ZZ (k) menos pronunciado. S AgAg (k) exhibe además dos rasgos distintivos: la aparición de un shoulder en k ~ 1 Å -1 , cuyo origen analizaremos en el siguiente capítulo, y un segundo máximo en k ~ 2.3 Å -1 . Este último está relacionado con la distancia entre el segundo y el tercer pico de g AgAg (r), y es el responsable de que S c (k) muestre en PIM1e un ancho pico principal con un triple rizado, que también está presente en los datos experimentales, aunque los resultados MD lo exageran. 104 Capítulo 3. Modelos de ión polarizable para AgBr fundido Desde un punto de vista dinámico, la polarización inducida en los aniones deja, de algún modo, más espacio libre para que los cationes puedan difundirse mejor a través de la compacta estructura de los aniones (más grandes que los cationes). Referencias [1] V. Bitrián and J. Trullàs J. Phys. Chem. B 110, 7490 (2006). [2] Ç. Tasseven, J. Trullàs, O. Alcaraz, M. Silbert and A. Giró J. Chem. Phys. 106, 7286 (1997). [3] P. Vashishta and A. Rahman Phys. Rev. Lett. 40, 1337 (1978). [4] P. Vashishta and A. Rahman, in Fast Ion Transport in Solids (p. 527), edited by P. Vashishta, J. N. Mundy and G. K. Shenoy, North-Holland, Amsterdam (1979). [5] J. Y. Derrien and J. Dupuy Phys. Chem. Liq. 5, 71 (1976). [6] M. Inui, S. Takeda, Y. Shirakawa, S. Tamaki, Y. Waseda and Y. Yamaguchi J. Phys. Soc. Japan 60, 3025 (1991). [7] M. Wilson, P. A. Madden and B. J. Costa-Cabral J. Phys. Chem. 100, 1227 (1996). [8] J. E. Mayer J. Chem. Phys. 1, 327 (1933). [9] J. Trullàs, O. Alcaraz, L. E. González and M. Silbert J. Phys. Chem. B 107, 282 (2003). [10] V. Bitrián, J. Trullàs, M. Silbert, T. Enosaki, Y. Kawakita and S. Takeda J. Chem. Phys. 125, 184510 (2006). [11] A. J. Stafford, M. Silbert, J. Trullàs and A. Giró J. Phys.: Condens. Matter 2, 6631 (1990). Las Figuras 2 y 3 en esta referencia deben ser intercambiadas para coincidir con los pies de figura. Además, en la Tabla 1 el valor de α X para AgI debería ser 6.52 Å 3 y el de H XX para CuI debería ser 399 en unidades de e 2 Å 6 = 14.4 eVÅ 7 . [12] M. Bucher Phys. Rev. B 30, 947 (1984). [13] J. R. Tessman, A. H. Kahn and W. Shockley Phys. Rev. 92, 890 (1953). [14] R. W. G Wyckoff, Crystal Structures, 2nd ed., John Wiley & Sons, New York (1963). [15] M. P. Tosi, in Solid State Physics: Advances in Research and Application (vol. 16, p. 1), edited by F. Seitz, D. Turnbull and H. Ehrenreich, Academic, New York (1964). [16] G. J. Janz, F. W. Dampier, G. R. Lakshminarayanan, P. K. Lorenz and R. P. T. Tomkins, Molten Salts, National Bureau of Standards Reference Data Series Nº 15, Washington DC (1968). [17] D. A. Keen, W. Hayes and R. L. McGreevy J. Phys.: Condens. Matter 2, 2773 (1990). Capítulo 3. Modelos de ión polarizable para AgBr fundido 105 [18] Y. Kawakita, T. Enosaki, S. Takeda and K. Maruyama J. Non-Cryst. Solids 353, 3035 (2007). [19] V. M. Nield, D. A. Keen, W. Hayes and R. L. McGreevy J. Phys.: Condens. Matter 4, 6703 (1992). [20] V. F. Sears Neutron News 3, 26 (1992). [21] I. A. Blech and B. L. Averbach Phys. Rev. 137, A1113 (1965). [22] Y. Shirakawa, M. Saito, S. Tamaki, M. Inui and S. Takeda J. Phys. Soc. Japan 60, 2678 (1991). [23] Y. Shirakawa, S. Tamaki, T. Usuki, K. Sugiyama and Y. Waseda J. Phys. Soc. Japan 63, 1814 (1994). [24] A. Di Cicco, M. Taglienti, M. Minicucci and A. Filipponi Phys. Rev. B 62, 12001 (2000). [25] O. Alcaraz and J. Trullàs J. Chem. Phys. 113, 10635 (2000). [26] J. Trullàs, A. Giró and M. Silbert J. Phys. Condens. Matter 2, 6643 (1990). [27] J. C. Poignet and M. J. Barbier Electroch. Acta 26, 1429 (1981). [28] J. Trullàs and J. A. Padró Phys. Rev. B 55, 12210 (1997). 112 Capítulo 4. Estudio de α αα α-AgI y AgI fundido Beeman, durante 10 6 pasos de integración (5 ns) en el equilibrio. Los detalles del método de Ewald usado para calcular los términos de largo alcance son los mismos que los especificados en [29]. Las simulaciones de α-AgI las realizamos con N = 500 iones, a una temperatura T = 573 K y una densidad iónica ρ N = 0.0306 Å -3 [8]. Esta temperatura es la misma a la cual Cava et al. [7] llevaron a cabo sus experimentos de difracción de neutrones. En las simulaciones de esta fase tomamos configuraciones iniciales en las que los iones de iodo estaban situados en los nudos de una red bcc con un parámetro de red (lado de la celda unidad) a = 5.075 Å igual al experimental [8], y los de plata en dos posiciones tetraédricas por celda unidad de α-AgI (o lo que es lo mismo, en 4 de los 24 sitios tetraédricos que aparecen en la Figura 4.1, ya que cada uno de ellos pertenece a dos celdas unidad contiguas). Para asegurarnos de que la posición inicial de los iones de plata no es relevante, hicimos diversas pruebas con otras configuraciones iniciales en las que los iones de plata estaban situados en posiciones trigonales u octaédricas, y comprobamos que efectivamente los resultados obtenidos en el equilibrio no variaban. Las simulaciones de AgI fundido las llevamos a cabo a T = 923 K y ρ N = 0.0281 Å -3 [33], el mismo estado termodinámico que en los experimentos de difracción de neutrones de Kawakita et al. [24]. Para obtener las configuraciones iniciales líquidas seguimos el mismo procedimiento comentado para AgBr en el apartado 3.3. Los resultados de AgI fundido que mostramos en esta tesis fueron obtenidos con N = 1000 iones. Pese a que en [2] utilizamos N = 500, posteriormente hicimos nuevas simulaciones con un millar de partículas para disminuir el mínimo valor de k accesible en el cálculo de los factores de estructura. El tiempo computacional requerido por cada uno de los tres modelos es bastante distinto. En las simulaciones con N = 500 iones, las tandas con PIM1e y PIM2sqµ tardan 1.5 y 12 veces más, respectivamente, que aquellas con RIM. El paso de N = 500 a N = 1000 comporta un aumento notable del tiempo de cálculo, que se multiplica por ~ 3 en el caso de RIM y PIM1e y por ~ 6 en PIM2sqµ. Así, las simulaciones con PIM2sqµ pasan a ser unas 24 veces más costosas que las de RIM, y 16 veces más que las de PIM1e. El hecho de que el incremento del tiempo al pasar de N = 500 a N = 1000 sea mayor en el caso de PIM2sqµ tiene que ver con el cálculo de los momentos Capítulo 4. Estudio de α αα α-AgI y AgI fundido 113 dipolares. En este último modelo, el método iterativo de prediccióncorrección no converge en una pequeña fracción de los pasos, en los cuales es necesario recurrir al método exacto de inversión de la matriz, que es mucho más lento (ver el apartado 2.3). Esta fracción representa el 0.1 % de los pasos totales si N = 500, y aumenta hasta el 2 % cuando N = 1000. Además, en el resto de pasos en los que el método iterativo sí converge, las iteraciones necesarias para obtener los momentos con una precisión ∆ = 10 -6 pasan a ser de alrededor de 100 si N = 500 a oscilar entre 100 y 300 si N = 1000. Estos números contrastan con las ~ 15 iteraciones necesarias en PIM1e. Además de los dos PIM de la Tabla 4.1, cuyos resultados analizaremos y compararemos con los de RIM en los apartados siguientes, también simulamos un PIM2sq con las mismas polarizabilidades que en PIM2sqµ y con λ q = 1.86 Å (y λ µ = 0). Como ya comentamos en los apartados 2.6 y 3.3, la polarización catastrófica se produce tras varias decenas de miles de pasos, de manera inesperada pero inevitable. 4.4 Estructura de α αα α-AgI Las simulaciones de los tres modelos (RIM, PIM1e y PIM2sqµ) reproducen el comportamiento superiónico de la fase α. En los tres casos se observa que los iones I - no se difunden y permanecen vibrando térmicamente alrededor de sus posiciones iniciales en la red bcc, con un desplazamiento cuadrático medio de unos 0.2 Å 2 . Por el contrario, los iones de plata exhiben un comportamiento propio de un líquido, difundiéndose a través de la estructura que forman los aniones. Las propiedades de transporte iónico obtenidas para α-AgI las comentaremos en el apartado 4.7. 4.4.1 Distribución de los cationes en α αα α-AgI En las simulaciones de los tres modelos se observa que los iones de plata tienden a estar situados la mayor parte del tiempo en las caras de las celdas cúbicas de la red bcc aniónica. En la Figura 4.2 se muestra el mapa de la densidad media de los cationes en estas caras (más precisamente, en el interior de una lámina de grosor a/16 centrada en la cara). Los mapas de densidad resultantes son muy parecidos en los tres modelos, por lo que sólo mostramos uno de ellos. En él se aprecian regiones de mayor probabilidad 114 Capítulo 4. Estudio de α αα α-AgI y AgI fundido de ocupación alrededor de las cuatro posiciones tetraédricas (indicadas con una T) y, en menor medida, alrededor de las posiciones trigonales (X) y octaédricas (O). Figura 4.2: Izquierda: densidad media de cationes en una cara de la celda unidad de α-AgI, obtenida en las simulaciones de RIM a 573 K. Los iones de iodo (I) están situados en los cuatro vértices. La densidad se ha normalizado respecto a su valor máximo (amarillo), que se da en las posiciones tetraédricas (T), y las posiciones en los ejes se indican en unidades de a. Los sitios trigonales y octaédricos vienen indicados por X y O respectivamente. Derecha: perfiles de densidad en el segmento que une dos posiciones tetraédricas y una trigonal (T-X-T); dos tetraédricas y una octaédrica (T-O-T); y un anión, una trigonal y una octaédrica (I-X-O). Si nos fijamos en los perfiles de densidad a la derecha de la Figura 4.2, vemos que en los segmentos T-X-T y T-O-T la ocupación en los sitios trigonales y octaédricos es mínima, lo que implica que los iones de plata no 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 .2 0 .4 0 .6 0 .8 1 .0 T X 0 .2 0 .4 0 .6 0 .8 1 .0 0 .0 0.2 0.4 0 .6 0 .8 XXO I I T T O I-X-O T-O-T T-X-T T T T T T X X X X O I I I I Capítulo 4. Estudio de α αα α-AgI y AgI fundido 115 residen en ellos sino que los utilizan como lugares de paso entre los tetraédricos. Este resultado es consistente con los estudios experimentales, que establecen que el movimiento difusivo de los cationes tiene lugar mediante saltos entre posiciones tetraédricas adyacentes, principalmente a través de las posiciones trigonales y en menor medida a través de las octaédricas [7,8,4,5]. Las posiciones trigonales son puntos de silla, ya que su densidad es mínima en el segmento T-X-T, pero máxima en el I-X-O; en cambio, las posiciones octaédricas son mínimos locales. 4.4.2 Funciones de distribución radial en α αα α-AgI Ni g II (r) ni g AgI (r) presentan diferencias significativas de un modelo a otro; las primeras son prácticamente idénticas y en las segundas sólo varía ligeramente la altura del primer pico. A la izquierda de la Figura 4.3 se muestran las funciones de distribución radial g ab (r) para PIM1e. En cambio sí se observan diferencias importantes entre las g AgAg (r) presentadas a la derecha para los tres modelos. Figura 4.3: Izquierda: funciones de distribución radial de α-AgI a 573 K obtenidas con PIM1e: g AgAg (r) (curva continua), g AgI (r) (curva de puntos) y g II (r) (curva de trazos). Derecha: g AgAg (r) para RIM (curva de puntos), PIM1e (curva continua), y PIM2sqµ (curva de trazos). Adviértase la diferencia entre las escalas verticales de las dos gráficas. PIM1e 0 1 2 3 4 g ab ( r ) 0 2 4 6 8 10 r (Å) 0 2 4 6 8 10 r (Å) 0.0 0.4 0.8 1.2 1.6 g AgAg ( r ) Ag-Ag RIM PIM2sq µ PIM1e Ag-I Ag-Ag I-I 116 Capítulo 4. Estudio de α αα α-AgI y AgI fundido La forma de g II (r) es la propia de una red bcc a altas temperaturas, en la cual la primera capa de 8 aniones que rodean a un anión, y la siguiente de 6, se confunden y dan lugar al primer pico de g II (r) en a√3/2 = 4.4 Å, con un número de coordinación igual a 14. Por su parte, el primer pico de g AgI (r) está en 2.7 Å, muy cercano a la distancia entre un anión y un catión en una de las posiciones tetraédrica más cercanas (a√5/4 = 2.8 Å), y su número de coordinación es 4. Por lo que respecta a las g AgAg (r), la diferencia más notable es que en PIM1e el primer pico está en 3.2 Å, desplazado a la izquierda respecto a los de RIM y PIM2sqµ, que se encuentran en 4.3 Å, casi la misma posición del primer pico de g II (r). Como discutimos en el apartado 3.4.1, este desplazamiento es debido al apantallamiento de la repulsión entre los iones de plata como consecuencia de la inducción de polarización en los aniones, que les permite acercarse más. La segunda diferencia es que la g AgAg (r) de PIM1e tiene un shoulder (en realidad un máximo muy poco pronunciado) en ~ 6 Å. En el caso de PIM2sqµ, aunque la posición del primer pico es muy cercana a la de RIM a causa de los efectos amortiguadores de la polarización de corto alcance, también se observa un shoulder aproximadamente en la misma posición del shoulder de PIM1e. La forma de g AgAg (r) se puede relacionar con la ocupación de las posiciones tetraédricas indicadas en la Figura 4.4. En los tres modelos existe una probabilidad nula de que dos iones de plata estén separados una distancia inferior a 1.8 Å, lo cual excluye la ocupación simultánea de dos posiciones tetraédricas vecinas. La posición en 4.3 Å del primer pico de g AgAg (r) en RIM y PIM2sqµ es aproximadamente la separación entre dos posiciones tetraédricas unidas por el vector (1/2,1/2,1/2) a, a√3/2 = 4.4 Å (por ejemplo, los sitios 0 y 2 de la Figura 4.4). De entre todas las maneras de distribuir dos cationes por celda unidad en los sitios tetraédricos, esta es la máxima separación posible. En cambio, el primer pico en 3.2 Å para PIM1e implica que dos iones de plata vecinos tienden a localizarse alrededor de sitios tetraédricos unidos por (1/2,1/4,1/4) a, o vectores equivalentes del mismo módulo, separados una distancia a√6/4 = 3.1 Å (por ejemplo, los sitios 0 y 1 de la Figura 4.4). Capítulo 4. Estudio de α αα α-AgI y AgI fundido 117    Figura 4.4: La distancia entre los sitios tetraédricos 0 y 1 es la típica entre dos iones de plata vecinos en PIM1e. La distancia entre 0 y 2 es la típica en RIM y PIM2sqµ. El análisis RMC (Reverse Monte Carlo) realizado por Nield et al. [8] a partir de datos experimentales de difracción de neutrones, así como las simulaciones ab initio MD llevadas a cabo por Wood y Marzari [34] también presentan un primer pico de g AgAg (r) alrededor de 3 Å. Sin embargo, Nield et al. obtuvieron además un segundo pico en aproximadamente la misma posición del primer pico de g AgAg (r) en RIM y PIM2sqµ, ~ 4.3 Å. El hecho de obtener dos picos de g ++ (r), uno a la misma distancia que el primero de g -- (r) y otro a una r menor, es una constante en las simulaciones RMC de haluros de plata y cobre del grupo de McGreevy. Por ejemplo, este rasgo fue obtenido por estos autores en sus trabajos sobre las fases sólidas de CuBr [35] y CuI [36] y sobre AgBr [37]. En el artículo acerca de CuI, los autores mostraron además que el pico a r menor no se reproduce en simulaciones MD con un modelo de ión rígido [38], lo que está en consonancia con nuestros resultados para α-AgI e indica que su origen puede estar relacionado con los efectos de la polarización. 4.4.3 Distribución de los ángulos entre enlaces en α αα α-AgI Otra manera de caracterizar la estructura es mediante la distribución de probabilidad de los ángulos entre enlaces mostrada en la Figura 4.5. Adoptamos para ello la definición de enlace descrita en 3.4.2. El ángulo I-I-I es el formado por los enlaces entre un ión central de iodo y dos vecinos de su misma especie, y Ag-I-Ag es el que forman los enlaces entre un ión central de iodo y dos iones de plata vecinos. Como era de esperar, la distribución del ángulo I-I-I es casi idéntica para los tres modelos, ya que en 118 Capítulo 4. Estudio de α αα α-AgI y AgI fundido los tres corresponde a una red bcc a temperaturas altas. Por el contrario, la distribución de Ag-I-Ag es diferente en cada modelo. La de PIM1e es máxima en ~ 70º, cerca de 66º, que es el ángulo que forman dos posiciones tetraédricas separadas por el vector (1/2,1/4,1/4) a, o uno equivalente, y un ión de iodo primer vecino de las dos (por ejemplo, los sitios 0 y 1 de la Figura 4.4 y el anión central). En cambio, las de RIM y PIM2sqµ presentan un máximo en ~ 100º y ~ 96º respectivamente. Estos ángulos son parecidos al de 101º que forman dos sitios tetraédricos separados por (1/2,1/2,1/2) a y un anión primer vecino de ambos (como los sitios 0 y 2 y el anión central de la Figura 4.4). Por tanto, estas distribuciones angulares confirman lo concluido anteriormente a partir de las funciones de g AgAg (r). Figura 4.5: Distribuciones de probabilidad de los ángulos φ formados por los enlaces entre un ión de iodo y (izquierda) dos iones de iodo vecinos, I-I-I, o (derecha) dos iones de plata vecinos, Ag-I-Ag, obtenidas para α-AgI a 573 K con RIM (curvas de puntos), PIM1e (curvas continuas) y PIM2sqµ (curvas de trazos). φ (rad) 0 .0 0.5 1.0 1.5 I-I-I φ (rad) Ag-I-Ag 3 π/4 0 π/4 π/2 π/4 π/2 3 π/4 π Capítulo 4. Estudio de α αα α-AgI y AgI fundido 119 4.5 Estructura de AgI fundido 4.5.1 Funciones de distribución radial Las funciones de distribución radial g ab (r) obtenidas para AgI fundido se muestran en la Figura 4.6. La forma de g II (r) ya no es la correspondiente a una red cristalina como en α-AgI, sino la típica de las g ab (r) de un líquido. La altura de los picos y la profundidad de los mínimos de g AgAg (r) y g AgI (r) disminuyen apreciablemente respecto a los de la fase α, una consecuencia lógica del aumento de la temperatura y la disminución de la densidad. La forma de las g ab (r) de los tres modelos es cualitativamente parecida a la descrita para AgBr fundido en el apartado 3.4.1 (pero con I en lugar de Br). g II (r) es casi idéntica en los tres casos, y las únicas diferencias apreciables entre las g AgI (r) atañen a la altura del primer pico, más alto en RIM, y a su posición: 2.65 Å en RIM, 2.70 Å en PIM1e y 2.80 Å en PIM2sqµ. Este ligero corrimiento del primer pico de g AgI (r) en los PIM respecto a RIM se explica por el mayor tamaño efectivo de los iones en el potencial VR0 (sin el término -P AgI /r 4 ), compensado en parte en PIM1e por la polarización inducida. Por su parte, las g AgAg (r) presentan una estructura mucho menos marcada que las g II (r), y muestran una elevada penetración de los cationes en la primera capa de aniones que rodea a otro catión. Esta penetración es mayor que en AgBr, hasta el punto que en los dos PIM la cola a r pequeñas de g AgAg (r) interseca la de g AgI (r). Como ocurre en AgBr fundido y α-AgI, y como consecuencia del apantallamiento de la repulsión electrostática entre cationes debido a la polarización inducida en los aniones, la g AgAg (r) de PIM1e tiene su primer pico desplazado hacia la izquierda, en 3.45 Å, corrido casi 1 Å respecto al de RIM y al primer pico de g II (r), que se encuentran en 4.35 Å. Además, como se aprecia en el recuadro interior de la Figura 4.6, la g AgAg (r) de PIM1e exhibe un segundo máximo en 5.85 Å, entre el primer y segundo pico de g AgI (r), en una posición muy cercana a la del shoulder de g AgAg (r) en la fase α. En conjunto, la g AgAg (r) de PIM1e deja de estar en fase con g II (r) y en oposición de fase con g AgI (r). El comportamiento de g AgAg (r) en PIM2sqµ es intermedio al de los otros dos modelos, ya que su primer pico está en 4.10 Å, la intensidad de sus oscilaciones posteriores es mayor que en PIM1e pero menor que en RIM, y exhibe un shoulder a r ~ 5.5 Å, cerca del segundo máximo en PIM1e. 120 Capítulo 4. Estudio de α αα α-AgI y AgI fundido Figura 4.6: Funciones de distribución radial de AgI fundido a 923 K para RIM (panel superior), PIM2sqµ (panel central) y PIM1e (panel inferior): g AgAg (r) (curvas continuas), g AgI (r) (curvas de puntos) y g II (r) (curvas de trazos). En los recuadros interiores se muestran con mayor detalle las g AgAg (r) para PIM2sqµ y PIM1e. 0 2 4 6 8 10 12 14 1 6 r (Å) 0 1 2 3 g ab ( r ) 1 2 3 g ab ( r ) 1 2 3 4 g ab ( r ) RIM PIM1e PIM2sq µ 3 5 7 9 1.0 1.2 3 5 7 9 1.0 1.2 Ag-Ag Ag-I I-I Capítulo 4. Estudio de α αα α-AgI y AgI fundido 121 Si se comparan las g ab (r) para AgI y AgBr fundidos se observa que los picos en AgI están desplazados a r mayores que en AgBr, puesto que el tamaño de los iones de iodo es mayor que el de los de bromo (y la densidad de partículas es menor en AgI). La única excepción es la g AgAg (r) de PIM1e, cuyo primer pico en AgI se encuentra a la izquierda del de AgBr. Este hecho se explica probablemente porque α I > α Br y el apantallamiento de la repulsión entre cationes es más intenso en AgI. 4.5.2 Factores de estructura El comportamiento cualitativo de los factores de estructura parciales S ab (k) representados en la Figura 4.7 es parecido al descrito en el apartado 3.4.3 para AgBr fundido. El pico o valle principal en el número de onda k M ~ 2π/λ ~ 1.7 Å -1 está asociado a la distancia λ ~ 3.7 Å entre oscilaciones consecutivas de las g ab (r). Mientras que los S II (k) son casi idénticos en los tres modelos, y los S AgI (k) muy parecidos, los S AgAg (k) presentan diferencias destacables. El pico de S AgAg (k) en k M es mucho más bajo en PIM1e que en RIM, como consecuencia de la estructura menos marcada de los cationes en PIM1e y del hecho que las oscilaciones de longitud de onda igual a 3.7 Å en su g AgAg (r) sólo se detectan a partir del tercer máximo. Además, el S AgAg (k) de PIM1e presenta un segundo máximo en k = 2.25 Å -1 , más ancho y ligeramente más alto que el de k M , asociado a la separación entre el segundo y el tercer máximo de su g AgAg (r), en 5.85 y 8.65 Å respectivamente. También en PIM2sqµ se observa un segundo máximo tras el pico principal, más bajo que el de PIM1e y ligeramente desplazado hacia la izquierda. Pero sin duda el rasgo más llamativo de S AgAg (k) en PIM1e es el prepico que exhibe en k ~ 1 Å -1 , que no se observa ni en RIM ni en PIM2sqµ. Tampoco está presente en las simulaciones de AgCl [29] y AgBr fundidos utilizando PIM1e. En las de AgBr aparece en su lugar un shoulder, que sólo se intuye en las de AgCl. En este punto conviene precisar que la diferencia entre el prepico de S AgAg (k) en el PIM1e de AgI y el shoulder en el de AgBr es que el primero tiene un máximo relativo (apenas perceptible, eso sí) y el segundo no. En el apartado 4.6 discutiremos el origen del prepico en AgI, que señala la existencia de una ordenación de medio alcance. En su trabajo sobre AgI fundido mediante simulaciones ab initio, 128 Capítulo 4. Estudio de α αα α-AgI y AgI fundido PIM1e. En el caso de las configuraciones de los aniones, vemos que la distribución es muy parecida en los dos modelos, y mucho menos dispersa que en el caso de los cationes. En las configuraciones catiónicas, D(R V ) es notablemente más ancha en PIM1e que en RIM, y su máximo está en un valor mayor de R V . Es decir, las configuraciones de los iones de plata en PIM1e tienen en promedio espacios vacíos más grandes que en RIM, y comprenden una mayor variedad en lo que se refiere al tamaño de los voids. La causa de esta diferencia es que en PIM1e, dado que la distancia entre cationes es menor, hay zonas de alta densidad catiónica, lo que implica a su vez la existencia de otras vacías, es decir de voids más grandes. Por lo tanto, el diferente comportamiento de D(R V ) en los dos modelos es consecuencia de la polarización inducida, y pone de manifiesto la inhomogeneidad de la disposición espacial de los cationes en PIM1e. Figura 4.10: Distribución de probabilidad D(R V ) del radio de los voids entre aniones (curvas finas) y cationes (curvas gruesas) en las simulaciones MD de AgI a 923 K con RIM (curvas de puntos) y PIM1e (continuas). Dado que el prepico aparece en S AgAg (k), lo que nos interesa es estudiar la distribución de los cationes. Por ello, a partir de ahora nos centramos en estudiar los voids de las configuraciones catiónicas, y analizamos las correlaciones entre ellos. De manera análoga a las funciones de distribución radial entre partículas g ab (r), podemos definir la función de distribución radial entre voids, g VV (r). En la Figura 4.11 se presentan las g VV (r) de los voids entre cationes obtenidas con RIM y PIM1e. Estas funciones presentan un brusco decaimiento inicial, en lugar de ser nulas en el origen como las 2 3 4 5 R V (Å) 0 1 2 D ( R V ) I (RIM) I (PIM1e) Ag (RIM) Ag (PIM1e) Capítulo 4. Estudio de α αα α-AgI y AgI fundido 129 g ab (r), debido a que los voids no tienen volumen excluido y por tanto no existe un límite inferior a la separación entre ellos. Tras este decaimiento g VV (r) exhibe una oscilación formada por un pronunciado mínimo y por un máximo, desplazados en PIM1e a r mayores que en RIM porque el tamaño de los voids es en promedio mayor en el primer modelo. Figura 4.11: Función de distribución radial g VV (r) de los voids entre iones de plata en las simulaciones MD de AgI a 923 K con RIM (curva de puntos) y PIM1e (curva continua). La información de la estructura de los voids catiónicos en el espacio recíproco viene dada por el factor de estructura void-void, S VV (k), que está relacionado con g VV (r) y se define análogamente a S ab (k), V V V V V V VV V VV 1 1 V 1 ( ) exp[ ·( )] 1 F[ ( )] N N i j i j S k i g r N ρ = = = − − = + ∑ ∑ k r r , (4.1) donde N V es el número de voids en cada configuración y ρ V su densidad media, r iV y r jV son las posiciones de los voids, y F indica la transformada espacial de Fourier. Como se observa en la Figura 4.12, S VV (k) presenta un primer pico muy marcado. El mayor tamaño de los voids en PIM1e se traduce en que el primer pico de su S VV (k) está desplazado hacia la izquierda respecto al de RIM. El S VV (k) de PIM1e exhibe su pico principal en k ~ 1.1 Å -1 , posición similar a la del prepico de S AgAg (k) en PIM1e. Es decir, que la escala de longitud característica asociada a la estructura de los voids entre 0 2 4 6 8 10 12 r (Å) 1 2 g VV ( r ) RIM PIM1e 130 Capítulo 4. Estudio de α αα α-AgI y AgI fundido cationes es la misma que la señalada por el prepico de S AgAg (k). Por tanto, este resultado parece confirmar que el prepico obtenido en las simulaciones de AgI fundido con PIM1e tiene su origen en la ocupación inhomogénea del espacio por parte de los iones de plata, que dejan grandes voids. La longitud característica de las correlaciones entre estos voids está presente también en las correlaciones entre iones de plata, y queda reflejada en el prepico en k ~ 1Å -1 del S AgAg (k) de PIM1e. Figura 4.12: Factores de estructura catión-catión S AgAg (k) (panel superior), y void-void S VV (k) (panel inferior) obtenidos en simulaciones MD de AgI a 923 K para RIM (curvas de puntos) y PIM1e (curvas continuas). Adviértase la diferencia de escala vertical entre ambas figuras. Por otro lado, el S VV (k) de RIM tiene su pico principal alrededor de 1.5 Å -1 . Esta posición es cercana a k M , de modo que la contribución de los voids en la estructura de los cationes en RIM no se distingue del pico 1 2 S AgAg ( k ) 0 2 4 S VV ( k ) 0123456 k (Å -1 ) RIM PIM1e Capítulo 4. Estudio de α αα α-AgI y AgI fundido 131 principal de S AgAg (k) y no da lugar a ningún prepico. Lo mismo ocurre si se calcula el factor de estructura S VV (k) de los voids de las configuraciones de aniones, y no de cationes como hasta ahora. Los S VV (k) de los voids aniónicos tienen un pico en k M , lo que explica que el prepico no se observe en S II (k) en ninguno de los dos modelos. Figura 4.13: Factores de estructura S AgAg (k) (panel superior) y S VV (k) (panel inferior izquierdo), y función de distribución radial g AgAg (r) (panel inferior derecho) obtenidos en simulaciones MD de AgI a 923 K con diferentes valores de la polarizabilidad de los iones de iodo, desde α I = 0 (RIM, curvas de puntos) hasta α I = 6.12 Å 3 (PIM1e, curvas continuas), pasando por los diversos valores indicados en la leyenda (curvas de colores). Para confirmar la relación que hemos establecido entre el prepico en S AgAg (k) y las correlaciones entre los voids catiónicos, y el papel que la polarización juega en ella, hemos realizado simulaciones de AgI fundido con diversos valores de la polarizabilidad α I , cuyos resultados para g AgAg (r), 1 2 0 2 4 012345 k (Å -1 ) α I = 3.45 Å 3 α I = 4.5 Å 3 α I = 5.0 Å 3 α I = 5.5 Å 3 α I = 6.12 Å 3 (PIM1e) α I = 0 Å 3 (RIM) 2 4 6 8 10 12 r (Å) 0. 0 0. 6 1. 2 S AgAg ( k ) S VV ( k ) g AgAg ( r ) 132 Capítulo 4. Estudio de α αα α-AgI y AgI fundido S AgAg (k) y S VV (k) presentamos en la Figura 4.13. Vemos que a medida que α I aumenta, el primer pico de g AgAg (r) y el pico principal de S VV (k) se desplazan hacia la izquierda. Conforme el pico de S VV (k) se aleja de k M y se acerca a 1 Å -1 , su efecto en S AgAg (k) comienza a hacerse patente y aparece un shoulder que crece con α I , hasta que en un valor entre 5.5 Å 3 y 6.12 Å 3 , el shoulder pasa a ser un máximo, es decir un prepico. La misma conclusión se extrae de simulaciones de AgCl fundido con diferentes valores de α Cl , que no mostramos aquí. Cuando α Cl = 3.45 Å 3 , el pico de S VV (k) es cercano a k M y, como se ve en la Figura 4.9, no existe prepico. Sin embargo, si la polarizabilidad se aumenta hasta 4.0 Å 3 , el pico de S VV (k) se encuentra en ~ 1.25 Å -1 y en esa misma posición aparece un prepico en S AgAg (k). La observación experimental de que la intensidad del prepico en los haluros de cobre fundidos aumenta en la secuencia CuCl-CuBr-CuI es consistente con estos resultados. 4.7 Transporte iónico en AgI Tabla 4.2: Coeficientes de autodifusión D a y conductividades iónicas σ calculados en simulaciones MD de α-AgI a 573 K y AgI fundido a 923 K con RIM, PIM2sqµ y PIM1e. α αα α-AgI (573 K) AgI fundido (923 K) RIM PIM2sqµ µµ µ PIM1e RIM PIM2sqµ µµ µ PIM1e D Ag (10 −5 cm 2 /s) 2.34 a 2.70 a 3.02 a 6.50 7.24 9.08 D I (10 −5 cm 2 /s) 0 0 0 2.22 1.88 1.75 σ (Ωcm) −1 2.2 b 2.1 b 1.9 b 3.4 b 2.8 b 2.8 b σ ΝΕ (Ωcm) −1 1.2 1.3 1.5 2.5 2.6 3.1 ∆ −0.87 c −0.53 c −0.27 c −0.36 −0.08 0.10 a D Ag experimental: D Ag (573 K) = 2.41×10 −5 cm 2 /s [56]. b σ experimentales: σ (573 K) = 1.97 (Ωcm) -1 [11], σ (923 K) = 2.51 (Ωcm) -1 [33]. c ∆ experimental: ∆ (573 K) = -0.6. Capítulo 4. Estudio de α αα α-AgI y AgI fundido 133 Figura 4.14: Desplazamientos cuadráticos medios Q Ag (t) (líneas superiores) y Q I (t) (líneas inferiores) de α-AgI a 573 K para RIM (líneas de puntos), PIM2sqµ (líneas de trazos) y PIM1e (líneas continuas). En la Tabla 4.2 se recogen los coeficientes de transporte iónico obtenidos en las simulaciones de α-AgI y AgI fundido con RIM, PIM1e y PIM2sqµ. En los tres modelos, en la fase superiónica α los aniones se mantienen alrededor de sus posiciones en la red bcc, por lo que su coeficiente de autodifusión D I = 0, y los cationes se difunden a través de la red aniónica, es decir D Ag ≠ 0. Los desplazamientos cuadráticos medios Q a (t) representados en la Figura 4.14 ilustran con claridad el diferente comportamiento de las dos especies. Mientras el valor de Q I (t) se mantiene constante e igual a ~ 0.2 Å 2 tras su subida inicial, la pendiente de Q Ag (t) es distinta de cero. El mayor coeficiente D Ag en α-AgI es el de PIM1e, lo cual está relacionado con el mayor número de posiciones tetraédricas accesibles que los iones de plata tienen en este modelo. Como hemos visto en el apartado 4.4, la polarización inducida permite que un ión de plata en PIM1e pueda acercarse más a un segundo catión que en los otros modelos, de tal modo que dispone de más sitios tetraédricos en los que situarse y, por tanto, puede difundirse con mayor facilidad. Es decir, la polarización inducida en los aniones abre nuevos canales de difusión para los cationes, aumentando en cierto modo el espacio libre que estos tienen para difundirse. Pese a que la estructura de PIM1e es más acorde con los datos experimentales de α-AgI, el valor de 0 1 2 3 t ( p s ) 0 1 2 Q a ( t ) (Å 2 ) RIM PIM2sq µ PIM1e Q Ag ( t ) Q I ( t ) 134 Capítulo 4. Estudio de α αα α-AgI y AgI fundido D Ag obtenido con RIM es más cercano al experimental. De todas formas, como ya dijimos en el capítulo anterior, lo que más nos interesa de estos resultados es comparar el comportamiento de cada modelo para valorar el efecto de la polarización, más que el valor concreto de los coeficientes de difusión, que depende decisivamente de detalles del potencial efectivo a parejas. En AgI fundido sigue habiendo un marcado contraste entre la difusividad de las dos especies iónicas. El cociente D + /D - varía entre tres y cinco en función del modelo, y es mayor que en AgBr fundido (ver Tabla 3.2), principalmente porque en AgI la diferencia de tamaños entre cationes y aniones es más acentuada, y en el caso de los PIM porque además α I > α Br . D Ag es mayor en los PIM que en RIM debido a que, de manera análoga a lo que sucede en α-AgI, los momentos dipolares inducidos en los aniones permiten que aumente el espacio libre por el que los cationes pueden moverse y favorecen su difusión. El hecho de que D Ag en PIM2sqµ sea intermedio a los de RIM y PIM1e, y no menor como ocurría en el PIM1sq de AgBr estudiado en el capítulo anterior, indica que la influencia de la polarización en los resultados de PIM2sqµ está menos amortiguada que en el PIM1sq de AgBr ( λ = λ q = λ µ es menor en PIM2sqµ), además de que α I > α Br . Por su parte, el orden de los valores de D I en los tres modelos es el inverso al de D Ag , aunque las diferencias son más pequeñas. El que D I sea menor en los PIM que en RIM no es un efecto de la polarización, sino que está relacionado con el distinto tamaño efectivo de los aniones en los potenciales de ión rígido utilizados, que es mayor en el VR0 de los PIM que en el VR de RIM. La forma de las funciones de autocorrelación de velocidades C Ag (t) y C I (t) es, tanto en α-AgI como en la fase líquida, la típica de los haluros de plata y cobre, y ya fue discutida en los apartados 1.4.1 y 3.5.1. Es por ello que no las mostramos en este caso. La conductividad iónica σ de los PIM es menor que la de RIM en las dos fases estudiadas, y muestra un mejor acuerdo con los resultados experimentales. En ambas fases, los valores de σ están ordenados de manera inversa a σ NE = ρ N e 2 (D Ag + D I )/2k B T. Es decir, cuanto mayor es la suma de los coeficientes de autodifusión, menor es σ . Este resultado indica la relevancia de las correlaciones entre iones distintos en el cálculo de σ , Capítulo 4. Estudio de α αα α-AgI y AgI fundido 135 cuantificada por el parámetro ∆ = 1σ / σ NE . Tanto en α-AgI como en la fase fundida, el valor de ∆ para RIM es el más negativo de los tres modelos. En el otro extremo está PIM1e, cuyo ∆ es el menos negativo en α-AgI y positivo en AgI fundido. Este último es el único caso en que σ < σ NE . La tendencia a valores mayores de ∆, observada también en las simulaciones de AgBr fundido, parece ser pues una característica constante en los modelos de ión polarizable. A partir de los valores experimentales de D Ag y σ en la fase α se deduce ∆ = -0.6, cercano al obtenido con PIM2sqµ. Para AgI fundido no existen datos experimentales de D Ag que permitan calcular ∆. Figura 4.15: Función δ ZZ (t) de AgI fundido a 923 K para RIM (curva de puntos), PIM2sqµ (curva de trazos) y PIM1e (curva continua). La función δ ZZ (t), que da cuenta de las correlaciones entre velocidades de iones distintos y está relacionada con ∆ a través de la ecuación (1.46), se muestra en la Figura 4.15 para los tres modelos de AgI fundido simulados. Al comparar estos resultados con los de sistemas similares como AgCl y AgBr fundidos, se observa que, independientemente del modelo, las oscilaciones de δ ZZ (t) se hacen cada vez menos marcadas en la secuencia AgCl-AgBr-AgI. El comportamiento de la δ ZZ (t) de RIM para AgI es bastante diferente al de los PIM, y también al de las δ ZZ (t) obtenidas para AgCl y AgBr fundidos. A diferencia de estas, no muestra un máximo entre -0.2 0.0 0.2 δ ZZ ( t ) 0.0 0.2 0.4 0.6 0.8 t (ps) RIM PIM2sq µ PIM1e 136 Capítulo 4. Estudio de α αα α-AgI y AgI fundido 0.1 y 0.2 ps sino un mínimo. En cambio, las δ ZZ (t) de los PIM sí recuerdan a las calculadas para AgCl y AgBr tanto con RIM como con PIM. Las funciones δ ZZ (t) obtenidas en las simulaciones de la fase α son bastante similares a las de la fase fundida, sólo que sus oscilaciones son ligeramente más pronunciadas. 4.8 Conclusiones En este capítulo hemos constatado que, para evitar la polarización catastrófica en simulaciones MD de AgI en las que ambas especies iónicas son polarizables, es necesario utilizar un modelo PIM2sqµ. Es decir, que al tener en cuenta la polarización inducida por las interacciones deformadoras de corto alcance, no sólo hay que amortiguar la contribución del campo eléctrico de las cargas, sino también la del campo eléctrico de los dipolos. Sin embargo, la simulación de este tipo de modelos es muy costosa como consecuencia de la ralentización del cálculo de los momentos dipolares. Sería interesante estudiar si la variante utilizada por Madden et al., explicada en el primer párrafo de la página 69, es menos costosa y ofrece mejores resultados. Además de PIM2sqµ, hemos simulado un PIM1e y el modelo de ión rígido en que ambos se basan (RIM). Los tres modelos reproducen el comportamiento superiónico de α-AgI, en el que los iones de plata muestran una clara preferencia por ocupar las posiciones tetraédricas en las caras de la red bcc formada por los aniones. En los tres casos se observa que el movimiento difusivo de los cationes se produce mediante saltos entre posiciones tetraédricas adyacentes, mayoritariamente a través de los sitios trigonales, lo que está de acuerdo con los resultados experimentales. Pero mientras que en RIM y PIM2sqµ la distancia típica entre cationes vecinos es la máxima que puede darse (~ 4.3 Å), en PIM1e se reduce hasta 3.2 Å, de modo que en este último pueden ocuparse pares de sitios que en RIM y PIM2sqµ están prohibidos, lo que favorece la difusión. Los resultados del análisis RMC de datos de difracción de neutrones y los de simulaciones ab initio parecen indicar que la estructura obtenida en PIM1e guarda un mayor acuerdo con la realidad. El acercamiento entre cationes vecinos también se observa en las simulaciones de AgI fundido con PIM1e y, en menor medida, con PIM2sqµ. Capítulo 4. Estudio de α αα α-AgI y AgI fundido 137 Su causa es, como vimos en el capítulo anterior para AgBr fundido, el apantallamiento de la repulsión electrostática entre ellos que se da como consecuencia de la inducción de polarización en los aniones. El hecho de que la separación entre cationes vecinos sea menor que entre aniones provoca que la g AgAg (r) de PIM1e deje de estar en fase con g II (r), lo que en el espacio recíproco se refleja en la menor altura del pico principal de S AgAg (k) en k M ~ 1.7 Å -1 . Además, aparece otro máximo en k ~ 2.25 Å -1 , de altura similar a la del pico principal. Pero sin duda la característica más llamativa del S AgAg (k) de AgI fundido es el prepico en k ~ 1 Å -1 que se obtiene con PIM1e, que revela la existencia de una escala de longitud característica de medio alcance, y que no se observa ni con RIM ni con PIM2sqµ. El prepico de S AgAg (k) se manifiesta también en el S c (k) de PIM1e, que es el único factor de estructura simulado que reproduce el prepico observado en los datos obtenidos mediante la difracción de neutrones. Así pues, las simulaciones permiten concluir que el prepico observado experimentalmente es debido a las correlaciones catióncatión y es un efecto de la polarización inducida en los iones de iodo. La aparición del prepico en PIM1e parece ser debida a la inhomogeneidad de la disposición espacial de los cationes, ya que la reducción de la separación entre ellos debida a la polarización provoca la aparición de zonas de alta densidad de cationes que coexisten con otras vacías, o voids, de gran tamaño. El orden de medio alcance en la estructura de los iones de plata, que resulta muy difícil de detectar en la función de distribución radial g AgAg (r), está presente de forma mucho más evidente en la estructura de los voids entre cationes. Así pues, el prepico de S AgAg (k) obtenido con PIM1e refleja la longitud característica de las correlaciones entre voids. El que la polarización inducida en los aniones esté en el origen del prepico concuerda con la observación de que la intensidad del prepico en los factores de estructura experimentales de los haluros de cobre aumenta con la polarizabilidad aniónica [57]. Como resultado de la polarización inducida, los iones de plata tienen más espacio libre para moverse, y los coeficientes de autodifusión de los cationes tanto en la fase α como en la líquida son mayores para PIM1e que para PIM2sqµ y RIM (en este orden). A pesar de ello, la conductividad iónica es más baja en los modelos de ión polarizable que en RIM, y se ajusta mejor al valor experimental. 240 Apéndice F. Derivación del teorema de fluctuación-disipación para χ L (k)