Tesis:
Monte Carlo simulations of crystal polymorphism in polymer-based systems
- Autor: HERRANZ FEITO, Miguel
- Título: Monte Carlo simulations of crystal polymorphism in polymer-based systems
- Fecha: 2023
- Materia:
- Escuela: E.T.S. DE INGENIEROS INDUSTRIALES
- Departamentos: INGENIERIA QUIMICA INDUSTRIAL Y DEL MEDIO AMBIENTE
- Acceso electrónico: https://oa.upm.es/73157/
- Director/a 1º: KARAGIANNIS, Nikolaos
- Resumen: Polymers are encountered in many aspects of our lives, varying from the proteins or nucleic acids of our cells to the synthetic ones that we use for packaging, construction, storage, clothing, communication, transportation among an endless list of materials, objects, and applications. Their outstanding properties render them essential for cutting-edge technologies as batteries, solar panels, smart materials, fuel cells, or drug delivery vessels. Even though their relevance in modern life is undisputed, there are still a lot of unresolved questions around them and open challenges in engineering novel polymer-based materials with enhanced characteristics.
Many of the polymer properties can be understood by gauging the behavior at the level of atoms and molecules. The way monomers and atoms are connected along the chain backbone along with the interactions between the difference species, can provide an explanation for the distinctly different macroscopic properties. However, studying such complex systems to establish a connection between behavior of atoms and macroscopic properties is far from trivial requiring advanced experimental or simulation techniques.
In many-body systems like the aforementioned ones based on polymers, the large population of atoms and molecules requires a description based on statistical physics and, as a consequence, macroscopic laws are necessarily of a statistical nature. Molecular simulations on macromolecular systems are based on this statistical description where macroscopic observables are obtained as averages over microscopic properties. In general, two main branches of molecular simulation exist: Molecular Dynamics (MD) and Monte Carlo (MC). While MD is necessary when time dependence and dynamics are important, it is computationally very expensive when applied on slow-evolving systems like the macromolecular ones. On the other hand, Monte Carlo algorithms, by being stochastic, lack dynamical information but can be very efficient for the rapid equilibration and robust sampling of the phase space.
Compared to experiments, molecular simulations allow the study of “what-if” questions at extreme conditions with minimal effort compared to normal ones, are cost-effective and pose no safety or environmental hazards. Naturally, the core of the simulation is based on a model, a collection of mathematical equations that represent the behaviour of a system or process. This model, is regularly, a result of a trade-off between the desired accuracy and the required computational time. Very complicated models can be accurate but also computationally intensive, while coarse-grained ones lack accuracy but can access longer time and length scales.
The current state of molecular simulations clearly does not allow them to replace the conventional experimental techniques in most of the existing applications. Rather, they act as a valuable complementary tool, that often acts a bridge that unites theoretical predictions with experimental observables. In specific cases it can also lead the design of hybrid materials and minimize the number of samples required to be synthesized, characterized and eventually disposed.
The present thesis focuses on one important aspect of polymers, which is self-organization, and in particular their phase behavior and crystallization. Crystallization in real polymers is usually associated with the formation of ordered regions, called lamellae, where chains fold and align. Process conditions, chain stiffness, confinement and or presence of nanofillers affect crystallization, usually leading to semi-crystalline polymers where lamellae co-exist with disordered regions. However, the current study focuses on simplified models of polymer representation, which are closer to (hard) colloidal systems than “realistic” or “traditional” ones.
In all physical systems to be presented here macromolecules are represented as freely-jointed, linear sequences of non-overlapping spheres. Bond lengths are considered as either tangent, within a numerical accuracy, or are characterized by the presence of gaps of varied width. All intra- and inter-molecular interactions are described by the hard sphere (HS) or the square well (SW) potential. The former corresponds to an athermal system, where phase behavior is dictated by an increase in the total entropy of the system, while the latter corresponds to short-range attraction and self-organization is dominated by the total free energy.
All simulation studies to be reported are carried out using the Simu-D simulator-descriptor. The software, which was developed in the course of the present thesis, consists of two main parts. First, the main part (simulator) is built around a new Monte Carlo suite for polymer-based systems under various conditions, including in the bulk, under confinement, on surfaces/interfaces and in nanocomposites. Second, the analysis part (descriptor) is based on the latest version of the Characteristic Crystallographic Element (CCE) norm, a refined metric to characterize local structure by comparing the local arrangement of monomers against specific reference crystals in two or three dimensions. This descriptor tool is essential for the precise detection of phase transitions, the characterization of the established morphologies, and to accurately distinguish between different polymorphs.
From the technical perspective this modelling and simulation tool aims to be a general, efficient, and user-friendly software for a wide variety of coarse-grained, polymer-based systems. Present applications of the tool include bulk polymers at different packing densities and chain lengths, nanocomposites with nanofillers in the form of cylinders and spheres of varied concentration and size, extremely confined and maximally packed assemblies in two and three dimensions, blends of monomers and polymers and macromolecules terminally grafted on surfaces and/or nanoparticles. The Monte Carlo protocol is built around diverse types of Monte Carlo algorithms especially designed for polymers, mainly conventional local moves (like reptation, end-mer rotation, flip etc) and more advanced chain-connectivity-altering moves. Furthermore, new group of moves, as identity-exchange or cluster-based ones, are incorporated here to fulfil the modelling needs of different systems. Extensive scripting supports the use of both parts of the Simu-D, to be tailored for the demands of the simulated system. For example, the mix of moves, including their types, attempt probabilities and number of trial configurations (for local moves) is determined automatically by comparing the present parameters (chain length, packing density, presence of confinement and/or nanoparticles) with optimized ones, existing in data libraries.
All simulations of the work follow the same methodology: it is comprised of three major steps. First, an initial configuration is generated, subject to the constraints imposed by the physical system to be studied, including the number and type of different species (chain monomers, free monomers (solvent), nanoparticles), intra- (molecular/chain architecture) and inter-connectivity (anchoring) of the species, packing density, level of confinement and its geometry, among other factors. Then, the Monte Carlo scheme is conducted on the order of hundreds to billions or even trillions of iterations (steps), first to reach the equilibrium state and then to sample representative equilibrium microstates, each one being as structurally uncorrelated from the others as possible. Usually, this main part is performed under isochoric conditions, although the Simu-D suite can also perform isobaric-isothermal simulations.
The last step is the post-simulation analysis of the data as collected in trajectories through very long MC simulations. This step implies conventional post-processing to quantify the short- and long-range structure of the polymer chains. More importantly, it entails the use of the CCE norm descriptor to precisely identify the local environment and to quantify its similarity to reference crystals or local symmetries. Further analysis includes the quantitative characterization of the Voronoi cell around each atom or particle of the system and how this is affected by the observed phase transition.
The simulation approach, as described above, is applied on two main and quite different physical systems. In the first case we study polymorphism and perfection in the athermal crystallization of freely jointed chains of hard spheres at packing densities above the melting transition. In the second case we study self-assembly of crystal polymorphs of attractive chains whose monomers interact with the square well potential, under dilute conditions.
The crystal nucleation and growth of definitely entangled polymers of hard spheres is studied using Simu-D in an unprecedentedly long MC simulation, spanning several years. The simulated systems comprise of 54 chains of average length 1000 monomers in the bulk, under periodic boundary conditions. The Hard Sphere (HS) model is the simplest model that incorporates the concept of excluded volume and treats monomers as impenetrable spheres without any short- or long-range interactions. Starting from a purely amorphous configuration, and under constant volume the Monte Carlo scheme drives the system to a final crystal of very high perfection. In the whole transition, a competition occurs between the hexagonal closed packed (HCP) and the face centered cubic (FCC) close packed crystals. In the past, in the majority of modelling studies, be on monomeric or polymeric hard spheres, established ordered structures corresponded to a defect-ridden, random hexagonal close packed (rHCP) polymorph, of dual HCP and FCC character, with a single (chains) or multiple (monomers) stacking directions. Here, we demonstrate that given enough computational time and a robust equilibration scheme, the polymeric system of long and entangled chains can grow towards a perfect FCC crystal, passing through various polymorphs in full accordance with the Ostwald´s rule of stages. The analysis of the conformational characteristics of the chains reveals a loss in the conformational entropy during transition to the crystal phase. However, it is more than compensated by an increase in the translational entropy, as revealed from the corresponding analysis of the spatial arrangement and accessible volume.
The second case study corresponds to the phase behavior of freely jointed chains of tangent spheres interacting with the square well (SW) potential. This interaction model represents monomers again as impenetrable spheres but now their pair-wise energetic (attractive) interactions is defined by two parameters: intensity, ε, and range, σ2 of attraction. The square well is perhaps the simplest model to conduct simulations on free-energy-driven phenomena and processes at the atomistic level. Starting from athermal freely jointed chains at very dilute conditions, the SW potential is activated by fixing the values of intensity and range of interaction. This provokes an immediate quench of polymer chains, which coil and get together forming cluster(s) due to the attractive interactions between their monomers. These clusters can eventually crystallize depending strongly on the applied intensity and range of interaction. We embark on a vast collection of simulations with different values of intensity and range of interaction to explore the complete phase behavior and the corresponding established cluster morphologies. On one hand, the intensity of interaction is found to be the key parameter that determines the obtention of amorphous, crystal of glassy domains. The glassy clusters are rich in fivefold symmetry sites, acting as crystal inhibitors, and effectively trap the polymers in a structure that possesses short- but lacks long-range order. On the other hand, the range of interaction is the parameter that determines the type of crystal polymorph we can obtain within the crystal domain. A simple geometrical model is proposed which can predict the theoretically stable morphologies ignoring entropic contributions as the ones imposed by chain connectivity. Still, the latter is incorporated by considering that tangency of bonded spheres leads to corresponding crystals with no gaps between successive sites. In two dimensions there is an excellent agreement between the theoretical predictions and simulation results on the appearance of triangular and square crystals and on the onset of the corresponding transitions between the dominant polymorphs as a function of attraction range. In three dimensions the results of the simulations in Simu-D reveal the appearance of close packed structures as FCC and HCP, but also non-compact crystals as BCC and HEX. Furthermore, simulations reveal the existence of Frank-Casper quasi-crystals, whose geometry does indeed place them as favourable structures in the SW phase diagram of chains. These ordered morphologies are further accompanied by two regions of clearly amorphous character, which are denoted here as “dead zones”. The present findings could be considered as a first step towards a more complicated computer-aided design of polymer crystals, based on hard colloids, through proper tuning of the range and intensity of their interactions, in a fashion similar to the pioneering concept of digital alchemy.
The present thesis is comprised by the published articles of the aforementioned works along with an overview of the theoretical framework, the employed methodology, an analysis of results not included in the publications, and a list of major conclusions and current future lines of research that stem from the present study.
RESUMEN
Los polímeros se encuentran en muchos aspectos de nuestras vidas, desde las proteínas o ácidos nucleicos de nuestras células hasta los sintéticos que usamos para embalajes, construcción, almacenamiento, ropa, comunicación o transporte, entre una lista interminable de materiales, objetos y aplicaciones. Sus sobresalientes propiedades los hacen esenciales para tecnologías de vanguardia como baterías, paneles solares, materiales inteligentes, celdas de combustible o encapsulamiento de medicamentos. Aunque su relevancia en la vida moderna es indiscutible, todavía hay muchas preguntas sin resolver a su alrededor y desafíos abiertos en la ingeniería de nuevos materiales basados en polímeros con características mejoradas.
Muchas de las propiedades de los polímeros se pueden entender midiendo el comportamiento a nivel de átomos y moléculas. La forma en que los monómeros y los átomos están conectados a lo largo de la columna vertebral de la cadena, junto con las interacciones entre las diferentes especies, puede proporcionar una explicación de las propiedades macroscópicas claramente diferentes. Sin embargo, estudiar sistemas tan complejos para establecer una conexión entre el comportamiento de los átomos y las propiedades macroscópicas está lejos de ser trivial y requiere técnicas experimentales o de simulación avanzadas.
En sistemas de múltiples cuerpos como los antes mencionados, la gran población de átomos y moléculas requiere una descripción basada en la física estadística y, en consecuencia, las leyes macroscópicas son de naturaleza estadística. Las simulaciones moleculares de sistemas macromoleculares se basan en esta descripción estadística donde los observables macroscópicos se obtienen como promedios sobre las propiedades microscópicas. En general, existen dos ramas principales de simulación molecular: Molecular Dynamics (MD) y Monte Carlo (MC). Si bien la MD es necesaria cuando la dependencia del tiempo y la dinámica son importantes, es computacionalmente muy costosa cuando se aplica en sistemas de evolución lenta como los macromoleculares. Por otro lado, los algoritmos de Monte Carlo, al ser estocásticos, carecen de información dinámica, pero pueden ser muy eficientes para el equilibrado rápido y muestreo robusto.
En comparación con los experimentos, las simulaciones moleculares permiten el estudio de preguntas-hipotéticas en condiciones extremas con un esfuerzo mínimo, son rentables económicamente y no plantean peligros para la seguridad o el medio ambiente. Naturalmente, el núcleo de la simulación se basa en un modelo, una colección de ecuaciones matemáticas que representan el comportamiento de un sistema o proceso. Este modelo, es regularmente, el resultado de un compromiso entre la precisión deseada y el tiempo computacional requerido. Los modelos muy complicados pueden ser precisos, pero también computacionalmente intensivos; mientras que los de grano grueso carecen de precisión, pero pueden acceder a escalas de tiempo y longitud más largas.
El estado actual de las simulaciones moleculares claramente no les permite reemplazar las técnicas experimentales convencionales. Más bien, actúan como una valiosa herramienta complementaria, que a menudo actúa como un puente que une las predicciones teóricas con los observables experimentales.
La presente tesis se centra en un aspecto importante de los polímeros, que es la autoorganización y, en particular, el comportamiento de fase y la cristalización. La cristalización en polímeros reales generalmente se asocia con la formación de regiones ordenadas, llamadas lamelas, donde las cadenas se pliegan y se alinean. Las condiciones del proceso, la rigidez de la cadena, el confinamiento o la presencia de nanorrellenos afectan la cristalización, lo que generalmente conduce a polímeros semicristalinos donde las lamelas, coexisten con regiones desordenadas. Sin embargo, el estudio actual se centra en un modelo simplificado de representación de polímeros, que está más cerca de los sistemas coloidales que de los "realistas".
En todos los sistemas físicos que se presentarán en este trabajo, las macromoléculas se representan como secuencias lineales de esferas que no se superponen y unidas libremente. Las longitudes de enlace se consideran tangentes, dentro de una precisión numérica, o se caracterizan por la presencia de tolerancias de enlace. Todas las interacciones intermoleculares e intermoleculares serán descritas por la esfera dura (HS) o el potencial de pozo cuadrado (SW). El primero corresponde a un sistema atérmico, donde el comportamiento de fase está dictado por un aumento en la entropía total del sistema, mientras que el segundo corresponde a la atracción y la autoorganización está dominada por un aumento en la energía libre total.
Todos los sistemas estudiados se realizan utilizando el simulador-descriptor Simu-D. El software, que se desarrolló en el curso de la presente tesis, consta de dos partes principales. Primero, una nueva suite Monte Carlo para sistemas basados en polímeros en diversas condiciones, incluido en masa, en confinamiento, en superficies/interfaces y en nanocompuestos. La segunda parte es la norma del elemento característico cristalográfico (CCE), una métrica refinada para caracterizar la estructura local al comparar la disposición local de los monómeros con cristales de referencia específicos en dos o tres dimensiones. Esta herramienta descriptora es esencial para la detección precisa de transiciones de fase, la caracterización de las morfologías establecidas y para distinguir diferentes polimorfismos con precisión.
Desde el punto de vista técnico, esta herramienta de modelado y simulación pretende ser un software general, eficiente y fácil de usar para una amplia variedad de sistemas basados en polímeros de grano grueso. Las aplicaciones actuales de la herramienta incluyen polímeros en masa con diferentes densidades de empaquetamiento y longitudes de cadena, nanocompuestos con nanorrellenos en forma de cilindros y esferas de concentración y tamaño variados, ensamblajes extremadamente confinados en dos y tres dimensiones, mezclas de monómeros y polímeros y macromoléculas injertadas terminalmente sobre superficies o sobre nanopartículas. El protocolo Monte Carlo se basa en diversos tipos de algoritmos Monte Carlo especialmente diseñados para polímeros, principalmente movimientos locales convencionales (como reptación, rotación final, volteo, etc.) y movimientos más avanzados que alteran la conectividad de la cadena. Además, aquí se incorporan nuevos grupos de movimientos, como el intercambio de identidad o basados en clústeres, para satisfacer las necesidades de modelado de diferentes sistemas. La utilización de scripts permite el uso de ambas partes de Simu-D, adaptándose al sistema simulado. Por ejemplo, la combinación de movimientos, incluidos sus tipos, probabilidades de intento y número de configuraciones de prueba (para movimientos locales) se determina automáticamente al comparar los parámetros presentes (longitud de la cadena, densidad de empaquetamiento, presencia de confinamiento y/o nanopartículas) con los optimizados, existentes en las bibliotecas de datos.
Todas las simulaciones del trabajo siguen la misma metodología que consta de tres pasos principales. Primero, se genera una configuración inicial, sujeta a las restricciones impuestas por el sistema físico a estudiar, incluyendo el número y tipo de diferentes especies (monómeros de cadena, monómeros libres (disolvente), nanopartículas), intra- (arquitectura molecular/cadena) e interconectividad (anclaje) de las especies, densidad de empaquetamiento, nivel de confinamiento y su geometría, entre otros factores.
En segundo lugar, el esquema de Monte Carlo lleva a cabo del orden de cientos a miles de millones de iteraciones (pasos), primero para alcanzar el estado de equilibrio y luego para muestrear microestados de equilibrio representativos, cada uno de los cuales está estructuralmente tan poco correlacionado con los demás como sea posible. Normalmente, esta parte se realiza en condiciones isócoras, aunque la suite Simu-D también puede realizar simulaciones isobáricas-isotérmicas.
El último paso es el análisis de los datos recopilados en las largas trayectorias de MC. Este paso implica un procesamiento posterior convencional para cuantificar la estructura de corto y largo alcance de las cadenas poliméricas. Pero lo que es más importante, implica el uso de la norma CCE para identificar con precisión el entorno local y cuantificar su similitud con los cristales de referencia o las simetrías locales.
El enfoque de simulación, como se describe anteriormente, se aplicará en dos sistemas físicos principales y bastante diferentes. En el primer caso estudiaremos el polimorfismo y la perfección en la cristalización atérmica de cadenas de esferas duras tangentes libremente unidas a densidades de empaquetamiento por encima de la transición de fusión. En el segundo caso estudiamos el autoensamblaje y cristalización de cadenas atractivas cuyos monómeros interactúan con el potencial del pozo cuadrado, en condiciones diluidas.
La nucleación y crecimiento de cristales en polímeros entrelazados de esferas duras se estudia utilizando Simu-D en una simulación MC sin precedentes, que abarca varios años. Los sistemas simulados se componen de 54 cadenas de longitud media de 1000 monómeros en masa, en condiciones de contorno periódicas. El modelo Hard Sphere (HS) es el modelo más simple que incorpora el concepto de volumen excluido y simplemente considera a los monómeros como esferas impenetrables sin interacciones de corto o largo alcance. Partiendo de una configuración puramente amorfa, y a volumen constante, el esquema de Monte Carlo conduce el sistema hacia un cristal final de altísima perfección. En toda la transición, se produce una competencia entre los cristales compactos hexagonales (HCP) y los cristales compactos cúbicos centrados en las caras (FCC). En el pasado, en la mayoría de los estudios de modelado, ya sea sobre esferas duras monoméricas o poliméricas, las estructuras ordenadas establecidas corresponden a un polimorfo hexagonal aleatorio (rHCP) lleno de defectos, de carácter dual HCP y FCC, con una sola o varias direcciones de apilamiento. Aquí, demostramos que, dado el tiempo computacional suficiente y un esquema de equilibrio robusto, el sistema polimérico de cadenas largas y entrelazadas puede crecer hacia un cristal FCC perfecto, pasando por varios polimorfos.
El análisis de las características conformacionales de las cadenas revela una pérdida en la entropía conformacional en la transición a la fase cristalina. Sin embargo, es más que compensado por un aumento en la entropía traslacional, como lo revela el análisis de la disposición espacial y el volumen accesible.
El segundo caso de estudio corresponde al comportamiento de fase de cadenas de esferas tangentes libremente articuladas que interactúan con el potencial del pozo cuadrado (SW). Este modelo de interacción representa nuevamente a los monómeros como esferas impenetrables, pero ahora su interacción energética, en forma de atracción, está definida por dos parámetros: intensidad, ε, y rango, σ2 de atracción. El pozo cuadrado es quizás el modelo más simple para realizar simulaciones de fenómenos y procesos impulsados por energía libre a nivel atomístico. Partiendo de cadenas atérmicas libremente articuladas en condiciones muy diluidas, el potencial SW se activa fijando los valores de intensidad y rango de interacción. Esto provoca un enfriamiento inmediato de las cadenas de polímeros, que se unen debido a la atracción formando grupos. Estos grupos eventualmente pueden cristalizar dependiendo fuertemente de la intensidad aplicada y el rango de interacción. Nos embarcamos en una amplia colección de simulaciones con diferentes valores de intensidad y rango de interacción para explorar el comportamiento de fase completo y las correspondientes morfologías de las agrupaciones establecidas. Por un lado, la intensidad de la interacción resulta ser el parámetro clave que determina la obtención de dominios cristalinos o vítreos amorfos. Los grupos vítreos son ricos en sitios de simetría quíntuple, que actúan como inhibidores de la cristalización y atrapan eficazmente los polímeros en una estructura que posee un orden de corto alcance, pero que carece de orden de largo alcance. Por otro lado, el rango de interacción es el parámetro que determina el tipo de cristal que podemos obtener dentro del dominio cristalino. Se propone un modelo geométrico simple que puede predecir las morfologías teóricamente estables cuando se desprecian los efectos de la entropía, incluida la conectividad de la cadena. En dos dimensiones existe una excelente concordancia entre las predicciones teóricas y los resultados de la simulación sobre la aparición de cristales triangulares y cuadrados y sobre el inicio de las transiciones correspondientes. En 3-D, los resultados de las simulaciones en Simu-D revelan la aparición de estructuras compactas como FCC y HCP, pero también cristales no compactos como BCC y HEX. Además, las simulaciones revelan la existencia de cuasi-cristales de Frank-Casper, cuya geometría los ubica como estructuras favorables en el diagrama de fase. Estas morfologías ordenadas se acompañan además de dos regiones de carácter claramente amorfo, que se denominan “zonas muertas”. Los presentes resultados motivan al diseño de cristales poliméricos mediante un ajuste efectivo del potencial de interacción.
La presente tesis está compuesta por los artículos publicados de los trabajos mencionados junto con una descripción general del marco teórico, la metodología empleada, un análisis de los resultados no incluidos en las publicaciones y una lista de las principales conclusiones y futuras líneas de investigación que podrían derivarse del presente estudio.