Determinación, mediante simulaciones computacionales, de la secuencia de apilado en estructuras compactas de capas
Determination of close packed polytype sequence by computational simulation
Pablo Serrano Alfaro1,2*,Ernesto Estevez-Rams 1,2
1Facultad de Física de la Universidad de la Habana, San Lazaro y L, CP-10400C. Habana, Cuba. pablo@fisica.uh.cu
2Institute of Materials Research and Engineering (IMRE). estevez@imre.uh.cu ]]>
*Autor para la correspondencia: pablo@fisica.uh.cu
RESUMEN
Se trabaja en la obtención de la secuencia de apilado en estructuras compactas de capas (CPs), a partir de información contenida en el patrón de difracción de rayos x. Este problema solo se había abordado anteriormente de forma analítica, y su solución exacta es intratable para secuencias reales. Con este fin se emplea un algoritmo de búsqueda meta-heurística para problemas de optimización global: simulated annealing. Se tomó como entrada del problema las funciones de correlación de capas, las cuales son extraídas del antes mencionado patrón de difracción en muestras policristalinas. Para reducir la complejidad computacional del problema se tienen en cuenta una serie de aspectos que pasan por áreas como la cristalografía teórica y la computación de alto rendimiento (generación de politipos). Se hace una revisión de magnitudes importantes dentro de este tema, como son, la longitud de correlación y la matriz de correlación. Se han validado los resultados obtenidos mediante secuencias sintéticas generadas en un amplio rango de tipos, y proporción de defectos planares. Se muestran partes de estos resultados, haciendo énfasis en algunos tipos de estructuras compactas de capas ampliamente conocidas como son las: fcc, hcp y dhcp.
Palabras clave: estructuras compactas de capas, simulated annealing, funciones de correlación de capas, longitud de correlación, polytypes
ABSTRACT
We aim at producing the close packed polytype sequence (CPs), based on information contained in the x-ray diffraction pattern. This problem has previously been addressed analytically (Estevez and Azanza, 2014), and its solution is intractable for real sequences. We use optimization algorithm based on meta-heuristic procedures, simulated annealing. The probability correlation functions, which can be extracted from the aforementioned diffraction pattern in polycrystalline samples, was taken as input to the problem. To reduce the computational complexity of the problem a number of consideration are taken into account, from areas such as theoretical crystallography and high-performance computing (polytypes generation). Important quantities are reviewed, as, the correlation length and the correlation matrix. We have validated the results obtained using synthetic sequences generated in a wide range of types, and proportion of planar defects. Partial results are shown, with emphasis on some types of close packed structures widely known such as: fcc, hcp and dhcp.
Key words: close packed, simulated annealing, correlation function, correlation length, polytypes
INTRODUCCIÓN
La obtención de la secuencia de apilado en estructuras compactas de capas (CPs por sus siglas en inglés) resulta importante en varios sentidos. Primeramente, es conocido como obtener el patrón de difracción a partir de la mencionada CPs, lo cual no es la dirección de este trabajo. Conociendo la CPs se pueden realizar otros tipos de análisis en los cuales nuestro grupo de investigación se encuentra trabajando. Entre ellos destacan el estudio de la complejidad estructural vinculada con la complejidad matemática y computacional que subyace en las diferentes codificaciones de estas estructuras.
En la literatura se ha abordado cómo influyen diversos tipos de defectos en los patrones de difracción de rayos x (Berliner y Werner, 1986; Estevez, Martinez, 2001; Estevez, Aragon, 2003; Estevez, Welzel, 2008). Sin embargo, el problema inverso no se encuentra documentado. La primera referencia a dicho problema aparece en Estevez y Azanza, 2014, donde se trata de forma analítica y se llegan a un sistema de ecuaciones que lo resuelven de forma exacta. La dificultad radica en que, desde el punto de vista práctico, dicha solución no es viable por el carácter inexacto de la data. No obstante, el mencionado trabajo tiene valor teórico ya que no solo predice la existencia de la solución, sino que también establece una relación univoca entre las funciones de correlación y los CPs cristalográficamente equivalentes.
La información de entrada en el problema abordado son las funciones de correlación de capas, las cuales se obtienen directamente del patrón de difracción (Estevez, Aragon, 2003). Por lo tanto, el objetivo que nos hemos planteado en este trabajo es la extracción de la CPs, a partir de las funciones de correlación de capas. Para esto utilizamos el algoritmo de búsqueda meta-heurística para problemas de optimización global, simulated annealing.
Estructuras compactas de capas. Notaciones
Es conocido que una estructura compacta consiste en capas hexagonales idénticas de esferas, apiladas en la dirección perpendicular a las mismas (Serrano y Estevez, 2012), de forma tal que, las esferas de una capa descansan en los intersticios de las capas colocadas por encima y por abajo. De esta manera resulta la mayor compacidad posible. Muchos cristales en la naturaleza presentan esta estructura, los casos más sencillos son: cúbica centrada en las caras(fcc), hexagonal compacta(hcp), y la doble hexagonal compacta(dhcp). Debido a ello, estas estructuras pueden ser descritas por las siguientes notaciones:
ABC: Respecto a un cierto origen de coordenadas, las capas que tengan un desplazamiento horizontal nulo se le designa la letra A, las que tengan un desplazamiento horizontal de 
son los vectores base de las capas hexagonales, una letra B, y las que tengan un desplazamiento horizontal de 2una letra C. La condición de compacidad impone que no pueda existir un desplazamiento horizontal nulo entre dos capas consecutivas, en otras palabras, no pueden existir dos letras iguales de forma consecutiva. Más detalles en Estevez, Martinez, 2001; Estevez, Aragon, 2003.
Hägg: Notación binaria en la cual se considera un desplazamiento positivo entre capas el vector , a este desplazamiento se le asigna la letra 1 (ej. AB, BC, CA), en el caso contrario la letra 0 (ej. BA, CB, AC).
Entre estas tres notaciones existe una relación biunívoca.
Funciones de correlación de capas
Las funciones de correlación de capas Ps
dan la probabilidad de encontrar dos capas separadas a una distancia perpendicular 
capas, con un desplazamiento horizontal entre ambas de s veces el vector . Nótese que solo existen tres funciones de correlación. P0 corresponde a la correlación entre capas del mismo tipo, P1 y P2
a la correlación entre capas con desplazamiento positivo y negativo entre ellas respectivamente.
A partir de la anterior definición podemos definir la matriz de correlación (CM) para una CPs de n capas, como una matriz de 3 filas y n columnas. Las tres filas corresponden a las tres funciones de correlación, y las n columnas corresponden a las n posibles separaciones verticales entre capas (desde 0 a n-1). Más detalles en Serrano y Estevez, 2012.
Longitud de correlación
Definamos al politipo como el bloque mínimo de capas a partir del cual se obtiene la estructura compacta periódica (ej. fcc: ABC, hcp: AB, dhcp: ABCB). Para estructuras compactas perfectas obtenidas mediante la repetición periódica de un politipo, las funciones de correlación son periódicas y tienen como período la longitud del politipo (Figura. 1). Más detalles en Serrano y Estevez, 2012.
En la mayor parte de los casos las estructuras cristalinas crecen con defectos, entiéndase por ello cualquier alteración en el apilado perfecto de capas, que pueden ser de diferentes tipos y naturaleza (Estevez, Welzel, 2008). Cuando esto sucede decimos que se pierde correlación, las funciones de correlación dejan de oscilar y tienden a valores determinados por las frecuencias de los tipos de capas A(Pa), B(Pb) y C(Pc) en el apilado, dados por las siguientes ecuaciones. Ver Figura. 2.
]]>
En estas ecuaciones c es la longitud de correlación. Esta magnitud es importe en los posteriores cálculos y análisis que se realizarán, y por lo tanto requiere de algunos comentarios acerca de su interpretación. Matemáticamente se define por Estevez y Welzel, 2008, como la longitud máxima para la cual las envolventes superiores e inferiores de las funciones de correlación colapsan en un mismo valor. En la práctica la solemos determinar estableciendo una ventana de ancho mucho menor que 1(amplitud máxima de la oscilación de las funciones de correlación) centrada en los valores dados por las Ec. 1, y determinando la longitud máxima para la cual las oscilaciones no sobrepasan los límites inferiores y superiores de dicha ventana. Entonces el significado que le damos es que, a partir de dicha longitud dentro de nuestro CPs, no somos capaces de predecir que capas nos encontraremos con un grado de certeza mayor que el conocimiento acerca de la proporción de capas que existe dentro de la secuencia. En otras palabras, es como si nos situáramos en una posición aleatoria dentro de una secuencia perfecta conocida, a la cual se le han introducido una cantidad de defectos en posiciones aleatorias, y solo tenemos la información de la capa que existe en el justo lugar donde nos posicionamos, entonces nos hacen la siguiente pregunta, ¿Qué capa encontrarás? capas antes o después? La respuesta que daríamos es la siguiente:
<< c, con casi total certeza podemos predecir que capa encontraremos, porque conocemos la secuencia base perfecta.<c, aún somos capaces de predecir con mayor certeza de la que brindan las proporciones de tipos de capas, aunque vamos perdiendo esta capacidad de predecir a medida que nos acercamos a .>c , solo podemos decir que con una probabilidad Pa encontraremos una capa A, y así análogamente para las capas B y C. Es decir, nuestra capacidad de predecir se reduce al mínimo dado nuestro conocimiento previo acerca de las proporciones.
SIMULACIÓN Y RESULTADOS
Dentro de la familia de métodos computacionales de búsqueda meta-heurística para problemas de optimización global, uno de los más utilizados en el campo de la Física es simulated annealing. Su nombre proviene del campo de la industria metalúrgica, y está relacionado con el proceso de recocido del acero y cerámicas, una técnica que consiste en calentar hasta una temperatura alta y luego enfriar lentamente el material para variar sus propiedades físicas. El método fue descrito independientemente por Scott Kirkpatrick, C. Daniel Gelatt y Mario P. Vecchi en 1983, y por Vlado Černý en 1985 y es una adaptación del algoritmo Metropolis-Hastings dentro del método de Montecarlo, utilizado en sus inicios para la minimización de la energía en los posibles estados dentro de sistemas termodinámicos. Su esquema general es el siguiente:
Se selecciona una propuesta de solución inicial Xi (CPs), una seudo-temperatura inicial y otra mínima para el sistema, y un esquema de enfriamiento. En nuestro caso además se capta como entrada 1, que no es más que la matriz de correlación de entrada.
Se escoge una función de costo F (CM, CM1); de tal forma que, mientras más pequeño sea su valor, más coincidencia hay con las funciones de correlación de entrada.
Se evalúa CM(Xi).
Se perturba Xi para obtener una propuesta vecina Xi+1.
Se evalúa CM(Xi+1)
Se evalúa 
Si E < 0 acepto Xi+1 como la nueva solución actualizada, en el otro caso acepto Xi+1 con una probabilidad 
Reduzco la seudo-temperatura del sistema.
Si la temperatura actual es mayor que la temperatura mínima regreso al punto 3, en el caso contrario finalizo la simulación y obtengo como salida la mejor solución encontrada.
Este algoritmo garantiza hallar siempre el mínimo global, si el sistema es ergódico y el tiempo de búsqueda es lo suficientemente grande. El tiempo de búsqueda está relacionado con el esquema de enfriamiento del sistema. El hecho de aceptar en ocasiones soluciones peores que las que se tenían previamente evita quedar atrapado en mínimos locales.
Con la Ecuación 1 podemos determinar el comportamiento de las funciones de correlación para >c, sin embrago, dado que nuestro problema es el inverso, estas ecuaciones no brindan una solución única para la determinación de las frecuencias de los tipos de capa. El problema de conocer de antemano las frecuencias de los tipos de capa es esencial a la hora de proponer una solución inicial en la simulación, ya que garantiza el cálculo de las funciones de correlación solo para valores menores que c. Esto hecho es fundamental para reducir la complejidad algorítmica en el cálculo de CM, que pasa de ser O(n2) a O(n).
En las simulaciones se han tenido en cuenta otros aspectos, siempre con el objetivo de disminuir el tiempo de cálculo. Entre ellos destacan la generación eficiente de politipos cristalográficamente no equivalentes (Estevez, Mojicar, 2008), que constituye una adaptación a los algoritmos de Cattel y Ruskey, 2000. Las siguientes figuras (4 y 5)muestran ejemplos de simulaciones en estructuras compactas.
]]>CONCLUSIONES
Se han estudiado hasta la fecha una serie de casos, todos generados sintéticamente y procurando ser exhaustivos en cuanto a la utilización de diferentes tipos de errores y la proporción con que estos son generados. En todos casos hemos obtenido buenos resultados en las simulaciones, aunque la práctica nos ha demostrado que es preciso ser cuidadoso en la selección de la propuesta inicial de la simulación. Esto es debido a que el espacio de búsqueda o de configuraciones de nuestro problema es muy extenso y con las dinámicas de perturbación de las soluciones que estamos aplicando actualmente no somos capaces de recorrer dicho espacio, haciendo que nuestro sistema no sea ergódico. Por lo tanto, una parte importante de la simulación radica en le selección de la propuesta inicial, la cual después es mejorada con el algoritmo de búsqueda.
REFERENCIAS BIBLIOGRÁFICAS
ESTEVEZ, E and AZANZA, C. On the non-recursive determination of close packed polytype sequence from the interlayer displacement count (2014).
R. BERLINER and S.A. WERNER, PHYS. Rev. B 34, 3586 (1986).
ESTEVEZ, E. MARTINEZ, A. PENTÓN, A and LORA, R. PHYSICAL REVIEW B, VOLUME 63, 054109 (2001).
ESTEVEZ, E. ARAGON, B. FUESS, H and PENTÓN, A. PHYSICAL REVIEW B 68, 064111 (2003).
ESTEVEZ, E. WELZEL, U. PENTÓN, A and MITTEMEIJER, E. Acta Cryst. A64, 537–548 (2008).
KIRKPATRICK, S. GELATT, C. D. VECCHI, M. P. (1983). Optimization by Simulated Annealing. Science (en inglés) 220 (4598): 671–680
SERRANO, P. ESTEVEZ, E. Tesis de Lic. Estudio de polytypes compactos por magnitudes derivadas de la entropía de Shannon (2012).
RODRIGUEZ, E. ESTEVEZ, E. Tesis de Lic. Estudio mediante simulaciones de Monte Carlo del diagrama de fases y la cinética de transformación en estructuras politípicas (2014).
CATTEL, K., RUSKEY, F., SAWADA, J., SERRA, M. & MIERS, C. R. (2000). J. Algorithms, 37, 267–282.
]]> ESTEVEZ, E. MOJICAR, J. Acta Cryst. (2008). A64, 529–536
Recibido: 15/06/2016
Aceptado: 10/10/2016