Master thesis: Computational study of biodegradable polymers with applications in 3D printing

Estudio computacional de pol´ımeros biodegradables con aplicaci´on en impresi´on 3D Sheintly Karla D´ıaz Gonz´alez Tutora: Petra Baˇcov´a Co-tutor: Sergio I. Molina Trabajo Fin de M´aster en Nanociencia y Tecnolog´ıa de Materiales Facultad de Ciencias, Puerto Real (C´adiz), curso 2024/2025 Resumen Este trabajo presenta un estudio computacional detallado sobre el ´acido polil´actico (PLA), un pol´ımero biodegradable que ha cobrado gran relevancia en ´areas como la impresi´on 3D y el envasado sostenible. Gracias a su origen renovable y su versatilidad en el procesamiento, el PLA se posiciona como un material fundamental en la b´usqueda de soluciones m´as sostenibles. Sin embargo, muchas de sus propiedades t´ermicas y mec´anicas dependen en gran medida de caracter´ısticas estructurales como el peso molecular y la estereoqu´ımica de sus cadenas, lo que resalta la importancia de entender a fondo la relaci´on entre la estructura y las propiedades. Para abordar esta complejidad, en este estudio se ha realizado un an´alisis sistem´atico basado en simulaciones de din´amica molecular, enfocado en examinar la temperatura de transici´on v´ıtrea (Tg), una propiedad crucial en el comportamiento de los pol´ımeros. Para ello, se eligi´o un campo de fuerza adaptado espec´ıficamente para representar con precisi´on homopol´ımeros y copol´ımeros de PLA, lo que permiti´o investigar la influencia del peso molecular y la composici´on estereoqu´ımica. Se modelaron seis sistemas diferentes y se llevaron a cabo procesos de enfriamiento a diversas velocidades, aplicando dos m´etodos de an´alisis del Tgy extrapolando los resultados para compararlos con datos experimentales recientes. Los resultados obtenidos muestran que hay diferencias en el comportamiento t´ermico dependiendo del tipo de sistema. En el caso de los homopol´ımeros tipo L, el Tgdepende del tama˜no de las cadenas de acuerdo con la predicci´on te´orica. En el caso de los copol´ımeros y el homopol´ımero tipo D, los resultados se alinean razonablemente bien con la ecuaci´on Flory-Fox, es decir, muestran que los sistemas que contienen tambi´en el mon´omero D tienen Tgm´as baja que los homopol´ımeros tipo L. Este enfoque ha permitido validar la metodolog´ıa utilizada y acercar de manera m´as precisa los resultados te´oricos a los valores experimentales, lo que refuerza la importancia de los estudios computacionales para investigar y predecir el comportamiento de materiales complejos. En resumen, este trabajo no solo enriquece nuestro entendimiento fundamental del PLA, sino que tambi´en ofrece herramientas esenciales para el dise˜no racional de pol´ımeros biodegradables con propiedades optimizadas para su uso en la industria. ii Abstract This work presents a detailed computational study on polylactic acid (PLA), a biodegradable polymer that has gained significant relevance in areas such as 3D printing and sustainable packaging. Thanks to its renewable origin and processing versatility, PLA is positioned as a fundamental material in the search for more sustainable solutions. However, many of its thermal and mechanical properties largely depend on structural characteristics like the molecular weight and stereochemistry of its chains, which highlights the importance of thoroughly understanding the structure-property relationship. To address this complexity, this study conducted a systematic analysis based on molecular dynamics simulations, focused on examining the glass transition temperature (Tg), a crucial property in polymer behavior. For this, a force field specifically adapted to accurately represent PLA homopolymers and copolymers was chosen, allowing for the investigation of the influence of molecular weight and stereochemical composition. Six different systems were modeled, and cooling processes were performed at various rates, applying two methods for Tganalysis and extrapolating the results for comparison with recent experimental data. The obtained results show that there are differences in thermal behavior depending on the system type. For L-type homopolymers, the Tgdepends on chain size in accordance with theoretical predictions. For copolymers and D-type homopolymers, the results align reasonably well with the Flory-Fox equation, indicating that systems, which also contain the D monomer, have a lower Tgthan L-type homopolymers. This approach has allowed for the validation of the used methodology and a more precise alignment of theoretical results with experimental values, which reinforces the importance of computational studies for investigating and predicting the behavior of complex materials. In summary, this work not only enriches our fundamental understanding of PLA but also offers essential tools for the rational design of biodegradable polymers with optimized properties for industrial use. iv ´ Indice general 1. Introducci´on 1 1.1. Materialespolim´ericos ................................. 1 1.1.1. Estructura ................................... 1 1.1.2. Propiedades................................... 2 1.2. Pr´acticas y pol´ımeros sostenibles . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1. ´ AcidoPolil´actico(PLA)............................ 4 1.2.2. Impresi´on3D.................................. 5 1.2.3. T´ecnicas computacionales en la ciencia de los pol´ımeros . . . . . . . . . . 6 1.2.4. Simulaci´on atom´ıstica de PLA: Estado del arte . . . . . . . . . . . . . . . 7 2. Objetivos 9 3. Materiales y M´etodos 11 3.1. Din´amicaMolecular .................................. 11 3.1.1. Ecuaciones de movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.2. Algoritmo.................................... 13 3.1.3. Sistemas estad´ısticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.4. Condiciones de contorno, convenci´on de im´agenes m´ınimas y radio de corte 16 3.2. Modelomolecular.................................... 19 3.2.1. Campodefuerza................................ 19 3.2.2. Tipode´atomos................................. 20 3.2.3. Detalles de los sistemas y par´ametros de simulaci´on . . . . . . . . . . . . 25 4. Resultados y Discusi´on 29 4.1. Ajustes del campo de fuerza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.2. Estimaci´on de la temperatura de transici´on v´ıtrea . . . . . . . . . . . . . . . . . 30 4.3. Influencia de la velocidad de enfriamiento en la estimaci´on de la Tg........ 34 4.4. Dependencia de la Tgen el peso molecular (Mn) .................. 37 vi ´ Indice general 4.5. Discusi´on sobre los datos te´oricos y los datos presentes en las fichas t´ecnicas del PLAparaimpresi´on3D ................................ 41 5. Conclusiones 43 6. Perspectivas de futuro 45 Bibliograf´ıa 45 Ap´endice Anexo 51 1. Introducci´on La relaci´on entre la estructura, las propiedades, el procesamiento y el rendimiento de los pol´ımeros, es fundamental para su estudio. La organizaci´on molecular en los pol´ımeros determina sus propiedades y comportamiento, lo que implica directamente en su rendimiento en diversas aplicaciones industriales. Mejorar estas propiedades depende en gran medida de como se gestionan estos factores, por lo que comprender su interrelaci´on es crucial para avanzar e innovar en este campo. En este trabajo, con el fin de comprender mejor esta relaci´on, se comienza con una revisi´on detallada de los materiales polim´ericos (Sec. 1.1). A continuaci´on, se tratar´an los pol´ımeros sostenibles (Sec. 1.2), haciendo un especial ´enfasis en el ´acido polil´actico (PLA), en su uso en la impresi´on 3D, las t´ecnicas computacionales aplicadas a la ciencia de los pol´ımeros y las simulaciones atom´ısticas para este material. 1.1. Materiales polim´ericos Los pol´ımeros son materiales macromoleculares. El t´ermino pol´ımeros proviene del griego “poly” y “mers” que significa “muchas partes” y se refiere a mol´eculas que consisten en muchas unidades elementales qu´ımicas, llamadas mon´omeros. Estos ´ultimos son unidades de repetici´on estructurales que est´an conectadas entre s´ı por enlaces covalentes [1]. 1.1.1. Estructura La estructura de un pol´ımero se genera durante la polimerizaci´on, el proceso por el cual las unidades elementales (mon´omeros) se unen covalentemente. El n´umero de mon´omeros en una mol´ecula se denomina grado de polimerizaci´on (N). La masa molar (Mn) de un pol´ımero es igual a su grado de polimerizaci´on por el peso molecular del mon´omero (Mmon) [1]: Mn=NMmon (1.1) Cuando en un pol´ımero determinado, todos los mon´omeros son id´enticos, este recibir´a el 8Cap´ıtulo 1. Introducci´on previos al reproducir correctamente propiedades estructurales y t´ermicas, resultando ideal para estudiar fen´omenos moleculares en materiales basados en PLA [16]. Sin embargo, es necesario destacar su principal desventaja, que es el uso de potenciales externos (tabulados), que carecen de flexibilidad, versatilidad y hacen que el proceso de c´alculo de fuerza sea m´as lento. En el estudio de Guseva, Lazutin y Vaselevskaya, se realizaron simulaciones de din´amica molecular del PLA amorfo y cristalino usando el campo de fuerza GAFF. Simularon sistemas con distintas longitudes de cadena (13, 54 y 150 mon´omeros) y contenido de unidades D(0-50 %). La Tg fue estimada mediante simulaciones directas ´unicamente para el homopol´ımero (PLLA), mientras que para los copol´ımeros con contenido de unidades D se utiliz´o una aproximaci´on emp´ırica basada en la ecuaci´on de Flory-Fox, la cual predice una disminuci´on de T g al aumentar la fracci´on de unidades D. Sin embargo, los valores simulados de Tg resultaron ser mayores que los datos experimentales y los predichos emp´ıricamente, lo que se atribuye a limitaciones del modelo y a diferentes tasas de enfriamiento [17]. El trabajo de Christofi, Baˇcov´a y Harmandaris presenta una metodolog´ıa computacional basada en redes neuronales U-net para reconstruir estructuras at´omicas a partir de modelos CG de pol´ımeros biodegradables. Su objetivo es facilitar y acelerar el proceso de equilibraci´on de los sistemas atom´ısticos. Utilizan el campo de fuerza PLAFF3 simplificado, y estudian PLLA, PDLA y copol´ımeros con longitudes de cadena de hasta 100 mon´omeros. Los sistemas creados por el algoritmo de aprendizaje autom´atico mostraron propiedades estructurales similares a las estructuras objetivas, reduciendo as´ı el tiempo de simulaci´on de forma significativa. Gracias a su versatilidad, esta herramienta abre la puerta a simulaciones de estructuras PLA m´as largas y complejas, lo que hasta ahora no hab´ıa sido posible debido a sus tiempos computacionales inviablemente largos. [18]. En el art´ıculo de Klajmon, Aulich, Lud´ık y ˇ Cervinka, las simulaciones de MD se realizaron utilizando un campo de fuerza PLAFF3. El rango de temperatura estudiadas fue de 500 a 1000 K, con sistemas compuestos por 40 cadenas y una longitud mayor de 5000 g/mol, cercanas a 10,000 g/mol. Para poder medir la T g, se aplicaron tasas de enfriamiento de 40 y 5 K/ns, obteni´endose un valor de T g cercano al experimental, con una desviaci´on de 20 K. La T g se estim´o mediante dos m´etodos: un ajuste hiperb´olico de los datos de equilibrio y el m´etodo de intersecci´on de equilibrio, ambos aplicados a las densidades simuladas en funci´on de la temperatura [15]. 2. Objetivos En el mundo de los pol´ımeros biodegradables, el ´acido polil´actico (PLA) brilla por su origen renovable y su creciente uso en campos como la impresi´on 3D, el envasado y la biomedicina. Sin embargo, su adopci´on a gran escala todav´ıa enfrenta algunos obst´aculos, ya que sus propiedades t´ermicas y mec´anicas son bastante sensibles y dependen de factores estructurales como el peso molecular y la estereoqu´ımica. Ante estas limitaciones, los m´etodos computacionales se presentan como una alternativa sostenible para explorar de manera controlada la relaci´on entre la estructura y las propiedades. El objetivo principal de este trabajo es desarrollar un enfoque computacional sistem´atico que permita analizar el comportamiento t´ermico del PLA, tanto en su forma de homopol´ımero como en diversas configuraciones de copol´ımeros. Para lograrlo, se propone recopilar y analizar la literatura cient´ıfica existente sobre la temperatura de transici´on v´ıtrea (Tg) del PLA, una propiedad clave desde perspectivas tanto cient´ıficas como industriales. Adem´as, se busca seleccionar y adaptar un campo de fuerza vers´atil que represente adecuadamente las caracter´ısticas estructurales del PLA, incluyendo su estereoqu´ımica, para simular sistemas con diferentes pesos moleculares, proporciones de mon´omeros L/D y velocidades de enfriamiento. Uno de los objetivos principales de este trabajo es realizar un an´alisis cuanlitativo sobre c´omo var´ıa la Tgen funci´on del peso molecular y la composici´on estereoqu´ımica de las cadenas polim´ericas. Esto nos permitir´a establecer una relaci´on clara entre la estructura y las propiedades t´ermicas. Adem´as, de validar los resultados obtenidos a trav´es de simulaciones de din´amica molecular, compar´andolos con datos en la literatura. Esto nos ayudar´a a evaluar la precisi´on del modelo utilizado. As´ı, el estudio no solo profundizar´a en nuestra comprensi´on del PLA, sino que tambi´en contribuir´a al desarrollo de herramientas computacionales que sean ´utiles para el dise˜no racional de materiales biodegradables con aplicaciones pr´acticas en la industria. Asimismo, este enfoque computacional funciona como un puente entre la teor´ıa y la pr´actica, ya que permite validar modelos utilizando datos simulados y facilita la interpretaci´on de los resultados experimentales. As´ı, las simulaciones ayudan a disminuir tanto el costo como la 10 Cap´ıtulo 2. Objetivos complejidad de los ensayos f´ısicos, ofreciendo una herramienta predictiva valiosa para el dise˜no y la optimizaci´on de materiales. 3. Materiales y M´etodos 3.1. Din´amica Molecular La Din´amica Molecular (MD) es un m´etodo computacional que hace posible simular la evoluci´on de los sistemas formados por m´ultiples ´atomos o mol´eculas, constituyendo, junto al desarrollo de la potencia de c´alculo en las ´ultimas d´ecadas, una t´ecnica computacional muy utilizada para estudiar las propiedades de equilibrio y las din´amicas de los sistemas de muchos cuerpos, hasta el punto de convertirse en la t´ecnica computacional m´as utilizada en el estudio de los materiales polim´ericos [19]. En t´erminos sencillos, es una forma de reproducir un sistema de la vida real, donde se puede seguir cada movimiento de cada parte del sistema, conociendo d´onde est´an, con qu´e velocidad se mueven los elementos, y otros factores importantes. Emerge como una herramienta esencial que todo investigador debe conocer y aplicar. Este m´etodo sigue el mismo enfoque de un experimento real, permitiendo realizar el proceso repetidamente, ajustando o conservando las condiciones iniciales seg´un se requiera. Adem´as, las simulaciones pueden interrumpirse y retomarse con facilidad, ofreciendo un control m´as preciso sobre las variables del sistema que en un laboratorio convencional [20]. La din´amica molecular investiga un amplio rango de sistemas fisicoqu´ımicos con una debida condici´on inicial, a partir de la interacci´on de un n´umero de part´ıculas, ´atomos o mol´eculas a trav´es de un campo de fuerzas o de potencial interat´omico. Con ello, se trata de determinar posiciones y velocidades, mediante la integraci´on num´erica de la ecuaci´on de Newton [21]. El proceso del c´alculo de las magnitudes necesarias en una simulaci´on se puede realizar conociendo la trayectoria de los grupos de ´atomos en el sistema. Dicha trayectoria corresponde a la salida principal de la simulaci´on y est´a dada por las posiciones y las velocidades, en cada uno de los pasos de la misma. El c´alculo de la trayectoria se realiza resolviendo las ecuaciones del movimiento cl´asico (las ecuaciones de Newton) que existen para cada mol´ecula o ´atomo del sistema, teniendo en cuenta las interacciones entre ellas. 12 Cap´ıtulo 3. Materiales y M´etodos 3.1.1. Ecuaciones de movimiento La determinaci´on actualizada de las posiciones y velocidades a medida que avanza el tiempo en un sistema conformado por distintos ´atomos se encuentra dictada por las ecuaciones de movimiento, que en el formalismo Hamiltoniano tienen la forma siguiente [19]: ˙qi=∂H ∂pi ,˙pi=−∂H ∂qi (3.1) donde el hamiltoniano est´a definido como: H(p,q) = X i ˙qipi−L. (3.2) Denotando qicomo las coordenadas generalizadas que describen la configuraci´on molecular y ˙qi sus derivadas temporales, piser´ıa el momento generalizado y ˙pisus derivadas. En el contexto de formalismo lagrangiano, el Lagrangiano Lse define como L(q,˙ q, t) = K(˙ q(t)) −U(q(t)), donde Krepresenta la energ´ıa cin´etica y Ula energ´ıa potencial del sistema. En situaciones donde la energ´ıa de potencial Ues independiente de tiempo [19], el Hamiltoniano Hrepresentar´ıa la energ´ıa total del sistema: H(p,q) = K(p) + U(q).(3.3) En el caso de coordenadas cartesianas, donde la energ´ıa cin´etica se define como K(˙ r) = 1 2m˙ r2 (como ˙rrepresentando las velocidades en las direcciones cartesianas) y la energ´ıa potencial se expresa como U(r), las ecuaciones de movimiento ( 3.1) para npart´ıculas denominadas como i= 1,2, ..., n, a trav´es de la ecuaci´on 3.3 se transforman de la siguiente manera: ˙ ri=pi mi =vi˙ pi=−∇riU(r) (3.4) donde rirepresenta la posici´on cartesiana de cada part´ıcula iy˙ ries la derivada temporal de ri, es decir vila velocidad del ´atomo iymisu masa. La dificultad de obtener la resoluci´on de las ecuaciones de forma anal´ıtica hace necesario optar por m´etodos num´ericos como el de las diferencias finitas. Estos m´etodos demandan un paso de tiempo discreto (∆t > 0) para calcular la posici´on y la velocidad en intervalos definidos. Por lo que se hacen necesarios integradores como el Velocity Verlet o el Leapfrog que son entre los m´as utilizados. Estos integradores son capaces de ofrecer valores que son coherentes y cercanos al resultado exacto, son sencillos de implementar y son perfectamente adaptables en situaciones y tipos de sistemas diferentes [22]. Velocity Verlet: destaca en simulaciones num´ericas por su capacidad de gestionar tanto la posici´on como la velocidad en un mismo intervalo temporal e integrar de forma precisa, 3.1. Din´amica Molecular 13 especialmente cuando se acoplan la temperatura y/o presi´on. En comparaci´on con integradores como el Leapfrog, resuelve la imprecisi´on en el tratamiento de la velocidad satisfactoriamente [22, 20]. Leapfrog: se basa en la utilizaci´on de coordenadas y velocidades en pasos o tiempos alternos y determina las velocidades en medio paso, en comparaci´on con el Velocity Verlet. Es especialmente ´util en casos donde se requiere conservar la energ´ıa del sistema a lo largo de la simulaci´on. Es un m´etodo de integraci´on muy preciso [23]. A fin de analizar la trayectoria obtenida y hallar las propiedades de la simulaci´on, se asume que el sistema es erg´odico. Esta es una propiedad que considera que los promedios que se obtienen a lo largo del tiempo en la evoluci´on de un sistema son equivalentes a los promedios que se obtienen sobre un conjunto estad´ıstico de estados iniciales, conocido como ensemble. Es un requerimiento fundamental en la mec´anica estad´ıstica del equilibrio, por ser uno de los propios fundamentos del modelo de equilibrio termodin´amico [20]. En la din´amica molecular su base se fundamenta en que, en lugar de realizar un promedio ponderado de diferentes muestras, se realiza un promedio de tiempo de una simulaci´on larga, es decir, la propiedad se mide por un cierto tiempo considerando que lo sea estad´ısticamente, como si se tuvieran varias muestras independientes. 3.1.2. Algoritmo El proceso de exportar la trayectoria y las variables relevantes en una simulaci´on se asemeja al procedimiento de un experimento real. Primero, se prepara el sistema bajo estudio, estableciendo las condiciones iniciales adecuadas. Luego, se conecta el sistema simulado con los algoritmos que permiten registrar las magnitudes de inter´es. A lo largo de la simulaci´on, estas cantidades se registran en intervalos definidos, con el prop´osito de calcular promedios estad´ısticos que caractericen el comportamiento del sistema. El procedimiento pertinente en simulaciones por ordenador es amplio y pr´acticamente uniforme para diversos tipos de simulaciones de MD. Se conoce como algoritmo MD (figura 3.1), y su estructura est´a basada en los siguientes pasos [20, 19]: 1. Inicialmente, es necesario conocer cierta configuraci´on, posiciones y velocidades de cada part´ıcula, y el potencial entre ellas, es decir, su interacci´on. La configuraci´on inicial puede partir de varios sitios, como una estructura cristalina y la velocidad se puede ajustar por ejemplo seg´un la distribuci´on Maxwell-Boltzmann (ecuaci´on 3.5). 2. Se calculan las fuerzas que ejercen su acci´on sobre cada part´ıcula en el sistema, teniendo en cuenta las condiciones externas (temperatura, presi´on...). 14 Cap´ıtulo 3. Materiales y M´etodos Inicialización del sistema Calcular y almacenar las cantidades importantes Presentar el resultado final Iteraciones K Figura 3.1: Representaci´on gr´afica del algoritmo MD. 3. Se selecciona un m´etodo num´erico para integrar las ecuaciones cl´asicas de movimiento, con el prop´osito de obtener las nuevas posiciones y velocidades. Luego se almacenan estos valores y se utilizan para reemplazar los iniciales. 4. Se realiza el c´alculo y almacenamiento de cantidades significativas, como la energ´ıa y la densidad, utilizando la posiciones y velocidades actuales. 5. Se exponen los valores finales de salida. En el primer paso, es necesario inicializar las velocidades, ya que para el tiempo cero, estas deben tener un valor que se inserta en el integrador. La elecci´on m´as com´un, sobre todo en los sistemas con la temperatura constante, es ajustar las velocidades seg´un una distribuci´on de Maxwell-Boltzmann dada por: P(v) = m 2πkBT3 2 e −mv2 2kBT(3.5) 3.1. Din´amica Molecular 15 donde vrepresenta la velocidad de la part´ıcula, msu masa, kBla constante de Boltzmann y Tla temperatura del sistema. Durante la simulaci´on, los pasos que corresponden desde el segundo hasta el cuarto, se repiten, tantas veces como sea necesario (digamos x). El valor entero xest´a condicionado por la naturaleza del sistema y las propiedades que se deseen explorar. La configuraci´on inicial y las condiciones correspondientes que definen el tipo de procedimiento se ajustan conforme al problema espec´ıfico bajo investigaci´on 3.1. 3.1.3. Sistemas estad´ısticos En las simulaciones de MD, se intenta representar un sistema grande usando un n´umero reducido de mol´eculas, por lo general menos de un mill´on. Aunque se espera que los resultados de la simulaci´on sean similares a los obtenidos en experimentos reales a gran escala, es importante asegurarse de que las simulaciones sean correctas a nivel molecular. Una forma de hacerlo es verificando que las trayectorias de las part´ıculas, sus posiciones y movimientos en el tiempo, se ajusten con precisi´on a lo que predice la teor´ıa de la mec´anica estad´ıstica. El sistema estad´ıstico “natural”, para la MD es el sistema microcan´onico NVE (n´umero de part´ıculas N, energ´ıa Ey volumen Vconstantes) que emplea un algoritmo bien aceptado y simple, aunque en la pr´actica resulta dif´ıcil su empleo, ya que requiere mantener constante la energ´ıa y en un sistema real no es constante [20]. Aunque el uso directo de la MD da lugar al NVE, a veces existe la necesidad de evolucionar un sistema molecular bajo condiciones espec´ıficas. La mayor´ıa de las cantidades que se desean calcular son en realidad sistemas de temperatura constante NVT (n´umero constante de part´ıculas, volumen y temperatura T), tambi´en llamado sistema can´onico donde el algoritmo que se utiliza es muy riguroso. El m´as parecido a un experimento com´un ser´ıa un sistema a presi´on y temperatura constantes como es el caso del sistema estad´ıstico NPT (n´umero constante de part´ıculas, presi´on Py temperatura), denominado tambi´en isob´arico-isot´ermico (Gibbs) [20]. Para mantener constantes la temperatura y la presi´on se utilizan los termostatos y bar´ostatos adecuados, respectivamente. En la simulaci´on es una forma de mantener ciertas condiciones externas constantes como en experimentos. Termostato Nos´e-Hoover Este procedimiento fue propuesto inicialmente por Nos´e y posteriormente adoptado por Hoover. Se trata del m´etodo que permite mantener la temperatura interna del sistema, el cual simula condiciones requeridas de equilibrio t´ermico (NVT, sistema can´onico). Ampl´ıa el sistema hamiltoniano introduciendo un nuevo grado de libertad con la finalidad de llevar el control y 16 Cap´ıtulo 3. Materiales y M´etodos la regulaci´on de la temperatura. Se puede explicar por medio de un t´ermino adicional que representa un ba˜no t´ermico que interact´ua con dicho sistema para regular las fluctuaciones de temperatura a trav´es de la transferencia de calor, asegurando que las fluctuaciones de temperatura sean controladas, por consiguiente, la temperatura se mantiene constante en el sistema durante la simulaci´on [24]. Termostato y bar´ostato Berendsen Dicho algoritmo se presentar´ıa como otra alternativa de llevar a cabo las simulaciones isot´ermicas y/o isob´arica de la MD, para lo cual se recurre a un lagrangiano extendido que acopla al sistema una temperatura y/o al ba˜no una presi´on. El efecto del algoritmo radica en que se va corrigiendo lentamente la desviaci´on de la temperatura suprimiendo las fluctuaciones de la energ´ıa cin´etica (relajar un sistema a la temperatura que queremos). En relaci´on con la presi´on, el algoritmo realiza una relajaci´on muy r´apida, sin oscilaciones, hacia una presi´on de referencia [25]. 3.1.4. Condiciones de contorno, convenci´on de im´agenes m´ınimas y radio de corte Mediante la simulaci´on por ordenador utilizando el m´etodo de din´amica molecular, las part´ıculas se introducen en una caja de simulaci´on, generalmente cuadrada y preferiblemente de dimensiones grandes para que haya suficiente espacio libre para el movimiento de las part´ıculas y efectos de tama˜no finito no puedan influir. Ya que el tama˜no de la caja est´a ligado de manera ´ıntima con el tiempo computacional no ser´a posible simular sistemas a macroescala. Se puede usar una caja de simulaci´on mucho m´as peque˜na, pero las mol´eculas estar´ıan en su inmensa mayor´ıa a un borde de la caja. La forma cl´asica de afrontar este problema o minimizar los efectos de borde en un sistema finito consiste en implementar condiciones de contorno peri´odicas [20] (PBC por sus siglas en ingl´es, Periodic Boundary Conditions). Bajo estas condiciones, se establece un sistema de red infinito al replicar la caja de simulaci´on en todo el espacio, estas r´eplicas se denominan im´agenes de caja. Cuando una part´ıcula abandona la caja, una de sus im´agenes ingresar´a por el lado opuesto en la misma direcci´on, tal y como se muestra en la figura 3.2. Las part´ıculas contenidas en la caja de simulaci´on se conservar´an y se puede pensar que el sistema carece de una superficie definida [20]. Seg´un el criterio de convenci´on de im´agenes m´ınimas (MIC por sus siglas en ingl´es, Minimum Image Convention) [12], cada part´ıcula interact´ua con la imagen de la part´ıcula que se encuentra m´as cercana a ella y no a trav´es de la caja. La distancia entre las part´ıculas en una direcci´on 3.1. Din´amica Molecular 17 Figura 3.2: Condiciones peri´odicas de contorno. La caja de simulaci´on (la celda principal) est´a pintada de verde, las l´ıneas discontinuas que se encuentran a su alrededor son sus copias. La flecha roja muestra como se realiza el movimiento de la part´ıcula verde hacia el circulo en l´ıneas discontinuas. Figura 3.3: Convenci´on m´ınima de imagen. La caja de color azul es la caja principal y la de color verde muestra la part´ıcula en color azul que interact´uan solo con la imagen peri´odica m´as cercana de las otras part´ıculas mostradas en verde. determinada, no puede superar la mitad de las longitudes de arista correspondientes de dicha caja. Es importante evitar doble conteo durante la evaluaci´on de la fuerza. Esta convenci´on se muestra en la figura 3.3. La distancia de corte se emplea para limitar el c´alculo de las interacciones a largas distancias 24 Cap´ıtulo 3. Materiales y M´etodos k j i l θijkl Figura 3.9: ´ Angulo diedro impropio entre i, j, k, l, donde el plano verde con transparencia pasa por j, k, l y el plano azul claro pasa a trav´es de i, j, k. El potencial de Lennard-Jones (ULJ ) modela la interacci´on Van der Waals, siendo atractivo a larga distancia y repulsivo a corta distancia entre pares de ´atomos. Por otro lado, las interacciones electrost´aticas entre ´atomos cargados se describen mediante el potencial de Coulomb (UC). En este estudio, las interacciones de no enlace entre ´atomos directamente conectados no se incluyen en el c´alculo si las part´ıculas est´an separadas por menos de 2 ´atomos. En otras palabras, s´olo ´atomos que est´an a distancias m´as grandes a lo largo de la cadena, que dos ´atomos, est´an incluidos en el c´alculo de los potenciales no enlazados. Las interacciones 1-4, que se refieren a las fuerzas entre ´atomos separados por tres enlaces consecutivos, fueron reescaladas seg´un McAliley [16] con un factor de escala de 0.5, reduciendo as´ı la intensidad de las interacciones. Interacci´on de Lennard-Jones El potencial de Lennard-Jones describe las interacciones entre pares de ´atomos en funci´on de la distancia rij entre ellos, utilizando par´ametros espec´ıficos que dependen de los tipos de ´atomos involucrados. Esta interacci´on est´a definida por: ULJ (rij)=4ϵij σij rij 12 −σij rij 6!(3.13) donde ϵij es la profundidad del pozo de energ´ıa potencial o energ´ıa de dispersi´on y σij la distancia a la que el potencial es cero. En nuestro caso, al parametrizar los potenciales de Lennard-Jones no enlazados, se determinan los valores espec´ıficos de σij yϵij para la interacci´on entre dos ´atomos iyj. Esto se logra combinando los par´ametros individuales σiyϵide cada ´atomo mediante promedios geom´etricos, 3.2. Modelo molecular 25 como se muestra en la ecuaci´on 3.14. Esta t´ecnica de combinaci´on de par´ametros es utilizada por el campo de fuerza OPLS [30]: σij = (σii σjj)1/2;ϵij = (ϵii ϵjj)1/2(3.14) La ecuaci´on 3.13 describe las interacciones atractivas de Van der Waals mediante el t´ermino negativo. Se adiciona el t´ermino repulsivo al potencial para evitar que sus dos ´atomos i, j se penetren mutuamente cuando su distancia es menor que la suma de sus radios at´omicos (σi 2+σj 2= σl). Para distancias rij mayores que σl, la fuerza entre los ´atomos iyjes atractiva, pero cuando rij es mucho mayor que σl, los ´atomos est´an demasiado separados para interactuar y la energ´ıa potencial de LJ se acerca a cero. Por otro lado, para distancias rij menores que σl, la fuerza se vuelve repulsiva, y cuando rij se acerca a cero, la energ´ıa potencial de LJ tiende a infinito debido a la superposici´on de los ´atomos. Interacci´on de Coulomb La interacci´on electrost´atica de Coulomb entre dos part´ıculas (´atomos) iyj, separadas por una distancia rij y con cargas QiyQjrespectivamente, se describe mediante la f´ormula dada por la ecuaci´on 3.15. En esta ecuaci´on, UC(rij) representa el potencial electrost´atico, ϵ0es la permitividad el´ectrica del vac´ıo, y ϵrson los relativos adimensionales de la permitividad del medio. UC(rij) = QiQj 4πϵ0ϵrrij (3.15) 3.2.3. Detalles de los sistemas y par´ametros de simulaci´on Se realizaron simulaciones de din´amica molecular en este trabajo, para determinar la temperatura de transici´on v´ıtrea (Tg), de diferentes sistemas de ´acido polil´actico (PLA) incluyendo homopol´ımeros y copol´ımeros de distinto contenido de cada tipo de mon´omero (L/D). Los sistemas simulados incluyeron homopol´ımeros de L-lactida (PLLA) con cadenas de 10, 30 y 100 mon´omeros (denominados PLLA10, PLLA30 y PLLA100, respectivamente), as´ı como un homopol´ımero de D-lactida (PDLA) con 100 mon´omeros. Tambi´en se estudiaron dos copol´ımeros con composiciones estereoqu´ımicas de 84/16 % y 45/55 % L/D, denominados copo16D y copo55D, respectivamente. Se simul´o el enfriamiento lineal controlado de sistemas fundidos desde 500 K hasta 200 K, aplicando diferentes velocidades de enfriamiento (desde 600 K/ns hasta 2.5 K/ns). Las velocidades espec´ıficas est´an asociadas a distintos tiempos de simulaci´on, desde 0.5 ns en los casos m´as r´apidos hasta 120 ns en los m´as lentos. El listado de los sistemas junto con las 26 Cap´ıtulo 3. Materiales y M´etodos Sistema NL/D Mn Velocidad de enfriamiento Tiempos de simulaci´on [ %] [g/mol] [K/ns] [ns] PLLA10 10 100/0 738.7 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 PLLA30 30 100/0 2179.9 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 PLLA100 100 100/0 7224.4 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 PDLA100 100 0/100 7224.4 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 copo16D 100 84/16 7224.4 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 copo55D 100 45/55 7224.4 600/300/100/40/25/10/5/2.5 0.5/1/3/7.5/12/30/60/120 Tabla 3.1: Caracterizaci´on de los sistemas de PLA simulados: n´umero de mon´omeros (N), composici´on estereoqu´ımica (L/D), masa molar (Mn), perfiles de enfriamiento (velocidades aplicadas en el rango de 200-500 K) y duraci´on de las simulaciones de din´amica molecular (MD). simulaciones realizadas se encuentra en la Tabla 3.1. Para modelar las interacciones a nivel at´omico, se utiliz´o un enfoque all-atom con el campo de fuerza PLAFF3, basado en una versi´on modificada del OPLS-AA adaptado espec´ıficamente para PLA [16]. Entre las modificaciones de este campo de fuerza en nuestro trabajo est´an: los par´ametros de enlaces, ´angulos y dihedros, con especial atenci´on a las interacciones estereoqu´ımicas entre unidades L y D. Adem´as, se incorpor´o un potencial tipo CMAP con el objetivo de mejorar la descripci´on de las conformaciones locales de la cadena polim´erica. Antes de someter los materiales a un enfriamiento lineal, los sistemas hab´ıan sido equilibrados previamente seg´un la metodolog´ıa publicada en la ref. [18]. Como el procedimiento de equilibraci´on no es parte de este trabajo, s´olo se resume brevemente: cada sistema consisti´o en 70 cadenas de PLA colocadas en una caja tridimensional con condiciones de contorno peri´odicas (PBC) en las tres direcciones espaciales. Para la generaci´on de cadenas de longitud deseada (30 o 100 mon´omeros), se parti´o de cadenas largas de 500 unidades que fueron recortadas y reequilibradas. En el caso de los estereois´omeros D o copol´ımeros, se modific´o la estereoqu´ımica mediante el intercambio del grupo metilo y el hidr´ogeno en el carbono quiral de algunas cadenas seleccionadas aleatoriamente. El protocolo de equilibraci´on se desarroll´o en varias etapas: eliminaci´on de heterogeneidades de densidad, mezcla a alta temperatura, enfriamiento progresivo y finalmente estabilizaci´on en condiciones objetivo. Este proceso es cr´ıtico debido a que la relajaci´on estructural en sistemas macromoleculares es lenta y depende fuertemente del peso molecular. Las simulaciones de enfriamiento realizadas en este trabajo se iniciaron partiendo desde con- 3.2. Modelo molecular 27 figuraciones fundidas pre-equilibradas a 500 K. Se realizaron en el ensamble NPT, manteniendo una presi´on y temperatura constante de 1 bar con el barostato Parrinello-Rahman, y controlando la temperatura entre 500 y 200 K con el termostato Nose-Hoover. Las velocidades at´omicas iniciales se generaron conforme a una distribuci´on de Maxwell-Boltzmann a 500 K. Se utiliz´o el integrador Leapfrog por su estabilidad en simulaciones de larga duraci´on, y los enlaces con ´atomos de hidr´ogeno restringidos mediante el algoritmo LINCS (H-bonds), permitiendo un paso de integraci´on de 1 fs. Las simulaciones moleculares se han realizado en GROMACS [26]. Se llevaron a cabo simulaciones de 120 y 60 ns en un cl´uster en Grecia, mientras que el resto se ejecut´o en el cl´uster de la Universidad de C´adiz (UCA). Las simulaciones m´as largas, que involucraron cadenas largas (100 mon´omeros) durante 120 ns, necesitaron alrededor de 10 d´ıas de c´alculo con 24 procesadores. Las simulaciones en el cl´uster local se realizaron con una GPU, cuyo rendimiento vari´o significativamente. El rendimiento m´aximo alcanzado fue de 17 ns/d´ıa. 28 Cap´ıtulo 3. Materiales y M´etodos 4. Resultados y Discusi´on 4.1. Ajustes del campo de fuerza Como se mencion´o en la secci´on 3.2.2, en t´erminos del campo de fuerza, los potenciales publicados en ref. [16] representan el estado del arte. Estos potenciales indican los par´ametros necesarios para simular los homopol´ımeros, que consisten en el mismo tipo de mon´omero. En concreto, en el caso de la rotaci´on alrededor de los enlaces de la cadena principal, implementan un potencial tabulado mostrado en la figura 4.1. Para representar adecuadamente la estructura y comportamiento del copol´ımero que incluye en su estructura ambos tipos de mon´omero, L y D, en este trabajo se ha implementado un ajuste espec´ıfico del potencial diedro de la cadena principal del pol´ımero. Para poder comparar de forma justa los resultados de los copol´ımeros y de los homopol´ımeros, los potenciales tabulados en los homopol´ımeros se reemplazaron tambi´en por la forma funcional del potencial diedro (l´ınea continua en la figura 4.1). Para describir de manera precisa la rotaci´on alrededor del enlace principal de los copol´ımeros, se ha utilizado un potencial torsional del tipo Ryckaert-Bellemans (v´ease la ecuaci´on 3.11). Este enfoque permite representar adecuadamente los perfiles de energ´ıa relacionados con la rotaci´on, especialmente en casos donde hay m´ultiples m´ınimos y barreras suaves. Este potencial se ha aplicado al diedro central del mon´omero, espec´ıficamente al que involucra los ´atomos Cα–OS– C–Cα(ve´ase figura 3.5), incluyendo sus combinaciones con el centro quiral Cαtipo D, para abarcar tanto homopol´ımeros como copol´ımeros. Los par´ametros ajustados a partir de perfiles de energ´ıa representativos son los siguientes: C0= 65.93390 kJ/mol, C1= 8.48321 kJ/mol, C2 = –75.20430 kJ/mol, y C3, C4yC5= 0.00000 kJ/mol. Estos valores se han a˜nadido de manera expl´ıcita en la secci´on de tipos de diedro en el archivo con los potenciales de enlace. Aparte de implementar la forma funcional del potencial diedro, para representar correctamente el acoplamiento estereoqu´ımico en los copol´ımeros, ha sido necesario diferenciar entre los centros quirales de los mon´omeros de tipo L y D. Para ello, se ha creado un nuevo tipo de ´atomo (opls 491d) que permite distinguir los entornos locales seg´un la quiralidad, adem´as de 30 Cap´ıtulo 4. Resultados y Discusi´on los tipos ya existentes que se mostraron en la secci´on anterior. Estos ´atomos se han incorporado manualmente en la secci´on de la estructura del archivo de topolog´ıa. Es importante destacar que para lograr el objetivo de este trabajo fue esencial reemplazar el potencial tabulado, debido a que el paquete de simulaci´on empleado, Gromacs, no permite el uso de dos potenciales tabulados distintos, lo cual ser´ıa necesario en el caso de copol´ımeros. Tambi´en, el uso de los potenciales tabulados resultar´ıa en una simulaci´on m´as lenta y menos precisa debido a la interpolaci´on de las tablas. Por lo tanto, se ha elegido una soluci´on m´as eficiente y robusta que se basa exclusivamente en par´ametros anal´ıticos, evitando el uso de tablas y asegurando la precisi´on y el rendimiento de la simulaci´on. −20 −10 0 10 20 30 40 50 60 70 80 −150 −100 −50 0 50 100 150 Distribución Ángulo [Grados] fit PDLA PLLA Figura 4.1: Distribuci´on de ´angulos dihedrales en sistemas de PLA. Potenciales tabulados para los dos homopol´ımeros est´an representados con los puntos, el ajuste usado en este trabajo con la l´ınea continua. 4.2. Estimaci´on de la temperatura de transici´on v´ıtrea Para estimar la Tgse usaron dos m´etodos: intersecci´on de equilibrio o l´ıneas secantes y tangente hiperb´olica. M´etodo de intersecci´on de equilibrio o l´ıneas secantes La temperatura de transici´on v´ıtrea (Tg) fue estimada a partir de la evoluci´on del volumen espec´ıfico (v) en funci´on de la temperatura (T) (v´ease ejemplos en la figura S1), mediante el m´etodo de intersecci´on de l´ıneas secantes. Esta t´ecnica consiste en identificar dos tramos lineales en la curva v(T): uno corresponde al estado m´ovil en las altas temperaturas y otro al estado v´ıtreo a bajas temperaturas. A cada tramo se le ajusta una funci´on lineal del tipo: f(T) = aT +b(estado v´ıtreo), g(T) = cT +d(estado m´ovil) 4.2. Estimaci´on de la temperatura de transici´on v´ıtrea 31 La Tgse define como el punto de intersecci´on entre ambas rectas, calculado mediante la igualdad: f(Tg) = g(Tg)⇒aTg+b=cTg+d(4.1) de donde se despeja: Tg=d−b a−c.(4.2) La figura 4.2 ilustra de manera gr´afica el proceso utilizado para estimar la Tg. En ella, se muestran los ajustes lineales realizados en las ´areas que corresponden al estado v´ıtreo y al estado m´ovil, as´ı como el punto donde ambas l´ıneas se cruzan, el cual se se˜nala con l´ıneas discontinuas. Se presentan dos casos representativos: PDLA 100 % D y el copol´ımero 55 % D. Un ejemplo adicional del procedimiento se puede ver en la figura S2 del anexo. Se utiliz´o una herramienta de visualizaci´on y ajuste de datos, concretamente el programa gr´afico de c´odigo abierto gnuplot, empleando el algoritmo de m´ınimos cuadrados no lineal de Marquardt-Levenberg. Este procedimiento permiti´o obtener los coeficientes (a,b,cyd) a partir de subconjuntos de datos seleccionados visualmente en las regiones lineales correspondientes. Es importante se˜nalar que este procedimiento no es exacto, ya que la elecci´on de los intervalos de ajuste puede introducir variabilidad en la estimaci´on de la Tgdebido a su dependencia del criterio del usuario. Para cuantificar la incertidumbre asociada al m´etodo, se relacionaron cinco repeticiones del proceso para cada sistema (homopol´ımeros y copol´ımeros), utilizando diferentes subconjuntos de datos dentro de las regiones lineales. Con las cinco estimaciones de Tg,i obtenidas, se calcul´o un promedio y su desviaci´on est´andar mediante las siguientes expresiones: ¯ Tg=1 n n X i=1 Tg,i (4.3) σ=v u u t 1 n−1 n X i=1 (Tg,i −¯ Tg)2(4.4) donde n= 5 es el n´umero de repeticiones. Teniendo en cuenta los seis sistemas de estudio, que incluyen cuatro homopol´ımeros y dos copol´ımeros, y considerando ocho velocidades de enfriamiento diferentes para cada uno, se realizaron cinco repeticiones independientes para cada combinaci´on de sistema y velocidad. En total, esto suma 6 ×8×5 = 240 estimaciones individuales de Tg, lo que permite llevar a cabo un an´alisis estad´ıstico s´olido y una evaluaci´on confiable de la incertidumbre asociada al m´etodo utilizado. Este enfoque es ampliamente utilizado debido a su simplicidad y a su correspondencia con m´etodos experimentales, como la calorimetr´ıa diferencial de barrido (DSC). No obstante, dadas 32 Cap´ıtulo 4. Resultados y Discusi´on sus limitaciones, se recurri´o adem´as a un m´etodo alternativo, basado en el ajuste de una funci´on tangente hiperb´olica, con el objetivo de obtener valores de Tgm´as precisos y con menor dispersi´on, como se detallar´a a continuaci´on. 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 200 250 300 350 400 450 500 (a) Volumen específico [cm3/g] Temperatura [K] 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 200 250 300 350 400 450 500 (b) Volumen específico [cm3/g] Temperatura [K] Figura 4.2: La estimaci´on de la Tgrealizada mediante el ajuste de l´ıneas secantes en dos de los casos analizados: (a) PDLA a 2.5 K/ns y (b) Copo55 a 5 K/ns. Los puntos en el gr´afico representan los datos obtenidos de la simulaci´on. Se trazan las l´ıneas que corresponden al estado v´ıtreo (l´ınea lila) y al estado m´ovil (l´ınea verde), y se determina Tgcomo el punto donde ambas l´ıneas se cruzan. Para el caso (a), se obtiene Tg≈394.3 K, mientras que para el caso (b), Tg≈ 441.2 K. M´etodo de tangente hiperb´olica Este m´etodo permite determinar la Tgcomo el punto de inflexi´on de la curva de la derivada, es decir, el centro entre los estados m´ovil y v´ıtreo del pol´ımero. La funci´on empleada para el ajuste tiene la forma: f(T) = b+Atanh T−Tg 2w(4.5) donde Arepresenta la amplitud del cambio en la derivada, wes un par´ametro relacionado con el ancho en la transici´on y bes el valor de desplazamiento vertical de la curva. Para aplicar este procedimiento, se parti´o de los datos de volumen espec´ıfico en funci´on de la temperatura obtenidos a partir de las simulaciones. Con el objetivo de reducir el ruido num´erico, se aplic´o un proceso de suavizado inicial utilizando un promedio m´ovil. Este suavizado se llev´o a cabo mediante la funci´on convolve de Python, empleando un filtro de 50 puntos y el modo valid, que evita los efectos de borde en los extremos de la se˜nal. A continuaci´on, se calcul´o la derivada num´erica del volumen espec´ıfico con respecto a la temperatura, y se volvi´o a suavizar esta nueva curva, generando as´ı una funci´on continua y m´as estable para el ajuste. 4.2. Estimaci´on de la temperatura de transici´on v´ıtrea 33 El ajuste de la funci´on tangente hiperb´olica a la curva derivada suavizada se realiz´o mediante m´ınimos cuadrados no lineales implementados en Python con la funci´on least squares. Esta funci´on tiene la ventaja de poder configurar los l´ımites de cada variable ajustada, lo que permite descartar los resultados no f´ısicos. Este enfoque permiti´o modelar con precisi´on la transici´on t´ermica, identificando directamente la Tgcomo el par´ametro correspondiente al punto medio de la funci´on ajustada, v´ease figura 4.3. Un ejemplo adicional del procedimiento se puede ver en la figura S3 del anexo. Para cuantificar la incertidumbre asociada a esta estimaci´on, se us´o el par´ametro wcomo una estimaci´on del ancho efectivo de la transici´on. Este valor se interpret´o como una medida del error en la estimaci´on de Tg. Este enfoque resulta especialmente ´util porque no depende tanto de la elecci´on manual de intervalos, lo que mejora la reproducibilidad de los resultados. Adem´as, ofrece una mayor sensibilidad para detectar peque˜nos cambios en las propiedades t´ermicas del sistema. En materiales simples, la transici´on v´ıtrea se manifiesta como un cambio bien definido de pendiente (escal´on) en la curva de densidad o volumen espec´ıfico, como en nuestro caso. Este comportamiento permite aplicar este m´etodo de estimaci´on, ya que asume una transici´on homog´enea, con un punto de inflexi´on claro que representa la Tg. No obstante, en sistemas m´as complejos como los copol´ımeros, este enfoque ha fallado debido a la heterogeneidad estructural del material. Al combinar regiones con diferentes movilidades: algunas que se relajan m´as r´apidamente y otras m´as lentamente, esta transici´on t´ermica suele ser m´as difusa. En experimentos, esto se refleja como un pico m´as ancho en la dependencia del flujo de calor frente a temperatura, mientras que en simulaciones puede dificultar la identificaci´on de un cambio claro de pendiente. Al analizar la curva de la derivada de los copol´ımeros, no se observa el comportamiento t´ıpico de dos regiones planas con una inflexi´on bien definida. En su lugar, se observa una curva mon´otonamente creciente, sin un punto de cambio claro, lo que complica la localizaci´on precisa de la Tg, como se puede ver en la figura 4.3b. Esto explica la discrepancia entre los valores estimados por distintos m´etodos y justifica la necesidad de adoptar un enfoque combinado. Teniendo en cuenta los fen´omenos f´ısicos descritos, la transici´on difusa observada en los copol´ımeros podr´ıa indicar que los segmentos correspondientes a los mon´omeros L y D probablemente tengan diferente movilidad. Dado que la movilidad segmentaria est´a estrechamente relacionada con la Tg, esto dar´ıa lugar a la diferencia de Tgentre los homopol´ımeros L y D. De hecho, este es el caso (v´ease la discusi´on a continuaci´on). 40 Cap´ıtulo 4. Resultados y Discusi´on Aunque el ajuste se bas´o ´unicamente en los datos para PLLA, estudiamos sistemas tambi´en para PDLA (XD= 1) y para copol´ımeros con fracciones intermedias de D-lactida. Esto nos permiti´o verificar cualitativamente que la tendencia predicha por la ecuaci´on se mantiene: Tg disminuye progresivamente al aumentar XD, siendo l´ogico que el valor m´as bajo corresponda al PDLA, un valor intermedio a los copol´ımeros, y el m´as alto al PLLA. Esta tendencia puede observarse cualitativamente tambi´en al comparar por pares los sistemas analizados, como se muestra en la figura S1, en el anexo. Aunque no tenemos suficientes datos para una validaci´on cuantitativa completa, este an´alisis coloca los resultados dentro de un modelo te´orico s´olido, lo que nos brinda una herramienta valiosa para comprender c´omo la composici´on influye en la transici´on v´ıtrea de los copol´ımeros de PLA. En la figura 4.5b, se observa que el copol´ımero con 55 % de D-lactida obtenido no se encuentra dentro del rango que se esperaba seg´un su composici´on. Dado su alto contenido de D-lactida, este sistema deber´ıa presentar una temperatura de transici´on v´ıtrea significativamente m´as baja. Sin embargo, se ubica en una regi´on de la gr´afica correspondiente a sistemas con menor contenido de D-lactida, lo que indica una desviaci´on considerable respecto a la predicci´on del modelo. Esta discrepancia pone en evidencia las limitaciones de los m´etodos utilizados para estimar Tg, y resalta la complejidad de interpretar la transici´on v´ıtrea en sistemas con microestructuras variables o no completamente caracterizadas. Una posible raz´on para la divergencia que observamos entre los valores simulados de Tgy la tendencia te´orica esperada seg´un el porcentaje de unidades D, es la estructura heterog´enea de los copol´ımeros. En los copol´ımeros que tienen un alto porcentaje de unidades D, como copo55, la estructura amorfa puede mostrar una distribuci´on m´as irregular de las cadenas. Esto puede resultar en que las zonas con alto porcentaje de D y L no est´en distribuidas de manera tan uniforme y, por lo tanto, se creen zonas de empaquetamiento heterog´eneo, provocando regiones de diferente movilidad y, por tanto, de transici´on v´ıtrea compleja. Tambi´en es importante tener en cuenta que cuanto m´as heterog´eneo sea el material estructural y din´amicamente, m´as dif´ıcil ser´a la estimaci´on de la transici´on v´ıtrea (v´ease la discusi´on relacionada con la figura 4.3) y, por lo tanto, estas dificultades t´ecnicas tambi´en pueden ser responsables de la desviaci´on. Este comportamiento sugiere que, m´as all´a de una simple correlaci´on con la composici´on, el valor de Tgen estos sistemas tambi´en est´a fuertemente influenciado por factores estructurales como la topolog´ıa de las cadenas y la distribuci´on de los extremos libres. Por lo tanto, los resultados obtenidos subrayan la importancia de considerar la arquitectura molecular completa, y no solo la fracci´on de comon´omero D, al interpretar el comportamiento t´ermico de los copol´ımeros en simulaciones. 4.5. Discusi´on sobre los datos te´oricos y los datos presentes en las fichas t´ecnicas del PLA para impresi´on 3D 41 280 300 320 340 360 380 400 420 0 5 10 15 20 25 30 35 40 45 50 (a) T0, Tg, Tg(exp) [K] Mn [kg/mol] Este trabajo Klonos et al. McAliley et al. Guseva et al. Glova et al. Klajmon et al. 260 280 300 320 340 360 380 400 0 2 4 6 8 (b) T0 [K] Mn [kg/mol] L10 L30 L100 D100 C16 C55 Figura 4.5: Relaci´on entre la temperatura de transici´on v´ıtrea Tgy el peso molecular promedio en n´umero (Mn). (a) Los valores obtenidos en este trabajo para PLLA (T0) se comparan con datos experimentales (Tg(exp)) (PLA con un contenido de L cerca de 95 %) y simulaciones de PLLA (sin extrapolar con WLF, Tg) reportadas en la literatura. Se incluye como referencia la curva original de Flory–Fox (l´ınea continua). (b) Comparaci´on de los valores de T0obtenidos en este estudio (representado por s´ımbolos) con la predicci´on de la ecuaci´on de Flory–Fox para PLLA con 0 % de D-l´actida: la l´ınea discontinua representa la curva original, mientras que la l´ınea continua muestra una versi´on ajustada por nosotros, modificando los par´ametros T∞ gy K para PLLA. 4.5. Discusi´on sobre los datos te´oricos y los datos presentes en las fichas t´ecnicas del PLA para impresi´on 3D En los materiales de PLA que se utilizan para la impresi´on 3D, la temperatura de transici´on v´ıtrea que se menciona en las fichas t´ecnicas puede variar debido a varios factores. Algunos de los 42 Cap´ıtulo 4. Resultados y Discusi´on m´as importantes son el peso molecular y, por ende, la polidispersidad, as´ı como la composici´on del material y la inclusi´on de aditivos. Algunos de estos aspectos han sido abordados en este trabajo. El tipo de PLA com´unmente usado en impresi´on 3D es el PLA3D850. Este material presenta una temperatura de transici´on v´ıtrea de 61.5±2.4◦C, equivalente a 334.7±2.4 K. [41] El PLA comercial es polidisperso, lo que significa que contiene cadenas de diferentes longitudes y, por ende, movilidades segmentarias variadas. Esta diversidad se traduce en una Tg amplia: dentro de una misma muestra, hay regiones que se ablandan a diferentes temperaturas. El pico de Tgresultante es menos claro y su posici´on representa un promedio ponderado entre los extremos, aunque el margen de error asociado rara vez se menciona. Adem´as, las cadenas de alto peso molecular juegan un papel crucial en el comportamiento t´ermico; a medida que aumenta su proporci´on, el Tgobservado tiende a acercarse al valor m´aximo te´orico T∞ g, que corresponde a cadenas infinitamente largas. El PLA que se utiliza en la impresi´on 3D no es un pol´ımero en su forma m´as pura, sino m´as bien una mezcla o copol´ımero que incluye diferentes proporciones de is´omeros L y D, junto con varios aditivos y modificadores. Esta variabilidad en su composici´on estereoqu´ımica tiene un impacto en la Tg, ya que el porcentaje de is´omero D y la presencia de copol´ımeros afectan la movilidad segmentaria y el grado de cristalinidad. Esto provoca que el Tgdel PLA comercial sea diferente de los valores t´ıpicos que se encuentran en el PLA puro. Los plastificantes, que se a˜naden como aditivos al PLA, reducen su temperatura de transici´on v´ıtrea y mejoran su ductilidad y elongaci´on al romperse. Sin embargo, si se utilizan en altas concentraciones, pueden afectar negativamente la resistencia mec´anica, ya que provocan una p´erdida de cohesi´on en la matriz polim´erica [42]. Estudiar sistemas monodispersos con composici´on bien definida facilita la comprensi´on de los fen´omenos que afectan la Tgy proporciona criterios para dise˜nar formulaciones de PLA con ventanas de procesamiento m´as estrechas, transiciones t´ermicas m´as definidas y un rendimiento m´as predecible en la impresi´on 3D. Aunque sintetizar estos sistemas modelo con un control preciso sobre el peso molecular y la composici´on puede ser complicado, las simulaciones ofrecen un control total sobre estas variables. Esto facilita un an´alisis m´as detallado de los factores que afectan la Tgy ayuda en el dise˜no racional de materiales optimizados para aplicaciones espec´ıficas. 5. Conclusiones En este estudio, se llevaron a cabo simulaciones de din´amica molecular para investigar la temperatura de transici´on v´ıtrea (Tg) del ´acido polil´actico (PLA), un pol´ımero de origen natural que ha captado mucho inter´es en aplicaciones como la impresi´on 3D. Se examin´o el impacto de tres factores clave: el peso molecular, la composici´on estereoqu´ımica (la proporci´on de mon´omeros L y D) y la velocidad de enfriamiento. Se simularon seis tipos de sistemas: tres homopol´ımeros de L-lactida (PLLA, con 10, 30 y 100 mon´omeros), un homopol´ımero de D-lactida (PDLA, 100 mon´omeros), y dos copol´ımeros con proporciones L/D de 84/16 % y 45/55 %, ambos con 100 mon´omeros. Para representar adecuadamente el comportamiento del PLA, espec´ıficamente para los copol´ımeros, se modific´o un campo de fuerza desarrollado por McAlileyet al. [16]. La modificaci´on empleada ten´ıa como objetivo mantener la descripci´on original de las rotaciones de la cadena central y, al mismo tiempo, mejorar la viabilidad de la simulaci´on. Cada sistema fue enfriado desde 500 K hasta 200 K a diferentes velocidades, que iban desde muy r´apidas (600 K/ns) hasta m´as lentas (2.5 K/ns). Durante este proceso, se registr´o el cambio de volumen espec´ıfico para identificar el punto de transici´on v´ıtrea. Con estos datos, se aplicaron dos m´etodos para calcular la Tg, y luego aplicando varios modelos te´oricos se compararon los resultados con los datos de la literatura. Tambi´en se analiz´o c´omo la Tgdepende del tama˜no de las cadenas utilizando la ecuaci´on de Flory-Fox [36]. En los homopol´ımeros PLLA, notamos que la Tgaumenta con la longitud de las cadenas dentro del rango que estudiamos. Este comportamiento tambi´en ha sido documentado en otros estudios de din´amica molecular, aunque en muchos casos no se hace una extrapolaci´on adecuada a las velocidades de enfriamiento experimentales. En nuestra investigaci´on, al aplicar esta extrapolaci´on, conseguimos estimaciones de la Tgque se acercan bastante a los valores experimentales que se han reportado en la literatura, lo que refuerza la validez y precisi´on del enfoque que hemos adoptado. En el caso de los copol´ımeros y el homopol´ımero PDLA, los resultados fueron, en general, coherentes con la ecuaci´on de Flory-Fox. Se observ´o una disminuci´on de la Tgal aumentar el 44 Cap´ıtulo 5. Conclusiones contenido de mon´omeros D, como sucedi´o con el copol´ımero que ten´ıa un 16 % de D (copo16D) y en el caso de PDLA. Sin embargo, el copol´ımero con un 55 % de D (copo55D) no sigui´o esta tendencia, mostrando una Tgm´as alta de lo que se esperaba. Esto podr´ıa deberse a una distribuci´on no aleatoria de los mon´omeros y/o a dificultades para estimar el punto de transici´on a partir de los datos de simulaci´on en el caso de sistemas din´amicamente heterog´eneos, como el de copo55D. Los resultados que hemos obtenido destacan la conexi´on tan estrecha entre la estructura estereoqu´ımica del PLA y su temperatura de transici´on v´ıtrea. Al cambiar la proporci´on de los mon´omeros L y D, se modifica la rigidez y el empaquetamiento de las cadenas, lo que impacta directamente en su comportamiento t´ermico. Este estudio resalta la importancia de los an´alisis computacionales, especialmente cuando se aplican de manera sistem´atica a sistemas modelo bien definidos. Las simulaciones de din´amica molecular nos permiten observar con gran detalle c´omo evoluciona la estructura de materiales complejos, como los pol´ımeros biodegradables. En nuestro caso, han sido clave para profundizar en la comprensi´on del PLA, logrando estimaciones de Tgque se confrontaron directamente con los datos te´oricos y experimentales y de esta manera sirvieron como puente entre estos dos enfoques. Este m´etodo computacional no solo enriquece nuestra comprensi´on te´orica, sino que tambi´en proporciona herramientas valiosas para el dise˜no racional de materiales con aplicaciones pr´acticas en la industria. 6. Perspectivas de futuro Los pol´ımeros biodegradables est´an ganando terreno como una alternativa fundamental frente a los pl´asticos convencionales, gracias a su origen renovable y su menor impacto en el medio ambiente. Sin embargo, para que puedan ser utilizados a gran escala, todav´ıa es necesario superar ciertos desaf´ıos relacionados con sus propiedades t´ermicas, mec´anicas y de procesabilidad. En este sentido, las t´ecnicas computacionales, como la din´amica molecular, se est´an consolidando como herramientas sostenibles y efectivas para investigar la relaci´on entre la estructura y las propiedades de estos materiales. Esto permite optimizar el dise˜no de nuevos pol´ımeros sin depender ´unicamente de ensayos experimentales. Al combinar estas metodolog´ıas con tecnolog´ıas como la impresi´on 3D, se abre la puerta a una fabricaci´on m´as inteligente y respetuosa con el medio ambiente. De cara al futuro, ampliar los estudios a copol´ımeros con composiciones m´as diversas podr´ıa ayudar a validar modelos te´oricos y a obtener par´ametros l´ımite como el T∞ g. Adem´as, contar con modelos validados nos brinda la oportunidad de explorar otras propiedades importantes, como la viscosidad, la difusividad o el comportamiento bajo deformaci´on, lo que ampliar´a nuestro conocimiento para aplicaciones industriales avanzadas. La investigaci´on en pol´ımeros biodegradables se est´a orientando hacia un enfoque interdisciplinario, donde la simulaci´on y la sostenibilidad jugar´an un papel clave en el desarrollo de materiales funcionales que respondan a los desaf´ıos del futuro. 46 Cap´ıtulo 6. Perspectivas de futuro Bibliograf´ıa [1] M. Rubinstein and R. H. Colby, Polymer physics. Oxford university press, 2003. [2] M. Cortizo, T. Oberti, and P. Peruzzo, Introducci´on a la s´ıntesis de pol´ımeros. Editorial de la Universidad Nacional de La Plata (EDULP), 2023. [3] J. E. Mark et al.,Physical properties of polymers handbook, vol. 1076. Springer, 2007. [4] D. I. Bower, An introduction to polymer physics. Cambridge University Press, 2002. [5] A. Samir, F. H. Ashour, A. A. Hakim, and M. Bassyouni, “Recent advances in biodegradable polymers for sustainable applications,” Npj Materials Degradation, vol. 6, no. 1, p. 68, 2022. [6] M. S. Kim, H. Chang, L. Zheng, Q. Yan, B. F. Pfleger, J. Klier, K. Nelson, E. L.-W. Majumder, and G. W. Huber, “A review of biodegradable plastics: chemistry, applications, properties, and future research needs,” Chemical Reviews, vol. 123, no. 16, pp. 9915–9939, 2023. [7] N.-A. A. B. Taib, M. R. Rahman, et al., “A review on poly lactic acid (pla) as a biodegradable polymer,” Polymer Bulletin, vol. 80, no. 2, pp. 1179–1213, 2023. [8] T. Casalini, F. Rossi, A. Castrovinci, and G. Perale, “A perspective on polylactic acid-based polymers use for nanoparticles synthesis and applications,” Frontiers in bioengineering and biotechnology, vol. 7, no. 259, 2019. [9] P. A. Klonos, N. D. Bikiaris, P. Barmpalexis, and A. Kyritsis, “Segmental mobility in linear polylactides of various molecular weights,” Polymer, vol. 305, p. 127177, 2024. [10] A. AL-Zaidi and F. Al-Gawhari, “Types of polymers using in 3D printing and their applications: a brief review,” EJTAS, vol. 1, pp. 978–85, 2023. [11] M. N. Andanje, J. W. Mwangi, B. R. Mose, and S. Carrara, “Biocompatible and biodegradable 3D printing from bioplastics: A review,” Polymers, vol. 15, no. 10, p. 2355, 2023. 48 Bibliograf´ıa [12] M. P. Allen et al., “Introduction to molecular dynamics simulation,” Computational soft matter: from synthetic polymers to proteins, vol. 23, no. 1, pp. 1–28, 2004. [13] W. G. Noid, “Perspective: Coarse-grained models for biomolecular systems,” The Journal of chemical physics, vol. 139, no. 9, 2013. [14] A. Soldera and N. Metatla, “Glass transition of polymers: Atomistic simulation versus experiments,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, vol. 74, no. 6, p. 061803, 2006. [15] M. Klajmon, V. Aulich, J. Ludik, and C. Cervinka, “Glass transition and structure of organic polymers from all-atom molecular simulations,” Industrial & Engineering Chemistry Research, vol. 62, no. 49, p. 2143721448, 2023. [16] J. H. McAliley and D. A. Bruce, “Development of force field parameters for molecular simulation of polylactide,” Journal of chemical theory and computation, vol. 7, no. 11, pp. 3756–3767, 2011. [17] D. V. Guseva, A. A. Lazutin, and V. V. Vasilevskaya, “Atomistic simulation of poly (lactic acid) of different regioregularity,” Polymer, vol. 221, p. 123577, 2021. [18] E. Christofi, P. Bacova, and V. A. Harmandaris, “Physics-informed deep learning approach for reintroducing atomic detail in coarse-grained configurations of multiple poly (lactic acid) stereoisomers,” Journal of Chemical Information and Modeling, vol. 64, no. 6, pp. 1853– 1867, 2024. [19] V. A. Harmandaris and V. G. Mavrantzas, “Molecular dynamics simulations of polymers,” in Simulation methods for polymers, pp. 178–218, CRC Press, 2004. [20] D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications. Elsevier, 2023. [21] J. Li, “Basic molecular dynamics,” in Handbook of Materials Modeling: Methods, pp. 565– 588, Springer, 2005. [22] H. Kumar and P. K. Maiti, Introduction to molecular dynamics simulation, pp. 161–197. Springer, 2011. [23] W. F. Van Gunsteren and H. J. Berendsen, “A leap-frog algorithm for stochastic dynamics,” Molecular Simulation, vol. 1, no. 3, pp. 173–185, 1988. Bibliograf´ıa 49 [24] S. Nos´e, “A molecular dynamics method for simulations in the canonical ensemble,” Molecular physics, vol. 100, no. 1, pp. 191–198, 2002. [25] H. J. Berendsen, J. v. Postma, W. F. Van Gunsteren, A. DiNola, and J. R. Haak, “Molecular dynamics with coupling to an external bath,” The Journal of chemical physics, vol. 81, no. 8, pp. 3684–3690, 1984. [26] M. J. Abraham, T. Murtola, R. Schulz, S. P´all, J. C. Smith, B. Hess, and E. Lindahl, “Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX, vol. 1, pp. 19–25, 2015. [27] C. Levinthal, “Molecular model-building by computer,” Scientific american, vol. 214, no. 6, pp. 42–53, 1966. [28] M. P. Oliveira, Y. M. Gon¸calves, S. K. Ol Gheta, S. R. Rieder, B. A. Horta, and P. H. Hunenberger, “Comparison of the united-and all-atom representations of (halo) alkanes based on two condensed-phase force fields optimized against the same experimental data set,” Journal of chemical theory and computation, vol. 18, no. 11, pp. 6757–6778, 2022. [29] M. A. Gonz´alez, “Force fields and molecular dynamics simulations,” ´ Ecole th´ematique de la Soci´et´e Fran¸caise de la Neutronique, vol. 12, pp. 169–200, 2011. [30] W. Jorgensen, D. Maxwell, and J. Tirado-Rives, “Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids,” Journal of the american chemical society, vol. 118, no. 45, pp. 11225–11236, 1996. [31] B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. a. Swaminathan, and M. Karplus, “CHARMM: a program for macromolecular energy, minimization, and dynamics calculations,” Journal of computational chemistry, vol. 4, no. 2, pp. 187–217, 1983. [32] M. S. Badar, S. Shamsi, J. Ahmed, and M. A. Alam, “Molecular dynamics simulations: concept, methods, and applications,” in Transdisciplinarity, pp. 131–151, Springer, 2022. [33] B. Hess, H. Bekker, H. J. Berendsen, and J. G. Fraaije, “LINCS: A linear constraint solver for molecular simulations,” Journal of computational chemistry, vol. 18, no. 12, pp. 1463– 1472, 1997. [34] M. L. Williams, R. F. Landel, and J. D. Ferry, “The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids,” Journal of the American Chemical society, vol. 77, no. 14, pp. 3701–3707, 1955.