Difracción de rayos X

Facultad de Física Grado en Física Difracción de rayos x: Resolución de estructuras cristalinas y problema de las fases. Trabajo de Fin de Grado Autor: Paula Llamas Martínez Tutor: Alberto Criado Vega Junio 2021 Índice 1. Objetivos 4 2. Introducción 5 3. Difracción por un cristal 7 3.1. Espacio recíproco y ecuaciones de Laue . . . . . . . . . . . . . . . . . . . . 10 4. Problema de las fases 11 5. Métodos del espacio directo 13 5.1. FuncióndePatterson.............................. 13 5.1.1. Líneas y secciones Harker . . . . . . . . . . . . . . . . . . . . . . . 16 5.1.2. Método del átomo pesado . . . . . . . . . . . . . . . . . . . . . . . 17 5.2. Reemplazoisomorfo............................... 18 5.2.1. Estructura no centrosimétrica . . . . . . . . . . . . . . . . . . . . . 19 5.3. Métodos de modificacion de la densidad electrónica . . . . . . . . . . . . . 22 5.3.1. Corrección de los errores de truncamiento en las series de Fourier . 23 5.3.2. Inversióndecarga............................ 26 6. Métodos directos 28 6.1. Invariantes y semi-invariantes de estructura . . . . . . . . . . . . . . . . . 30 6.2. Métodos probabilísticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6.2.1. Teoría de la representación . . . . . . . . . . . . . . . . . . . . . . . 33 6.3. Procedimiento para la obtención de fases mediante métodos directos. . . . 36 7. Resumen y conclusiones 40 Anexos 41 A. Teorema de convolución 41 B. Factor de estructura normalizado 41 C. Probabilidad del triplete 42 1. Objetivos El principal objetivo de la física es el conocimiento y explicación de los fenómenos del universo que nos rodea. El universo está formado por materiales, que a su vez están compuestos por átomos. El estudio de estos materiales es fundamental tanto como para explicar diversos fenómenos, así como para estudiar sus diversas aplicaciones útiles para el ser humano. En este documento nos vamos a centrar en la estructura de aquellos materiales cuyos átomos se disponen de manera ordenada y periódica en el espacio, los denominados cristales. Debido a su particular estructura, el estudio de los cristales se puede reducir a una pequeña porción denominada celda unidad, que tiene la propiedad de que su repetición en el espacio define todo el sólido. Así, estudiando la distribución en la celda unidad de los distintos átomos que conforman la estructura, podremos resolver el cristal. Para esto, es necesario realizar experimentos de difracción utilizando rayos X como emisión incidente. El resultado va a ser un patrón de difracción a partir del cual vamos a poder hallar las dimensiones de la celda unidad y analizar las simetrías existentes en el sólido. Sin embargo, para resolver completamente la estructura, es necesario hallar el factor de estructura, proceso que conlleva una mayor dificultad, debido a lo que se denomina el problema de las fases El objetivo de este proyecto es la realización de una revisión bibliográfica a través de diversos autores y publicaciones cristalográficas para recopilar los principales métodos de resolución del problema de las fases. Veremos que estos se dividen en dos tipos, dependiendo de si su aplicación es en el espacio real o en el recíproco. Estudiaremos los distintos métodos y procedimientos de cada tipo, así como las principales relaciones matemáticas necesarias. 4 2. Introducción A finales de 1895, Wilhelm Conrad Röntgen, físico alemán, estaba realizando experimentos con un tubo de rayos catódicos en la Universidad de Würzburg (Alemania) cuando observó que, al hacer atravesar el tubo (previamente cubierto con un cartón negro para eliminar la luz visible) una corriente de alta tensión, se apreciaba una luz fluorescente en una pantalla de platino-cianuro de bario [1, 2]. Investigó este efecto, colocando la pantalla a distintas distancias e interponiendo objetos entre ésta y el tubo, y descubrió que esta nueva radiación podía atravesar objetos de cierto espesor, y que era posible visualizar un espectro de los huesos y músculos de su mano en la pantalla. Algunos metales como el plomo o el platino también producían sombras en la pantalla que eran capaz de registrarse de manera gráfica en láminas fotográficas. Röntgen produjó así la que se considera la primera radiografía, exponiendo durante 15 minutos la mano de su mujer y registrando la imagen en una placa fotográfica. Por su carácter desconocido, a esta nueva radiación se la denominó rayos X, y por su descubrimiento, Röntgen recibió el Premio Nobel de Física en 1901. Pocos años después, se creó en la universidad de Munich un grupo de investigación para el estudio de los rayos X. Uno de sus objetivos, entre otros, era la confirmación de su carácter ondulatorio a partir de su longitud de onda. Para esto, había que hacer pasar esta radiación por una rejilla cuya distancia de separación fuera comparable a su longitud de onda, y estudiar el patrón de difracción resultante. El problema era la dificultad existente para fabricar una rejilla con las dimensiones necesarias, ya que los experimentos demostraban que la longitud de onda de los rayos X era muy pequeña. La solución vino de mano del físico alemán Max von Laue, quién sugirió que se utilizara un cristal como rejilla. Esto confirmaría, además de la naturaleza ondulatoria de los rayos X, la presunción de la disposición ordenada y periódica de los átomos en un sólido cristalino. Este experimento fue realizado con éxito por Walter Friedrich y Paul Knipping, los cuales colocaron un cristal de sulfato de cobre frente a un haz de rayos X colimado. Los resultados se verían en una placa fotográfica colocada detrás del cristal. El resultado del experimento (fig. 1) muestra una mancha central, correspondiente a la incidencia directa del rayo, rodeada de otras pequeñas manchas de menor tamaño, que corresponderían al 5 haz difractado, confirmandóse así ambas hipótesis. Esto hizo que en 1914, Laue obtuviese el Premio Nobel de Física. Figura 1: Imagen de la primera difracción de rayos X. Firma de Max von Laue [3]. Una vez demostrada la periodicidad de las posiciones atómicas en un cristal, el siguiente paso era lograr obtener estas posiciones. En 1913, William Henry Bragg y su hijo William Lawrence Bragg propusieron un modelo para hallar la estructura a partir del patrón de difracción y presentaron la famosa ley de Bragg, con la cual se obtenían las condiciones de difracción. 2dsin(θ) = nλ (1) Como los cristales eran periódicos, estipularon que los rayos X se reflejaban en ’planos’ formados por átomos, separados una distancia d(fig. 2). La diferencia de camino recorrido entre las ondas reflejadas en los puntos B y D vendría dada por 2dsen(θ), de manera que para que la interferencia de ambas sea positiva, debía ser igual a un número entero de longitudes de onda. Así, consiguieron resolver estructuras sencillas como NaCl, ZnS (blenda), o la del diamante, consiguiendo el Premio Nobel de Física en 1915. Nació así la disciplina física que se conoce como cristalografía de Rayos X. Actualmente, existen numerosas técnicas de difracción y métodos para determinar las estructuras de los materiales que nos rodean. Se han obtenido las estructuras de más de un millón de compuestos, desde moléculas 6 sencillas hasta macromoléculas [4, 5]. Figura 2: Reflexión de radiación por dos planos del cristal que están a una distancia d. La dirección de la radiación incidente viene representada por ~s0, y la de la reflejada por ~s [6]. El proceso de resolución de una estructura cristalina consta de cuatro pasos. El primero de ellos es la obtención del espectro de difracción, a partir del cual se puede hallar la celda unidad del cristal, siendo este el segundo paso. Posteriormente se analizarían las simetrías existentes con el objetivo de hallar el grupo espacial del cristal. El último paso sería la obtención de las posiciones atómicas, quedando totalmente resuelta la estructura. Los tres primeros pasos son relativamente sencillos, pero el último presenta una mayor complejidad debida al problema de las fases. A lo largo de esta memoria, vamos a analizar este problema, así como las distintos métodos para su resolución. 3. Difracción por un cristal Tenemos que un cristal actúa ante radiación X incidente como una malla tridimensional difractora, de manera que cuando se produce la interacción entre el cristal y el haz, obtenemos un patrón de difracción acorde con la disposición de los átomos existentes en la muestra. A partir de su estudio, es posible hallar la estructura del cristal. En primer lugar, vamos a analizar la interferencia de las ondas difractadas. De este análisis vamos a poder eliminar el tiempo ya que el patrón de difracción es estacionario. 7 También vamos a aproximar a ondas planas la radiación incidente, ya que el foco emisor de los rayos X va a estar lo suficientemente alejado del cristal. Así, comenzamos estudiando una situación simple en la que tenemos dos centros dispersores (fig. 1) OyO0separados una distancia |~r|=rsobre los que incide un haz de rayos de longitud de onda λ. Figura 3: Comportamiento de dos centros dispersores OyO0, relacionados por ~r ante radiación incidente. La dirección de la radiación incidente viene representada por ~s0, y la de la reflejada por ~s [6]. Los rayos incidentes tienen la siguiente forma, según la aproximación previamente considerada: ~u =Aei2π ~s0·~r (2) El haz incide con una dirección ~s0y es dispersado con una dirección ~s, de manera que la diferencia de fase entre las ondas difractadas por ambos puntos, siendo ~r =~ OO0, es: δ= 22π λsin(θ)r=2π λ(~s0−~s)·~r =2π λ~ r∗·~r (3) Así, si tomamos como cero la fase de la onda difractada por O, tenemos que la difractada por O0,~vO0, estará descrita por: ~vO0=fO0ei2π~ r∗·~r (4) Si nos alejamos de este caso puntual y extrapolamos este resultado a un caso en el que haya Ncentros dispersores, como un cristal, obtenemos que el resultado de la dispersión se puede recoger en: F(~ r∗) = N X j=1 fjei2π~ r∗·~rj=|F|eiϕ (5) 8 Donde F(~ r∗)se conoce como factor de estructura, y fjes el factor de dispersión o factor de forma atómica. fjexpresa la capacidad que tiene un átomo aislado para dispersar la onda que incide sobre él. Se puede expresar como: f= 4πZ∞ 0ρ(r)r2sin 2πrr∗ 2πrr∗dr (6) Donde ρ(r)es la densidad electrónica. El factor sin 2πrr∗ 2πrr∗modela la interferencia entre los electrones del átomo [6, 7]. Por otro lado, una red tridimensional infinita se puede representar matemáticamente como: L(~r) = +∞ X u,v,w=−∞ δ(~r −~ru,v,w)(7) Observamos una función nula en todos los puntos excepto en aquellos en los que ~r =~ru,v,w, siendo ~ru,v,w =u~a +v~ b+w~c. Los vectores ~a,~ by~c serían los vectores que conforman una celda de la red, y como ésta es periódica, u,v,w∈Z. Si ahora definimos ρM(~r)como la densidad electrónica de una celda unidad del cristal, tendremos que la densidad electrónica del cristal completo vendría dada por la convolución de L(~r)yρM(~r). ρ(~r) = L(~ r)∗ρM(~r)(8) En nuestro caso, vamos a tener que los centros dispersores van a ser los electrones del sólido (nube electrónica), de manera que introducimos la densidad electrónica del sólido en la ec. 5, que se transforma en: F(~ r∗) = ZVρ(~r)exp(2πi~ r∗~r)d~ r∗(9) Obtener la densidad electrónica implica obtener las posiciones atómicas, por lo cual fijamos así el objetivo final de la cristalografía. Si nos fijamos en la expresión anterior, observamos que F(~ r∗)coincide con la transformada de Fourier de ρ(~r),T[ρ(~r)] por lo que: ρ(~r) = ZV∗F(~ r∗)exp(−2πi~ r∗~r)d~ r∗ =1 V ∞ X h,k,l=−∞ Fhklexp[−2πi(hx +ky +lz)] =T−1[F(~ r∗)] (10) Mediante la transformada inversa de Fourier del factor de estructura, podremos obtener la función de densidad electrónica [6, 8, 9]. En la siguiente imagen podemos observar lo que significa físicamente lo descrito anteriormente. 9 En algunos casos también va a ser necesario eliminar el pico del origen, por la posibilidad de que éste oculte algún máximo con ~u pequeño. Para esto, volvemos a la expresión para el factor de estructura ec. 13, y suprimimos el término que corresponde a la contribución individual de cada átomo PN j=1 f2 j, dejando el término de interacción con los átomos vecinos. |F0 ~ h|2=|F~ h|2− N X j=1 f2 j= 2 N X j>k fjfkcos2π~ h(~rj−~rk)(21) 5.1.1. Líneas y secciones Harker Vamos a introducir ahora los conceptos de líneas y secciones Harker, que nos van a permitir analizar las simetrías presentes en el cristal a través de la función de Patterson [6, 12]. En 1935, David Harker, estudiando diagramas de Patterson, descubrió que existían zonas en el mapa en las cuales se agrupaban los máximos. Esto era debido a los vectores interatómicos relacionados entre sí por las simetrías presentes en el cristal. Dependiendo del elemento de simetría considerado, vamos a tener líneas o secciones en las cuales vamos a observar esta acumulación de máximos. Por ejemplo, las posiciones x,y,zyx,¯y,z, que son posiciones equivalentes en el grupo espacial Pm, en el cual tenemos un plano de reflexión como elemento de simetría, van a dar lugar a la línea Harker (0,2y, 0): u=x−x= 0; v=y−¯y= 2y;w=z−z= 0 Como en el método de Patterson los vectores interatómicos van a ser trasladados al mismo origen, la mayor parte de las traslaciones no van a poder ser identificadas. Así, mediante este métodos, vamos a poder identificar sólo 24 grupos espaciales de los 230 existentes. En la siguiente tabla se muestran las lineas y secciones Harker para los elementos de simetría más comunes [6, 7]: 16 Elemento de simetría Línea/sección Harper Eje 2 ka,b,c (0,v,w);(u,0,w);(u,v,0) Eje 21ka,b,c (1/2,v,w);(u,1/2,w);(u,v,1/2) Plano m ⊥a,b,c (u,0,0);(0,v,0);(0,0,w) Plano de deslizamiento a ⊥b,c (0,v,1/2);(1/2,0,w) Plano de deslizamiento b ⊥a,c (u,1/2,0);(0,1/2,w) Plano de deslizamiento c ⊥a,b (u,0,1/2);(0,v,1/2) Cuadro 1: Principales elementos de simetría y líneas y secciones Harker asociados. 5.1.2. Método del átomo pesado Tenemos que la altura de los máximos del mapa de Patterson es proporcional al producto de los números atómicos Z de los átomos involucrados. Existen situaciones en las cuales vamos a observar diferencias relativamente grandes entre estas alturas, lo que va a llevar a una simplificación del estudio. Este es el caso en el cual la estructura contenga un número límitado de átomos pesados. Así, estos átomos van a suponer una mayor contribución a la altura de los picos, por lo cual se va a poder identificar sus coordenadas con más facilidad. Por supuesto, también existen límites. Si la estructura tiene un número demasiado alto de átomos pesados, también vamos a encontrar nuevas dificultades en el diagrama de Patterson, ya que su contribución va a anular la de los átomos más pequeños. Así, el valor óptimo estipulado es cuando el cociente entre el sumatorio de los cuadrados del número atómico de los átomos pesados (Z2 p) y ligeros (Z2 l) es igual a 1. PZ2 p PZ2 l = 1 (22) Esta información va a ser sumamente útil a la hora de crear el modelo inicial, ya que los valores de la fase ψc ~ hde |Fc ~ h|van a ser bastante cercanos a los reales. Este modelo va a ser posteriormente computarizado y mejorado en ciclos mediante el método de síntesis sucesivas de Fourier. Este método funciona de la siguiente manera: Realizando una síntesis de Fourier al |Fc ~ h|calculado a partir de las posiciones de los átomos pesados, vamos a obtener nuevas posiciones atómicas. Con las coordenadas de estas nuevas coordenadas vamos a mejorar el modelo, que a su vez nos va a descubrir más posiciones, y así, sucesivamente. Es decir, nos encontramos ante un proceso iterativo que solo termina cuando todas las coordenadas atómicas son conocidas [6, 7]. 17 5.2. Reemplazo isomorfo Figura 7: Moléculas isomorfas, que difieren entre sí solo en que el grupo metilo en la primera (arriba) es reemplazado por un átomo de Cloro en la segunda (abajo) [7]. Este método [7, 10, 13, 14] hace uso del de Patterson para resolver estructuras isomorfas. Estas son aquellas que tienen una estructura molecular muy similar y que solo difieren en algún o algunos compuestos. Por ejemplo, las moléculas de 2-amino-4-metilo- 6-cloropirimidina y 2-amino-4,6-dicloropirimidina, en la que el (CH3) de la primera es sustituido por Cl en la segunda. Una de las características que más nos van a interesar de estas moléculas es que la estructura cristalina que forman también va a ser muy similar, teniendo la misma celda unidad y con dimensiones muy parecidas. Por ejemplo en el caso de las moléculas citadas previamente, ambas pertenecen al grupo espacial P21/a, con: a b c β Metil 16,43 ˚ A4,00 ˚ A10,31 ˚ A109º Cloro 16,45 ˚ A3,85 ˚ A10,28 ˚ A108º Cuadro 2: Dimensiones de la celda unidad para una molécula de 2-amino-4-metilo-6- cloropirimidina y otra de 2-amino-4,6-dicloropirimidina [6]. Cabe esperar por tanto, que en un experimento de difracción, obtengamos también un patrón muy similar para ambas moléculas. 18 Con lo anterior en mente, vamos a imaginar que tenemos dos estructuras isomorfas que difieren en que nátomos tipo A de la primera son sustituidos por nátomos de tipo B en la segunda. El factor de estructura de ambas será el mismo salvo por un término que se refiera a los nátomos que tienen distintos (fAofB), es decir: (F~ h)A=K~ h+fA n X j=1 exp(2πi~ h~rj) (F~ h)B=K~ h+fB n X j=1 exp(2πi~ h~rj) (23) Si restamos ambas expresiones, y denotando el sumatorio sobre los nátomos ’reemplazados’ como Ch, llegamos a la siguiente expresión. (F~ h)A−(F~ h)B= (fA−fB)Ch(24) |F~ hA|y|F~ hB|,fAyfBson conocidos. Si la posición del átomo reemplazado también es conocida (por ejemplo si su posición está restringida por simetrías, etc), podemos hallar Ch, de manera que encontraríamos (F~ h)Ay(F~ h)B. En el caso de estructuras centrosimétricas, ya habríamos encontrado la fase, ya que el factor de estructura para estos casos es real. Así, si F=|F|,φ= 0 y si F=−|F|,φ=π. En el caso en el que nuestra estructura sea centrosimétrica, pero no conozcamos Ch, debemos utilizar el método de Patterson para encontrar la fase. En primer lugar, debemos suponer que, en la mayoría de casos, el hecho de que (F~ h)Ay(F~ h)Btengan signos opuestos es poco habitual, ya que solo sería posible si el átomo reemplazado fuese un centro dispersor muy potente. Tendríamos entonces que (F~ h)A−(F~ h)B=±(|F~ hA|−|F~ hB|). De esta manera, lo único que podemos hallar a partir de la ec.24 es el módulo de Ch, pero no su signo. Sin embargo, la función de Patterson de la diferencia entre los factores de estructura resultará en un mapa en el cual solo encontraremos picos en posiciones relacionadas con los átomos reemplazados, ya que el resto de picos, al ser iguales entre ambas estructuras, se cancelarán entre sí. La interpretación de este mapa será más sencilla, de manera que podremos encontrar las posiciones atómicas de átomo reemplazado, pudiendo calcular Ch. 5.2.1. Estructura no centrosimétrica En primer lugar, suponemos que nuestra celda solo tiene un átomo reemplazado. Si elegimos este átomo como el origen de la celda unidad, el sumatorio de la ec.23 resulta 19 ser 1, (ya que la exponencial estaría elevada a 0 y nsería 1), luego podríamos reescribirla de la siguiente manera: (F~ h)A=K~ h+fA (F~ h)B=K~ h+fB (25) Estas ecuaciones, gráficamente, tendrían la forma de dos triángulos, uno de lados |(F~ h)A|, K~ hyfA, y otro de lados |(F~ h)B|,K~ hyfB, como se puede observar en la fig.8 a). Si conociesemos los tres lados de este triángulo, podríamos hallar las respectivas fases, que vendrían dadas por el ángulo formado por el lado |(F~ h)|yfen cada caso. Sin embargo, no conocemos K~ h, luego necesitaremos utilizar otro método. Figura 8: a) Relación gráfica de los factores de estructura de dos moléculas isomorfas centrosimétricas con fyK. b) Método de construcción de los triángulos anteriores para hallar la fase [7]. Sobre una recta, tomamosun punto Pa partir del cual vamos a definir dos segmentos OP YO0Pde longitudes fAyfB, respectivamente. Trazamos a continuación dos circunferencias de centros OyO0y radios |(F~ h)A|y|(F~ h)B|. Estas dos circunferencias van a cortarse en dos puntos QyR, de manera que la unión de los finales de los segmentos OP YO0Pcon este punto dará como resultado dos pares de los triángulos que mencionabamos, obteniendo así para cada uno dos valores de fase, de misma magnitud pero de distinto signo, fig.8 b). Lo que sí podemos deducir es la parte real de la fase, ya que sería la proyección sobre la recta formada por los puntos P,OyO0, que formarían la recta real. Mediante el método de sintesis de Fourier vamos a crear un modelo utilizando la parte real de la fase, obteniendo un mapa en el cual vamos a tener una estructura centrosimétrica. El mapa de Fourier tiene el doble de picos que átomos y para cada par de 20 picos relacionados centrosimétricamente, se debe seleccionar uno y rechazar el otro. Sin embargo, conociendo las propiedades quimicas del cristal, que nos dan una idea inicial de la distribución estructural de éste, podemos seleccionar los picos adecuados y resolver la estructura. Ahora nos trasladamos a la situación en la que hay nátomos reemplazados. En este caso, no tendríamos una expresión tan simple como la de la ec. 25, si no que volveríamos a la ec. 23. Recogiendo el término del sumatorio en un término χ, lo que tendríamos es: (F~ h)A=K~ h+χA (F~ h)B=K~ h+χB (26) Si en este caso contruimos el mapa de Patterson producto de la diferencia entre los factores de estructura, tendremos picos debidos a los vectores de los átomos reemplazados consigo mismos, Z2 A−Z2 By a la interacción de estos con el resto, (ZA−ZB)Zj. Si comparamos ahora con los mapas de Patterson individuales, podremos encontrar la contribución de los nátomos AyB, encontrando así χAyχB. El procedimiento para hallar la fase, sería análogo al anterior, con la diferencia de que en este caso, los segmentos OP yO0P, que en este caso tendrían de módulo χAyχB, no caerían sobre la recta real, viniendo la fase de esta contribución, φχ,~ hdada por el ángulo que forman estos segmentos con ella (fig.9 a)). Figura 9: a) Relación gráfica de los factores de estructura de dos moléculas isomorfas no centrosimétricas con χyK. b) Método de construcción de los triángulos anteriores para hallar la fase [7]. En este caso, la fase de K~ hvendría dada por φχ±∆φ~ h. Utilizando argumentos trigo- 21 nométricos, ya que estamos trabajando con triángulos, podemos hallar ∆φ~ h. Así, si nos fijamos en la fig.9 b), cos ψ=(|χB,~ h|−|χA,~ h|)2+|FB,~ h|2−|FA,~ h|2 2|FB,~ h|(|χB,~ h|−|χA,~ h|)(27) |K~ h|2=|FB,~ h|2+|χB,~ h|2−2|FB,~ h||χB,~ h|cos ψ(28) Así, ∆φ~ h= arc cos |FB,~ h|2−|K~ h|2−|χB,~ h|2 2|K~ h||χB,~ h|!(29) Como antes se ha mencionado, tendríamos dos soluciones posibles para la fase de K~ h. Si calculamos una síntesis de Fourier utilizando ambos valores: δ(~r) = 1 VX ~ h|K~ h|ncos(2π~r ·~ h−φχ−∆φ~ h) + cos(2π~r ·~ h−φχ+ ∆φ~ h)o(30) Cada término tendrá una distinta contribución, de manera que el correcto creará una imagen de la estructura mientras que el otro formará un ruido de fondo. Así, podremos identificar el correcto y hallar la fase. 5.3. Métodos de modificacion de la densidad electrónica Una vez hemos obtenido una estimación de la fase, ya podemos hacer un modelo inicial de la densidad electrónica. Habitualmente, suelen tener distintos errores, siendo el más importante de ellos una pobre resolución, de manera que se pierde información acerca de los picos menos intensos. Es interesante destacar que los mapas están constituidos por pixeles, de manera que se asemejan a una malla de espaciado muy pequeño. Va a ser ese espaciado el que va a determinar la calidad del mapa, de manera que si este es relativamente pequeño, tendremos uno de alta resolución, mientras que si es más grande estaremos hablando de una baja resolución. Otro problema es que en estos mapas iniciales, suelen existir regiones negativas, cuando la densidad electrónica siempre es mayor que cero [6, 15, 16, 17]. Con el objetivo de mejorar esta función de densidad, se propusieron métodos denominados como métodos de modificación de la densidad electrónica, que consisten en realizar ciclos que contienen los siguientes procesos : 22 1. Modificar la función de densidad electrónica. 2. Hallar su transformada de Fourier para obtener las fases y amplitudes. 3. Sustituir las amplitudes por las experimentales o por las calculadas mediante métodos directos y calcular el factor de estructura. 4. Obtener su transformada inversa de Fourier para obtener así la densidad electrónica. Dependiendo del tipo de moléculas con las que resultan mas eficientes, podemos distinguir dos tipos dentro de estos métodos. Si queremos resolver estructuras orgánicas, se utiliza un método que consiste en corregir los errores de truncamiento de las series de Fourier. Si la estructura es inorgánica, el método de inversión de carga es el óptimo. Antes de comenzar a describir los distintos procesos, es útil distinguir entre los modos de obtención de las fases iniciales. Así, si las utilizadas en el proceso no han sido sometidas a ningún proceso previo, se denominan ab initio. Si no tienen esta distinción, significa que han sido obtenidas a partir de otro proceso, por ejemplo el que hemos estudiado previamente, el reemplazo isomorfo. 5.3.1. Corrección de los errores de truncamiento en las series de Fourier Los errores de truncamiento son aquellos debidos al hecho de que estamos trabajando con aproximaciones, en vez de con datos o procedimientos matemáticos exactos. La idea es la siguiente, con el método de síntesis sucesivas de Fourier lo que obtenemos como resultado es una aproximación del factor de estructura del material a partir del cual hallamos la densidad electrónica. Debido a las limitaciones físicas dadas por el propio experimento, al realizar las series de Fourier se produce una acumulación de errores, de manera que las posiciones atómicas van a estar desplazadas de su lugar real y rodeadas de una serie de ondas positivas y negativas que serían como un fondo de ruido. Hay dos estrategias para solucionar esto: Free-lunch algorithm. Este método [6, 18, 19] se utiliza sobre todo en la resolución estructural de proteínas. Como se ha mencionado anteriormente, una baja resolución es uno de los problemas a los 23 que nos encontramos a la hora de analizar nuestro mapa de densidad electrónica. Esto hace que se pierdan contribuciones al factor de estructura (amplitud y fase) de estos picos más pequeños que no es posible distinguir. Utilizando un algoritmo basado en la extrapolación de los módulos y fases que faltan, cuyo objetivo será el de reducir el error de las reflexiones medidas, obtendremos un mapa con mayor resolución que nos permita distinguir todas las posiciones atómicas. El nombre free-lunch se debe a que no es necesario realizar ningún otro experimento o medida aparte para utilizarlo. Para explicar su funcionamiento, utilizaremos una estructura simple, en la que tendremos una celda unidad con a= 10˚ A, que contiene dos átomos de Mg, uno de O, otro N y dos de C. La función de densidad ρ(x)correcta, viene dada por la curva negra en la fig. 10 y suponemos que ha sido obtenida a partir de unos datos con una resolución de 0,5˚ A. Figura 10: Curvas de densidad electrónica para una estructura formada por dos átomos de Mg, uno de O, otro N y dos de C. La línea negra representa la función correcta, construida a partir de un experimento de difracción con una resolución de 0,5˚ A. La roja es la construida a partir de un experimento con una resolución de 1,7˚ A, que ha sido posteriormente mejorada por: un algoritmo de modificación de la densidad electrónica, (línea azul); un algoritmo de modificación de la densidad electrónica combinado con un algoritmo de extrapolación de los máximos no observados (línea verde) [18]. Nos vamos al caso en el que tengamos una peor resolución, de 1,7˚ A. A partir de estos 24 datos, la curva que vamos a obtener es la de color rojo en la figura anterior. Observamos regiones en las que ρ(x)<0, los picos están menos pronunciados y se solapan entre sí, y se ha perdido el segundo átomo de C, C2. El método a utilizar será la combinación de dos algoritmos de manera cíclica, el primero de ellos se va a encargar de modificar ρ(x) de acuerdo con la siguiente función, donde jes el número de ciclo: ρ(x) =              0si ρ(x)<0 ρ(x)1,3si ρ(x)>0, j ≤3 ρ(x)si ρ(x)>0, j > 3 (31) La función de la primera de las condiciones es eliminar las zonas negativas, y de la segunda, destacar los máximos. Sometiendo la curva roja a 15 iteracciones ρ→ϕ→ρ, obtenemos la curva azul. Seguimos teniendo las mismas dificultades, con la única diferencia de que se han eliminado casi totalmente los tramos negativos. Si comenzamos de nuevo el proceso de 15 ciclos, introduciendo el algoritmo NMRE (acrónimo de non-measured reflection extrapolation), que actúa en el medio ciclo ρ→ϕintroduciendo máximos extrapolados además de los que ya teníamos, obtenemos la curva verde de la fig. 10. Como podemos observar, esta curva si que identifica correctamente las posiciones atómica, pero tiene como desventaja la producción de mínimos falsos, aunque de baja amplitud. Los módulos del factor de estructura extrapolados vienen dados por la ec. 32, donde |E~ h(cal)|se calcula mediante la inversión del mapa de densidad. También se le asignan distintos pesos a los máximos dependiendo de si son los observados o los extrapolados. |E~ h(ext)|2= 1 + σ2 A(|E~ h(cal)|2−1) (32) Dependiendo de si la fase se ha obtenido ab initio o no, la aplicación del método varía. Para estructuras en las que las fases se han obtenido ab initio, se suele realizar un número alto de ciclos, entre 50 y 70. Si la fase ha sido refinada por un tratamiento anterior, se suelen realizar unos 30. Resolución por corrección. Debido al caracter finito de las dimensiones del espacio recíproco que utilizamos en la práctica (idealmente sería un espacio infinito), tenemos que la función de densidad real 25 Si conocemos distintos pares ϕ~ kj, ϕ~ h−~ kjque se refieren al mismo ϕ~ h, tendremos que la probabilidad total vendrá dada por el productorio de cada distribución r. Lo podremos expresar como, P(ϕ~ h) = r Y j=1 Pj(ϕ~ h) = Aexp   r X j=1 G~ h~ kjcos (ϕ~ h−ϕ~ kj−ϕ~ h−~ kj) (56) donde Aes un factor de normalización. Utilizando argumentos trigonométricos, podemos expresar la probabilidad como, P(ϕ~ h) = Aexp[α~ hcos (ϕ~ h−β~ h)] (57) Donde hemos definido ωj=ϕ~ kj+ϕ~ h−~ kj, α~ h=    r X j=1 G~ h~ kjcos ωj  1/2 +  r X j=1 G~ h~ kjsin ωj  1/2  (58) y tan β~ h=Pr j=1 G~ h~ kjsin ωj Pr j=1 G~ h~ kjcos ωj (59) Esta última expresión se denomina como la fórmula de la tangente y juega un rol muy importante en la determinación de la fase. De hecho, nos da el valor más probable de ϕ~ h, ϕ~ h≃β~ h, con una varianza que depende de α~ h. Para el caso de estructuras centrosimétricas, tendríamos, P(+) = 1 2+1 2tanh(σ3σ−3/2 2|E~ kE~ h−~ kE~ h|)(60) Al igual que anteriormente, consideramos el caso en el que tengamos varios pares S(~ k), S(~ h− ~ k)que se refieran al mismo S(−~ h), de manera que S(~ h)≃S(~ kj)S(~ h−~ kj)j= 1,2, ..., r (61) Lo descrito anteriormente conforma la teoría básica de probabilidad, la cual se ha ido desarrollando y mejorando. Aunque lo que hemos visto está aplicado al triplete, nos da una idea general del método a seguir. Así, el objetivo es encontrar una distribución de probabilidad que relacione las fases invariantes o relaciones invariantes entre ellas, Φ, con las amplitudes observadas. P(Φ|{|E|})(62) 32 El símbolo | indica que todas las amplitudes de {|E|} son conocidas. El principio anterior lo podemos aplicar al caso del cuarteto de dos maneras distinas. La primera de ellas vino de mano de Simerska en 1965. Dada la relación Qentre las fases, se da que el valor más probable de esta es 0, Q=ϕ~ h+ϕ~ k+ϕ~ l+ϕ−~ h−~ k−~ l= 0.(63) El problema de este método es que la varianza, que depende inversamente de B, B=2 N|E~ hE~ kE~ lE~ h+~ k+~ l|(64) tiene un valor muy alto debido al factor 1/N, lo que hace que sea poco fiable y por lo tanto no tenga interés práctico. Por otro lado, podemos definir el cuarteto como la suma de dos tripletes, uno formado por los primero términos de la ec. 63, menos un término denominado cruzado que puede escogerse entre ϕ~ h+~ k,ϕ~ h+~ lyϕ~ k+~ l; y otro segundo triplete formado por los dos últimos términos de la ec. 63 más el término cruzado escogido, T1=ϕ~ h+ϕ~ k−ϕ~ h+~ ky T2=ϕ~ l+ϕ−~ h−~ k−~ l+ϕ~ h+~ k. Si el término |E~ h+~ k|es grande, ambos términos serán cero, y obtendremos Q≃0con una mayor fiabilidad. Por el contrario, también se puede demostrar que si los valores de las magnitudes cruzadas son bajos, Q≃π. Tenemos entonces que las cuatro fases del cuarteto dependen de siete términos, de manera que los estudios de probabilidad se centraron en establecer una teoría que incluyese los tres términos adicionales. Este estudio fue realizado por Hauptman en 1975, quien obtuvo la distribución de probabilidad de Qen función de estos siete factores, estableciendo el principio de ’vecindario’ (Neighborhood Principle), y posteriormente fue refinada y generalizada por Giacovazzo en lo que se conoce como teoría de la representación. 6.2.1. Teoría de la representación En primer lugar, es interesante entender el concepto de ’neighborhood’ introducido Hauptman en su estudio del cuarteto [23]. La idea es que si utilizamos sólo los factores básicos |E~ h|,|E~ k|,|E~ l|y|E~ h+~ k+~ l|obtenemos una estimación de la fase con una fiabilidad menor. Si ahora sumamos a los factores anteriores, los tres factores cruzados, se tiene que la estimación de la fase mejora, siendo la varianza menor que en el caso anterior. 33 Hauptman entonces dedujo que, si añadiamos más factores a lo anterior, obtendríamos un valor de la fase aún mejor. Imaginamos cada conjunto de factores añadido como una capa que va rodeando al anterior, obteniendo una estructura en forma de ’vecindarios anidados’ (nested neighborhoods). Podemos observarlo en la fig. 12, Figura 12: Esquema de capas en la teoría de la representación. En la primera capa que rodea a la fase observamos los factores principales. En la segunda capa se recogen los tres términos cruzados [23] . . Giaccovazzo generalizó este principio, proponiendo su teoría de la representación, que proporciona reglas para identificar las relaciones entre fases haciendo uso de la simetría del cristal. Destaca también la importancia de clasificar las reflexiones en ’sets’ o ’capas’, al igual que Hauptman, aunque ambas clasificaciones no coinciden entre sí. En primer lugar, se define la primera representación de un invariante de estructura Φ, formada por el propio invariante más el conjunto de invariantes que difieren del original en un ángulo constante producto de las operaciones de simetría. Dentro de esta primera representación, se introduce también el concepto de primera capa fásica, compuesta por los factores |E|que intervienen en la primera representación. Como ejemplo, vamos a utilizar el caso del cuarteto. Suponiendo que no hubiese ninguna operación de simetría, la primera representación estaría formada por Q, y la primera capa fásica por las siete magnitudes |E|ya conocidas. En el caso de que si exista una operación de simetría C= (R, ~ T), tal que deja alguna 34 reflexión cruzada invariante, (~ h+~ k)R= (~ h+~ k), tendríamos: Q0=ϕ~ hR +ϕ~ kR +ϕ~ l+ϕ−~ h−~ k−~ l(65) Que sería equivalente a Q,Q0=Q+ 2π(~ h+~ k)~ T, pero con nuevas magnitudes cruzadas |E~ hR+~ l|y|E~ kR+~ l|, teniendo así la primera capa nueve términos y obteniendo una mejor estimación de Q. A continuación introducimos la segunda representación. Para un triplete, está formada por el conjunto de quintetos dados por: C=T+ϕ~ H−ϕ~ H= (ϕ~ h+ϕ~ k+ϕ−~ h−~ k+) + ϕ~ H−ϕ~ H(66) Donde ~ Hes cualquier vector de la red recíproca que se refiera a una reflexión alta. Nos encontramos así que Ces equivalente a T, con la única diferencia de que la primera va a depender aparte de las tres |E|del triplete, de |E~ H|y de las seis magnitudes cruzadas, |E~ h±~ H|,|E~ k±~ H|y|E~ l±~ H|, estando así la segunda capa fásica formada por 10 elementos. Si en vez de seleccionar tan solo un vector ~ H, seleccionamos M, tendremos que la segunda representación estará formada por Mquintetos, y la segunda capa fásica estará formada por 10 ×Melementos. Para el caso de los semi-invariantes de estructura, también podemos definir estas representaciones y capas fásicas. En el caso general, consideramos un semi-invariante con nfases, Φ = ϕ~ h1+ϕ~ h1+... +ϕ~ hn(67) Debemos encontrar una fase que ϕ~ Hy dos operaciones de simetría CiyCjque hagan que ψ1=ϕ~ h01+ϕ~ h01+... +ϕ~ h0n+ϕ~ HRi−ϕ~ HRj(68) sea un invariante de estructura, siendo ~ h0=R~ h. Así, podemos estimar el valor de ψ1y, a partir de este, calcular Φ, al estar relacionadas por una operación de simetría. Tenemos entonces que el conjunto de los invariantes ψ1y las magnitudes básicas y cruzadas asociadas formarán la primera representación y primera capa fásica. 35 Para el caso de la segunda, volvemos a trabajar con el caso del triplete. Suponemos que tenemos la relación entre fases: Φ = ϕ~ h+ϕ~ k+ϕ~ l(69) La segunda representación va a estar formada por el conjunto de quintetos, ψ2=ϕ~ h+ϕ~ k+ϕ~ l+ϕ~ H−ϕ~ H(70) y la segunda capa de fases, por las magnitudes básicas y cruzadas que intervengan en la estimación de los quintetos. La importancia de esta teoría es que, una vez encontramos los invariantes y la jerarquía de estructuras, es posible hallar una distribución de probabilidad para las fases. Esta es una ardua tarea y, actualmente, solo ha sido éxitosa para los casos más simples. 6.3. Procedimiento para la obtención de fases mediante métodos directos. A continuación vamos a describir cualitativamente los pasos que se suelen seguir a la hora de determinar las fases [6]. 1. Una vez tenemos los módulos de los factores de estructura, lo que se hace en primer lugar es normalizarlos, ya que es fundamental para su posterior análisis estadístico. Para esto, es necesario hacer uso de la ec. 79 que encontramos en el anexo B. 2. El segundo paso es buscar invariantes o semi-invariantes de fase. En concreto, nos vamos a centrar en la búsqueda de tripletes, ya que es lo más práctico. Nos interesa por lo tanto fijarnos en las reflexiones ~ hque posean un |E|alto, para posteriormente encontrar todos los pares ~ ky~ h−~ kque cumplan la misma condición de |E|. Se debe tener en cuenta el papel de la simetría en este proceso, que nos facilita la búsqueda. Sea ~ Riel conjunto de matrices de rotación del grupo espacial, a partir de una reflexión ~ h, vamos a buscar las combinaciones con ±~ k~ Ricompletando el triplete con (en el caso de que la reflexiones tengan un valor alto) ~ h±~ k~ Ri. 36 3. Elección de origen y enantioforma: Una vez tenemos los invariantes y semi-invariantes de estructura, debemos hallar las fases ϕen cuestión. El problema es que, mientras que las magnitudes |E|yΦcon las que hemos estado trabajando hasta ahora, son independientess del origen de la celda unidad, y no difieren entre una estructura y su enantioforma, las fases individuales si que dependen de estos dos factores. En el caso de las dos enantioformas, se tiene que para estructura no centrosimétricas, tienen signos opuestos para todas las fases. Es por lo tanto necesario, elegir un origen y una de las enantioformas. Los puntos propuestos, denominados puntos permitidos, para la elección de origen vienen determinados por la simetría del cristal. En el caso de estructuras centrosimétricas, estos puntos coinciden con los centros de inversión. Para distinguir el origen adecuado de entre todos los puntos permitidos, debemos hacer uso de fases previamente estimadas. 4. Elección de un conjunto óptimo de fases iniciales: Este paso es el más crucial en lo que respecta al procedimiento completo. En el estudio del triplete vimos que, teniendo dos valores de fase, podiamos estimar una tercera. Si teníamos varios pares de fases para la misma reflexión ~ h, podíamos calcular la fase haciendo uso de la fórmula de la tangente, ec. 59. El problema es que aún teniendo un gran número de tripletes, aún es necesario un mayor número de fases para comenzar el proceso. Por tanto, debemos obtener valores de otras fases extras. Utilizando la fórmula de la tangente en cadena a partir de los tripletes iniciales, podemos estimar los valores de otras fases. Es por esto que la elección de estos tripletes iniciales es muy importante, ya que si esta estimación inicial no es buena, el resto de fases calculadas a partir de estos van a tener errores muy altos. Una medida de la bondad de la estimación de las reflexiones del triplete venía dada por α~ h. Tenemos que a partir de hcos Φji=I1(Gj) I0(Gj)=D1(Gj)podemos medir el valor estimado de α~ hpara cada reflexión. hα~ hi= r X j=1 GjD1(Gj)(71) Mediante el llamado proceso de convergencia, vamos a seleccionar las reflexiones aptas para pertenecer al conjunto de valores iniciales. Este es un proceso paso a paso en el que se elimina, cada vez, la reflexión con hα~ himenor, así como todas 37 las relaciones de fase en la cual está involucrada. Se calculan otra vez las hα~ hiy se repite el proceso hasta obtener un conjunto inicial que contiene las reflexiones más fiables. Es importante analizar las hα~ hique van a ser eliminadas, ya que si esta interviene en la elección del origen, debe permanecer. 5. Una vez ya tenemos un conjunto de fases iniciales, podemos calcular el resto mediante principalmente los siguientes métodos: Método de adición simbólica. Asignamos a las fases conocidas un valor simbólico a partir de los cuales vamos a definir el resto mediante relaciones probabilísticas. Es útil solamente en el caso de estructuras centrosimétricas. Método de multisolución. La diferencia con el caso anterior es que en vez de símbolos, asociamos valores aproximados a las fases iniciales. Empíricamente se ha demostrado que este procedimiento tiene éxito aún con errores iniciales de 40-50º. El conjunto de fases iniciales estará compuesto por las fases determinadas por el método de convergencia, por las fases que fijan el origen, en el caso de que sea una estructura no centrosimétrica, una fase que fije la enantioforma y, por último, un número límitado de fases. Si tenemos ng reflexiones generales y ns reflexiones especiales, tendremos un total de, n= 4ng ×2ns (72) combinaciones, de las cuales una será la correcta. Sin embargo, hay que aplicar el proceso de determinación de fases a todas ellas. Enfoque aleatorio del problema de las fases. Los dos métodos anteriores se utilizaban sobre todo en los inicios de los métodos directos. Actualmente se utiliza un procedimiento en el cual asignamos valores aleatorios entre 0 y 360◦, que van a ser introducidas en la fórmula de la tangente. Vamos a tener reflexiones con un valor de |E|alto a los cuales se le van a asignar conjuntos de fases aleatorias que posteriormente, utilizando el método de convergencia van a ser refinadas. La popularidad de este método se debe a dos factores. El primero de ellos es que la solución correcta se obtiene rápidamente y sin la necesidad de utilizar un número muy alto de fases aleatorias iniciales. El segundo es que no es necesario definir fases iniciales que fijen el origen y la enantioforma, ya 38 que todas las fases utilizadas son aleatorias. 6. El siguiente paso es calcular el mapa de densidad electrónica a partir de las fases y|E|obtenidos. El problema es que como el número de reflexiones |E|es limitado debido a los procesos de convergencia en los cuales solo utilizamos los óptimos, han de ser sometidos a técnicas de modificación como las estudiadas en los métodos del espacio directo. 7. Figuras de mérito: Las figuras de mérito son funciones que califican la bondad de nuestro modelo. Ya que en los métodos anteriores vamos a obtener varias soluciones (multisolución), estas figuras se utilizan para descartar la mayor parte de las soluciones. Existen varias funciones con esta función, aunque actualmente se utilizan las dos siguientes: Coeficiente de correlación. CC Fom. CC =(h|Eobs||Ecali|−h|Eobs|ih|Ecal|i) (h|E2 obs|i−h|Eobs|i2)1/2(h|E2 cal|i−h|Ecal|i2)1/2(73) |Ecal|es el factor de estructura normalizado calculado a partir las partes positivas del mapa de densidad electrónica. El valor de CC varía entre 0 y 1, y mientras más se acerque al 1, mejor es la estimación. RAT Fom RAT =CCw h|E2 cal|i (74) En este caso, |Ecal|se calcula usando solo el 30 %de las reflexiones con menor valor de |E|.CCw representa la correlación entre |Eobs|yhcos αi, siendo α= ϕN−ϕP.ϕNse refiere a las fases considerando la contribución al factor de estructura de la estructura completa y ϕPcontando solo las contribuciones de las posiciones atómicas conocidas. 39 7. Resumen y conclusiones Desde principios del siglo pasado, que fue cuando nació, la cristalografía ha sido uno de los más importantes campos de estudio de la Física. Es una disciplina que ha tenido un desarrollo paulatino, ya que incluso a día de hoy, un siglo después, el problema de las fases sigue siendo un reto, requiriendo investigaciones con el objetivo de implementar nuevos métodos para su resolución. Una de las principales diferencias con el resto de ramas físicas es que el problema de las fases no tiene aparentemente solución analítica, algo que en principio incomoda a la Física. Sin embargo, a lo largo de este documento hemos visto como una combinación de la estadística, de la tecnología y los ordenadores, y de por supuesto, los propios principios físicos que rigen el comportamiento de los sólidos, es posible encontrar métodos por distintas vías para llegar a la solución correcta del problema de las fases. Así, hemos estudiado por un lado, los métodos del espacio directo, que comenzaron con el estudio de la función de Patterson. El método de Patterson es sumamente importante ya que, además de ser la primera solución del problema de las fases que no fuese por ensayo y error, dió pie a otros métodos también del espacio directo como el método de reemplazo isomorfo o el método de modificación de la densidad electrónica. Por otro lado, también hemos estudiado los métodos directos, que son aquellos que operan en el espacio recíproco. En primer lugar hemos examinado los principios y la teoría que van a regir estos métodos, para así contruir un método con lo pasos necesario para la obtención de las fases. Lo interesante de ambos tipos es que también se pueden combinar, dando lugar a lo que se denomina métodos del espacio dual, usando tanto métodos de modificación de la densidad electrónica como los procesos directos de estimación de fases. Cabe esperar que, de acuerdo con el desarrollo tecnológico, y con nuevas contribuciones teóricas, estos métodos sean mejorados y optimizados en un futuro e incluso se descubran nuevos métodos. Personalmente, no descartaría que se encontrase un solución exacta, ya que en numerosas ocasiones la física ha superado sus propios límites, teniendo que reestructurar teorías completas y cambiando todo lo que creíamos saber. 40 Anexos A. Teorema de convolución La convolución de dos funciones ρ(~r)yg(~r)se define como, C(~u) = ρ(~r)∗g(~r) = Zρ(~r)g(~u −~r)d~r (75) Donde los límites de integración dependen de los limites de integración. Se cumple que, ρ(~r)∗g(~r) = g(~r)∗ρ(~r)(76) El teorema de convolución nos dice que, T[ρ(~r)∗g(~r)] = T(77) Aplicando esto a la función densidad de estados, llegamos a la ec. 17. B. Factor de estructura normalizado El factor de estructura normalizado es una corrección al factor de estructura en la que se suponen todos los átomos puntuales en reposo (es decir, elimina las vibraciones térmicas de los átomos) y en una escala absoluta que facilita su uso y sus aplicaciones. Se calcula como: |E~ h|=|F~ h| √Σ (78) Donde es un factor que depende de la simetría de la red recíproca y, Σ=PN j=1 f2 j. Otra manera de calcular el factor de estructura normalizado a partir de las amplitudes observadas es: |E~ h|= |F~ h|2 obs K0exp(−2Bs2)PN j=1 f2 j!1 2 (79) Donde s= sin θ/λ, K’ es el recíproco del factor de escala y B el factor de temperatura, que modela el efecto del movimiento vibratorio de los átomos debido a su energía. 41