Master thesis: Designing a bottom-up computational approach for sustainable polymers

Dise˜no de un enfoque computacional ascendente para los pol´ımeros sostenibles Gonzalo Gonz´alez Huarte Tutora: Petra Baˇcov´a Trabajo Fin de M´aster en Nanociencia y Tecnolog´ıa de Materiales Facultad de Ciencias, Puerto Real (C´adiz), curso 2023/2024 Resumen Este trabajo presenta un estudio detallado de la policaprolactona (PCL), un pol´ımero biodegradable con aplicaciones potenciales en biomedicina y fabricaci´on aditiva. El objetivo principal es adaptar metodolog´ıas computacionales para analizar la relaci´on entre la estructura, propiedades y rendimiento de este pol´ımero, utilizando simulaciones de din´amica molecular. El enfoque computacional ascendente, basado en simulaciones atom´ısticas, permite entender c´omo las interacciones a nivel molecular influyen en las propiedades macrosc´opicas. Se detalla el desarrollo de un campo de fuerzas optimizado para estudiar la PCL en distintos entornos (vac´ıo, disoluci´on y fundido), con especial ´enfasis en el comportamiento estructural y din´amico de las cadenas polim´ericas. Para entender este comportamiento, se han calculado dos propiedades estructurales fundamentales en pol´ımeros: el radio de giro, Rgy la distancia entre extremos, Re. Tambi´en se han obtenido resultados cuantitativos de algunas propiedades claves para el desarrollo de los futuros modelos de los sistemas en el fundido que son: la densidad y el tiempo de relajaci´on. Los datos han sido comparados y verificados con los datos publicados en la literatura. El comportamiento de las cadenas de distinto peso molecular est´a de acuerdo con los modelos te´oricos. Este estudio pretende establecer una base para futuros trabajos que busquen optimizar el uso de biopol´ımeros en aplicaciones sostenibles, como la impresi´on 3D, mejorando su rendimiento y compatibilidad con procesos industriales. ii Abstract This work presents a detailed study of polycaprolactone (PCL), a biodegradable polymer with potential applications in biomedicine and additive manufacturing. The main objective is to adapt computational methodologies to analyse the relationship between the structure, properties and performance of this polymer, using molecular dynamics simulations. The bottom-up approach, based on atomistic simulations, allows us to understand the influence of the interactions at the molecular level on macroscopic properties. An optimised force field to study PCL in different environments (vacuum, dissolution and melt) is developed, with special emphasis on the structural and dynamic behaviour of the studied polymer systems. To understand this behaviour, two fundamental structural properties of polymers have been calculated: the radius of gyration, Rgand the end-to-end distance, Re. Quantitative results have also been obtained for some properties, essential for the development of future models of the systems in melt: the density and the terminal relaxation time. The data have been compared and verified with data published in literature. The behaviour of the chains with different molecular weights is in agreement with the theoretical models. This study aims to establish a basis for future work, in order to optimise the use of biopolymers in sustainable applications, such as 3D printing, by improving their performance and compatibility with industrial processes. iv ´ Indice general 1. Introducci´on 1 1.1. Materialespolim´ericos ................................. 1 1.2. Pol´ımeros biodegradables y pr´acticas sostenibles . . . . . . . . . . . . . . . . . . 2 1.2.1. Fabricaci´onaditiva............................... 4 1.2.2. Policaprolactona ................................ 4 1.3. T´ecnicas computacionales en la ciencia de los pol´ımeros . . . . . . . . . . . . . . 5 1.3.1. Simulaciones atom´ısticas con PCL: estado del arte . . . . . . . . . . . . . 7 2. Objetivos 9 3. Materiales y M´etodos 11 3.1. Din´amicaMolecular .................................. 11 3.1.1. Ecuaciones de movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.2. Algoritmo.................................... 12 3.1.3. Sistemas estad´ısticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.4. Condiciones de contorno y convenci´on de im´agenes m´ınimas . . . . . . . . 15 3.2. Sistemamodelo..................................... 17 3.2.1. Campodefuerzas ............................... 18 3.2.2. Interacciones enlazantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.3. Interacciones no enlazantes . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.4. Preparaci´on de las muestras y detalles de simulaci´on . . . . . . . . . . . . 23 4. Resultados y Discusi´on 27 4.1. Desarrollo del campo de fuerzas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2. Propiedades estructurales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2.1. Inspecciones visuales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2.2. Cuantificaci´on de las inspecciones visuales . . . . . . . . . . . . . . . . . . 30 4.2.3. Evoluci´on de las propiedades estructurales en el tiempo . . . . . . . . . . 32 vi ´ Indice general 4.2.4. Evoluci´on de las propiedades estructurales frente a N . . . . . . . . . . . . 33 4.2.5. Comparaci´on del valor de Rg/Reexperimental con el te´orico . . . . . . . 37 4.3. Propiedadesdin´amicas................................. 38 5. Conclusiones 43 6. Perspectivas de futuro 45 Bibliograf´ıa 47 Ap´endice Anexo 51 1. Introducci´on En cuanto a pol´ımeros se refiere, uno de los aspectos m´as importantes a considerar es encontrar un equilibrio entre la relaci´on estructura-propiedades-comportamiento. Ser capaz de conferir una determinada estructura al pol´ımero conlleva la aparici´on de unas propiedades espec´ıficas que pueden ser la clave para desempe˜nar la funci´on deseada. El objetivo final es comprender en su totalidad dicha relaci´on, para ello, introduciremos en esta secci´on diversos conceptos relacionados primero con la estructura de los materiales polim´ericos (Sec. 1.1), luego acerca de sus propiedades y sobre su comportamiento en servicio esperado desde un punto de vista de aplicaci´on industrial (Sec 1.2). En el presente trabajo, el objetivo principal es elucidar algunos aspectos acerca de esta relaci´on a trav´es de m´etodos computacionales, que ser´an introducidos en la Sec 1.3. 1.1. Materiales polim´ericos El t´ermino pol´ımero es utilizado frecuentemente para referirnos a una clase de macromol´eculas que consisten en un conjunto de unidades qu´ımicas, regularmente repetidas, unidas covalentemente entre s´ı formando una cadena. Esa unidad qu´ımica repetida, recibe el nombre de mon´omero. Si solo existe un tipo de mon´omero se llamar´a homopol´ımero, en cambio, si hay dos o m´as tipos de mon´omeros, se llamar´a heteropol´ımero. Normalmente, las cadenas polim´ericas se suelen representar de la siguiente manera: −(CxHy)N−, siendo Nla cantidad de mon´omeros que alberga la cadena, tambi´en llamado el grado de polimerizaci´on. Si una cadena polim´erica contiene un n´umero reducido de mon´omeros (generalmente menos de 20), se dice que es un olig´omero. La longitud de la cadena, y en consecuencia su peso molecular, es uno de los factores claves a la hora de justificar las propiedades de los pol´ımeros [1]. Esta estructura polim´erica se genera mediante un proceso llamado polimerizaci´on, donde ocurre la uni´on de mon´omeros mediante enlaces covalentes. La masa molar, M, del pol´ımero es directamente proporcional al n´umero de mon´omeros presentes en la cadena N: M=NMmon,(1.1) 8Cap´ıtulo 1. Introducci´on polim´erica de PCL de longitud variable con una membrana de fosfol´ıpido. Al hacer una comparativa de los valores de radio de giro, Rg, entre agua (un disolvente malo) y la membrana POPC (1-palmitoil-2-oleoil-sn-glicero-3-fosfocolina), han llegado a la conclusi´on que el valor de Rges superior en POPC. Las cadenas en POPC tienen configuraciones extendidas y adem´as, los valores medios de Rgincrementan al aumentar el n´umero de mon´omeros en la cadena polim´erica [14]. Yungerman et al. utilizaron dos tipos de los campos de fuerzas, OPLS [15] y PCFF [16], para investigar las propiedades de la PCL usando diacrilatos como reactivos funcionales en los grupos terminales. [17] De manera similar al art´ıculo [14], tambi´en han optado por las simulaciones de din´amica molecular con todos los ´atomos. Sus resultados fueron comparados con datos experimentales conocidos y asunciones te´oricas para verificar su aplicabilidad y de tal manera encontrar el campo de fuerza adecuado para la descripci´on at´omica de PCL. En sus resultados, se pudo concluir que es indiferente el uso de OPLS o PCFF ya que no presenta una diferencia estructural significativa. Es decir, no hay variaci´on en el tama˜no de las cadenas caracterizado por Rgen los sistemas diluidos ni en fundido [17]. Dada la versatilidad de OPLS y su uso en muchas simulaciones de biomol´eculas, este hecho abre v´ıas para emplear OPLS para investigar las interacciones de pol´ımeros como PLC con una gran variedad de sistemas moleculares. Di Pasquale et al. publicaron dos trabajos consecutivos [18, 19] sobre el comportamiento de PCL en diferentes disolventes, con el objetivo de observar la precipitaci´on de nanopart´ıculas polim´ericas. Para estudiar la precipitaci´on en reactores de chorro de impacto confinado (CIJRs, en ingl´es), en [18] han combinado el modelo computacional de fluido din´amico (CFD) con las simulaciones de din´amica molecular. Desde el punto de vista de inter´es del presente trabajo, se observaron los valores estimados para el radio de giro en dos disolventes: agua y acetona. A la vista de sus resultados obtenidos, es posible llegar a la conclusi´on de que existe una influencia notoria en el tipo de disolvente empleado, ya que los valores de Rgen acetona son significativamente mayores que en agua. Adem´as, se observa una relaci´on directamente proporcional entre Rgy el n´umero de mon´omeros de la cadena [18]. Yihan Nie elabor´o una tesis doctoral con el fin de estudiar la viabilidad del uso de PCL en la impresi´on 3D a˜nadiendo nanorellenos hechos de grafeno para potenciar sus propiedades mec´anicas durante la impresi´on [11]. Para lograrlo, estudi´o la reolog´ıa de la PCL y de sus compuestos de grafeno usando una simulaci´on de din´amica molecular fuera del equilibrio, para demostrar la influencia del tama˜no de los nanorellenos en la viscosidad del compuesto. Manteniendo Rgconstante, la viscosidad del pol´ımero variar´a de diferentes maneras seg´un si el tama˜no del nanorelleno sea superior o inferior a dicho valor [11]. 2. Objetivos En el ´ambito de los pol´ımeros biodegradables, la policaprolactona (PCL) destaca como uno de los materiales m´as populares debido a su amplia aplicabilidad en la fabricaci´on aditiva y en el campo de la biomedicina, concretamente en el uso de pr´otesis. Sin embargo, aunque actualmente se encuentra disponible comercialmente, la sustituci´on de los pol´ımeros sint´eticos por este material sostenible no es del todo viable por sus deficientes propiedades mec´anicas y su limitada estabilidad t´ermica. El n´ucleo de este estudio lo constituyen las t´ecnicas computacionales, resultando ser una v´ıa de estudio m´as respetuosa con el medio ambiente en comparaci´on con los m´etodos tradicionales de experimentaci´on. El objetivo principal de este trabajo es adaptar y extender las metodolog´ıas de simulaci´on com´unmente aplicadas a pol´ımeros sint´eticos, para abordar la complejidad inherente de los pol´ımeros biodegradables, en este caso de la policaprolactona. Con el fin de abordar los objetivos propuestos, relacionados con la correlaci´on entre estructura, propiedades y rendimiento de este pol´ımero, se propone un estudio computacional sistem´atico ascendente. El estudio est´a enfocado en la descripci´on detallada de la estructura de pol´ımeros biodegradables, empezando por el desarrollo de un campo de fuerzas, seguido por la verificaci´on del modelo y la creaci´on de nuevas simulaciones atom´ısticas detalladas. Ya que el objetivo final es extender el modelo presentado a modelos experimentalmente relevantes y de tal manera dilucidar el comportamiento de la PCL durante la fabricaci´on aditiva, el trabajo se centra en los sistemas en el fundido y en calcular y recolectar los datos claves para construir una base s´olida de los futuros trabajos computacionales. 10 Cap´ıtulo 2. Objetivos 3. Materiales y M´etodos 3.1. Din´amica Molecular Las simulaciones de Din´amica Molecular (Molecular Dynamics o MD, en ingl´es) son un m´etodo computacional capaz de predecir c´omo evoluciona el sistema en el tiempo. En un sistema de este tipo, el movimiento del sistema es cl´asico, es decir, sigue las ecuaciones de movimiento de Newton. Esta predicci´on nos ayuda a comprender las propiedades de los sistemas moleculares en funci´on de su estructura y teniendo en cuenta las interacciones microsc´opicas entre ellas [6]. 3.1.1. Ecuaciones de movimiento Nuestro sistema simulado consiste en un n´umero determinado de ´atomos y/o mol´eculas denotado como n. Estas part´ıculas interact´uan entre s´ı mediante un campo de fuerzas, tambi´en llamado potencial interat´omico. El objetivo de toda simulaci´on MD es determinar la posici´on y la velocidad de todas y cada una de las part´ıculas en un momento dado. Dado que las part´ıculas interact´uan entre s´ı, su movimiento est´a correlacionado y, por tanto, este problema a resolver se denomina problema de ncuerpos. Para ello hay que resolver, paso a paso, las ecuaciones del movimiento. El movimiento del sistema a trav´es de la fase se rige por las ecuaciones de movimiento hamiltonianas, donde el hamiltoniano (3.1), H, representa la energ´ıa del sistema [6]. H(p,r) = K(p) + U(r).(3.1) Hamiltoniano tiene dos componentes, la energ´ıa cin´etica K(˙ r) y la energ´ıa potencial, escrita como U(r). Las ecuaciones de movimiento para npart´ıculas indicadas como i= 1,2, ..., n, se encuentran en las ecuaciones 3.2. ˙ ri=pi mi =vi˙ pi=−∇riU(r),(3.2) donde rirepresenta la posici´on cartesiana y piel momento de cada part´ıcula i,˙ rirepresenta la derivada temporal de riy representa la velocidad de la part´ıcula (vi), y midenota su masa. 12 Cap´ıtulo 3. Materiales y M´etodos Como resultado de los ajustes realizados, el movimiento de cada part´ıcula del sistema sigue una trayectoria determinada. La trayectoria es la evoluci´on de cada posici´on y velocidad (momento) del sistema a lo largo del tiempo de la simulaci´on y representa la salida principal en la simulaci´on. Esta informaci´on podr´ıa obtenerse resolviendo las ecuaciones de movimiento en cada paso de la simulaci´on, lo cual es anal´ıticamente dif´ıcil y tedioso de resolver. Por eso se ha utilizado un integrador, un operador matem´atico que ayuda a resolver estas ecuaciones utilizando diferencias finitas con un paso de tiempo que se ha denotado como ∆t. Es una aproximaci´on a una serie de Taylor. Los integradores m´as famosos son el algoritmo de Verlet y el Leap Frog [20]. Para comprender el comportamiento de todo el sistema y poder medir sus propiedades empleando la trayectoria obtenida, se ha usado el concepto de ergodicidad de la mec´anica estad´ıstica. La ergodicidad puede entenderse como una propiedad del sistema. Se basa en la hip´otesis de que si se observa la evoluci´on de un objeto individual de un sistema a lo largo del tiempo, al cabo de cierto tiempo, se puede ver todos los estados estad´ısticamente posibles de ese objeto individual, y que bajo esta hip´otesis esa observaci´on puede ser representativa de todo un conjunto, es decir, de m´ultiples objetos del sistema. La hip´otesis puede entenderse pr´acticamente de tal manera que la media a largo plazo es igual a la media del conjunto en el l´ımite de un n´umero infinito de miembros en el conjunto [20]. 3.1.2. Algoritmo Existen muchos tipos diferentes de simulaciones MD, sin embargo, cada una de ellas sigue un conjunto sencillo de pasos, que se puede denominar algoritmo. Los pasos a realizar en una simulaci´on MD t´ıpica son: 1. Seleccionar un modelo de interacci´on. 2. Seleccionar condiciones de contorno. 3. Seleccionar condiciones iniciales. 4. Seleccionar sistema estad´ıstico (NVE, NVT, NPT...). 5. Llevar a cabo la simulaci´on hasta que se alcance el equilibrio (dependiente de la propiedad). 6. Realizar una simulaci´on de producci´on para obtener medias termodin´amicas, posici´on y velocidades. 7. Analizar los resultados por post-procesado. 3.1. Din´amica Molecular 13 El paso 1 se explica en detalle en la Sec. 3.2.1, el paso 2 se menciona en la Sec. 3.1.4, los pasos 3 y 4 se desarrollan en la Sec. 3.1.3. Los pasos 5-7 son el n´ucleo principal de la simulaci´on y se explican en detalle en el cap´ıtulo 4. Es necesario comprender que el paso 5 no s´olo depende de las propiedades, sino tambi´en del sistema. Por ello, depende de otras condiciones experimentales como la temperatura [20]. El algoritmo que se utilizar´a para configurar la simulaci´on se muestra en la Figura 3.1 Posiciones iniciales y velocidades Calcular fuerzas Actualizar posiciones y velocidades Calcular y almacenar Cantidades básicas Mostrar resultados (posición, temperatura...) Bucle en los time steps Figura 3.1: Ilustraci´on esquem´atica de los pasos seguidos en una simulaci´on din´amica molecular t´ıpica. El bucle se repite muchas veces, m´as concretamente, tantas veces como el n´umero de pasos temporales (time steps) seleccionados para la simulaci´on. La selecci´on depende de factores como la propiedad investigada y/o el tipo de sistema de simulaci´on. Para fijar las posiciones iniciales, se suele partir de las coordenadas con la distribuci´on espacial de los ´atomos dada por el campo de fuerza seleccionado, que se ha establecido previamente mediante c´alculos cu´anticos y se ha verificado mediante mediciones experimentales. Para inicializar las velocidades, se suelen generar ajust´andolas a una distribuci´on de MaxwellBoltzmann [6] dada por la ecuaci´on 3.3. P(v) = m 2πkBT3 2 e −mv2 2kBT,(3.3) 14 Cap´ıtulo 3. Materiales y M´etodos donde ves la velocidad de la part´ıcula, mrepresenta su masa, kBes la constante de Boltzmann y Tla temperatura del sistema. Esta elecci´on ser´ıa muy conveniente sobre todo en los sistemas con temperatura constante (v´ease la secci´on siguiente), ya que asegura que la simulaci´on comienza con la temperatura deseada. 3.1.3. Sistemas estad´ısticos En el cuarto paso del algoritmo empleado es necesario seleccionar un sistema estad´ıstico, que es un t´ermino utilizado en mec´anica estad´ıstica para referirse a un conjunto espec´ıfico de estados de un sistema en condiciones seleccionadas. Un sistema termodin´amico proporciona una forma de deducir las propiedades termodin´amicas de un sistema a trav´es de las leyes de la mec´anica cl´asica y cu´antica. [20] La idea principal es que diferentes sistemas estad´ısticos representan sistemas en equilibrio con diferentes caracter´ısticas estad´ısticas particulares. El sistema estad´ıstico que se elija depende del problema espec´ıfico que se intente resolver y de las condiciones macrosc´opicas en las que se quiera simular el sistema [20]. En otras palabras, establecer las restricciones macrosc´opicas que definen cada sistema estad´ıstico es una forma de imitar las condiciones experimentales en la simulaci´on microsc´opica. Existen determinados sistemas estad´ısticos caracterizados por variables constantes. Las m´as comunes son: n´umero de part´ıculas (N, o en el presente trabajo n), volumen (V), energ´ıa (E), presi´on (P) y temperatura (T). De acuerdo con esto, los sistemas de mec´anica estad´ıstica m´as comunes se caracterizan como: NVE: sistema microcan´onico NVT: sistema can´onico NPT: sistema isot´ermico-isob´arico El sistema microcan´onico no es muy pr´actico, ya que requiere mantener constante la energ´ıa total. Por otro lado, los conjuntos isot´ermico e isob´arico reflejan condiciones experimentales t´ıpicas. En este caso se utilizan termostato y barostato para controlar la temperatura y la presi´on, respectivamente. Ambos no son m´as que m´etodos num´ericos aplicados para modificar las ecuaciones de movimiento. Uno de los termostatos m´as famosos es el termostato Nos´e-Hoover, que introduce un nuevo t´ermino (grado de libertad, s) en el Lagrangiano. Este nuevo par´ametro desempe˜na el papel de un ba˜no de calor para amortiguar las desviaciones de temperatura con respecto al nivel deseado. Un punto importante en este an´alisis es que el conjunto de ecuaciones de movimiento 3.1. Din´amica Molecular 15 es ´unico, lo que significa que ninguna otra ecuaci´on de la misma forma puede conducir a una distribuci´on can´onica. Esto hace que el termostato Nos´e-Hoover sea apto para la producci´on que se ejecuta [21]. Otra opci´on para series de producci´on es el barostato de Parrinello-Rahman, que permite cambios no s´olo en el tama˜no, sino tambi´en en la forma de la caja de simulaci´on. Esto es interesante en la simulaci´on de s´olidos tales como pol´ımeros cristalinos, ya que permite cambios de fase en la simulaci´on que implican cambios en las dimensiones y ´angulos de la celda unitaria [22]. Otro muy conocido es el termostato y barostato de Berendsen. Berendsen propuso una forma m´as sencilla de realizar simulaciones MD isot´ermicas y/o isob´aricas sin un Lagrangiano extendido. Acopla el sistema a un ba˜no externo de calor y/o presi´on con el grado de acoplamiento definido por una constante de acoplamiento. El m´etodo es m´as f´acil de programar pero es impreciso debido a que la densidad de probabilidad del espacio de fase descrito no se ajusta a un sistema estad´ıstico espec´ıfico. Sin embargo, este m´etodo es muy eficaz para la equilibraci´on del sistema [23]. 3.1.4. Condiciones de contorno y convenci´on de im´agenes m´ınimas En el mundo real, un sistema est´a compuesto por un n´umero incontable de ´atomos del orden del n´umero de Avogadro: NAv = 6,023 ×1023. Simular un n´umero tan elevado de part´ıculas es muy costoso desde el punto de vista computacional. Por lo tanto, para imitar el comportamiento de un sistema m´as grande, se aproxima a un sistema m´as peque˜no y se aplican las condiciones peri´odicas de contorno. Para ello, se divide todo el espacio del sistema en peque˜nas cajas y se estudia una de ellas, que ser´a la “caja de simulaci´on” que se usar´a. Las dem´as cajas de alrededor no son m´as que copias peri´odicas de la caja de simulaci´on elegida. Para entender c´omo funciona esta aproximaci´on, se ilustra en la Figura 3.2. Las condiciones peri´odicas de contorno (P eriodic Boundary Conditions o PBC, en ingl´es) definen las posiciones de las part´ıculas en un sistema con n´umero constante de part´ıculas de tal forma que ninguna de ellas puede salir de la caja. Si una part´ıcula sale de la caja en una direcci´on, su imagen peri´odica realiza el mismo movimiento pero entrando en la caja. As´ı que de esta manera, la caja siempre contiene el mismo n´umero de part´ıculas [6]. La aplicaci´on del PBC tambi´en provoca algunos problemas computacionales. Al estudiar las interacciones de una part´ıcula en la caja de simulaci´on elegida con el resto de part´ıculas, es posible que una interacci´on con una part´ıcula y al mismo tiempo con su imagen peri´odica est´en ocurriendo al mismo tiempo. En este caso, se producir´ıa un artefacto llamado doble conteo. Por eso se adopta la convenci´on de im´agenes m´ınimas, de modo que cada part´ıcula interact´ua con 16 Cap´ıtulo 3. Materiales y M´etodos Figura 3.2: Condiciones de contorno peri´odicas. La caja de simulaci´on empleada en este trabajo est´a representada en magenta, las l´ıneas discontinuas son las cajas circundantes que son s´olo copias de la caja de simulaci´on elegida. la part´ıcula o imagen m´as cercana en la matriz peri´odica. La distancia entre las part´ıculas en la direcci´on dada no puede ser mayor que la mitad de la caja [6]. Figura 3.3: Ilustraci´on esquem´atica de la convenci´on de im´agenes m´ınimas. El recuadro magenta representa la caja de simulaci´on elegida, mientras que el recuadro gris representa la convenci´on de imagen m´ınima elegida. 3.2. Sistema modelo 17 La simulaci´on MD es un m´etodo computacional, y es necesario aplicar t´ecnicas que reduzcan el coste computacional sin perder precisi´on. Los c´alculos m´as costosos de la simulaci´on son los c´alculos de fuerzas, ya que tienen en cuenta cada interacci´on entre part´ıculas. Por ello, para disminuir este coste, se define una distancia de corte, tambi´en llamada radio de corte (rc) [22]. El c´alculo del potencial de pares s´olo se tendr´a en cuenta para aquellas part´ıculas cuya distancia mutua sea menor que el radio definido. Teniendo en cuenta la convenci´on de im´agenes m´ınimas explicada anteriormente, la distancia de corte no puede ser mayor que la mitad de la caja y su valor real viene dado por el campo de fuerza seleccionado. r c Figura 3.4: Ilustraci´on de la distancia de corte. El radio de corte, rc, est´a representado por el c´ırculo rojo y en magenta se ilustra la caja de simulaci´on. 3.2. Sistema modelo En este trabajo se han empleado simulaciones atom´ısticas ya que permiten obtener resultados cuantitativos. En cuanto a las escalas de longitud y tiempo, las simulaciones atom´ısticas operan a escalas de tiempo del orden de los 10−12 hasta 10−6segundos y a escalas de longitud de 10−9 hasta 10−7metros (v´ease Figura 1.3). En otras palabras, los procesos que se estudian son en su mayor´ıa a escala nanom´etrica e incluyen, entre otros, procesos como el enlace de hidr´ogeno, la relajaci´on segmentaria, etc. Para la mol´ecula estudiada de policaprolactona (PCL), se considera el modelo que contiene hidr´ogenos expl´ıcitos, llamado all−atom en ingl´es. Los sistemas fueron simulados por el paquete de simulaci´on de c´odigo abierto Gromacs [24] y por lo tanto en lo que sigue se usar´an las formas 24 Cap´ıtulo 3. Materiales y M´etodos mon´omero central para conseguir la estructura totalmente extendida con el n´umero deseado de mon´omeros. El campo de fuerzas del tr´ımero se dise˜n´o a mano teniendo en cuenta la conectividad de los ´atomos entre s´ı, as´ı como sus ´angulos. Seguidamente, se cre´o la topolog´ıa con la matriz de conectividad para el pol´ımero estudiado, es decir, la lista de los ´atomos teniendo en cuenta las interacciones de 2 cuerpos, 3 cuerpos y 4 cuerpos, definiendo tambi´en el tipo de cada interacci´on. El dise˜no de las simulaciones se realiz´o de la siguiente manera: primero se simularon las cadenas en vac´ıo, que es un disolvente muy malo, para causar el colapso de la cadena y poder colocarla en una caja de simulaci´on m´as peque˜na que la que corresponder´ıa a la cadena totalmente extendida. Luego se prepararon los sistemas en disoluci´on y al final los sistemas en el fundido. A continuaci´on, se resumen los detalles de la preparaci´on de cada sistema: Vac´ıo: Una cadena polim´erica totalmente extendida se coloca en una caja de simulaci´on. Por tener la extensi´on m´axima, la cadena es muy larga y requiere un gran tama˜no de dicha caja. La cadena se centr´o en una caja de simulaci´on lo suficientemente grande para comprender la mol´ecula completa sin mol´eculas de disolvente. Se aplic´o una din´amica estoc´astica para evitar un estado inm´ovil de la cadena polim´erica y un mayor radio de corte rcde 2 nm [28]. Despu´es de cierto tiempo la cadena adquiere un estado de bobina (coil, en ingl´es). El tiempo de la simulaci´on se puede observar en el Cuadro 3.1. La simulaci´on se detuvo poco despu´es de obtener ese estado colapsado mencionado anteriormente. Disoluci´on: Se us´o el resultado obtenido de la simulaci´on en el vac´ıo y se coloc´o esa mol´ecula en el centro de una caja de simulaci´on lo suficientemente grande, ajustando su tama˜no para albergar el n´umero deseado de mol´eculas de disolvente. En este caso se trata de agua seg´un el modelo SPC [29] usado en el art´ıculo de Di Pasquale et al. [18]. El n´umero de mol´eculas de agua se ajust´o para obtener una fracci´on de peso de PCL del 6 %. La fracci´on de peso de cada sistema en disoluci´on es la misma. Se realizaron dos simulaciones, la primera con un paso de tiempo de 0.5 fs debido a algunos solapamientos y situaciones inestables mientras se ejecutaba la simulaci´on, la segunda con el habitual de 1 fs. Se us´o el mismo termostato y barostato que en la referencia [18], en concreto el termostato y barostato de Berendsen, con una temperatura de simulaci´on de 300K y una presi´on de 1.01325 bares en ambas simulaciones. Fundido: El estado inicial se prepar´o a partir de las cadenas parcialmente colapsadas. Se colocaron 70 cadenas polim´ericas de manera aleatoria en una caja de simulaci´on, realizando una simulaci´on breve de 50 bares de presi´on. Este proceso de preparaci´on es similar al proceso publicado en la referencia [30] y se usa para eliminar los huecos entre cadenas y 3.2. Sistema modelo 25 obtener una densidad homog´enea. A continuaci´on, se equilibraron los sistemas de forma similar como la publicada en la referencia [10]. La duraci´on de las simulaciones que se usaron para la recopilaci´on de los datos junto con otros detalles se encuentra resumida en el Cuadro 3.1. Detalles de los sistemas simulados Tipo N´umero Mw Nc T t de simulaci´on de mon´omeros [g/mol] [K] [ns] 10 1143.5 1 300 100 30 3426.3 1 300 100 Vac´ıo 50 5709.2 1 300/600 80/15 100 11416.4 1 300/600 123/7 125 14270.0 1 300/600 212/6 10 1143.5 1 300 1/100 30 3426.3 1 300 1/100 Disoluci´on 50 5709.2 1 300 1/100 100 11416.4 1 300 1/100 125 14270.0 1 300 1/100 10 1143.5 70 500 500 Fundido 30 3426.3 70 500 500 50 5709.2 70 500 500 100 11416.4 70 500 500 Cuadro 3.1: Detalles de las simulaciones realizadas. Donde Mw hace referencia al peso molecular de la cadena (Molecular weight, en ingl´es), Nc es el n´umero de cadenas en la simulaci´on, Tla temperatura de la simulaci´on en Kelvin y tel tiempo de la simulaci´on en nanosegundos. 26 Cap´ıtulo 3. Materiales y M´etodos 4. Resultados y Discusi´on 4.1. Desarrollo del campo de fuerzas El primer paso del dise˜no de un enfoque ascendente consisti´o en el desarrollo del campo de fuerzas describiendo las fuerzas interat´omicas existentes. Para ello, se ha tomado como referencia trabajos publicados anteriormente [25, 31]. El tr´ımero creado con el software AVOGADRO [27] ha servido para identificar diferentes tipos de ´atomos, tal y como se muestra en la Figura 4.1. Adem´as, se tom´o como referencia el trabajo [17]. En dicho trabajo, la mol´ecula de PCL termina con grupos vinilo, a diferencia de la mol´ecula generada en el presente trabajo cuyos grupos terminales son un ´acido carbox´ılico e hidroxilo, por lo que la simulaci´on se llev´o a cabo en la forma que aparece tambi´en durante la s´ıntesis. La notaci´on de los ´atomos se resume y explica en el Anexo en el Cuadro S1. 270 268 269 267 136 140 157 154 155 4671 4662 4651 156 490 469 176 Figura 4.1: Estructura inicial del tr´ımero usada para el desarrollo del campo de fuerzas y creada por AVOGADRO. Los n´umeros etiquetan los tipos de ´atomos seg´un la notaci´on usada de forma convencional por el OPLS. [15]. 28 Cap´ıtulo 4. Resultados y Discusi´on El conjunto final de par´ametros aplicados para el modelado del potencial enlazante y no enlazante se enumera en los Cuadros S2 y S3 en el Anexo. 4.2. Propiedades estructurales 4.2.1. Inspecciones visuales Se comenz´o el an´alisis de datos mediante la inspecci´on visual de los sistemas simulados. Para crear im´agenes de las simulaciones usamos el programa Visual Molecular Dynamics (VMD) [32]. (a) (b) (c) (d) Figura 4.2: (a): Captura de la configuraci´on inicial totalmente extendida de la cadena de 30 mon´omeros de PCL. (b): Captura de la configuraci´on semi colapsada de la cadena polim´erica. (c): Captura de la configuraci´on totalmente colapsada de la cadena de 30 mon´omeros tras la ejecuci´on de la simulaci´on en el vac´ıo. (d): Captura de la representaci´on de la caja de simulaci´on con la cadena simulada de PCL en su interior. Los colores asignados a los distintos ´atomos son ci´an para el carbono, gris para el hidr´ogeno y rojo para el ox´ıgeno. Inicialmente, se realiz´o una simulaci´on usando vac´ıo como medio. Se parti´o de una cadena de PCL totalmente extendida, empezando a colapsarse a medida que avanza la simulaci´on hasta llegar a una configuraci´on totalmente colapsada (“coil”). El proceso se puede apreciar siguiendo las capturas en la Figura 4.2 en el orden (a),(b) y (c). Esto se debe a que el vac´ıo se considera un “disolvente muy malo”. La cadena polim´erica se contrae porque busca interactuar consigo 4.2. Propiedades estructurales 29 misma reduciendo su superficie en vez de interactuar con el medio que, debido a que se encuentra en el vac´ıo, no existen mol´eculas con las que interactuar. La cadena individual en una caja de simulaci´on en el vac´ıo est´a mostrada en la Figura 4.2(d). El cambio de dimensiones de la cadena durante la simulaci´on en vac´ıo es apreciable tambi´en por las dimensiones de la caja, inicialmente es ligeramente m´as grande que la cadena totalmente extendida pero, al colapsar la mol´ecula, el tama˜no de la caja en comparaci´on con la cadena se ve muy superior. (a) (b) Figura 4.3: (a): Captura de la configuraci´on colapsada de la cadena de 30 mon´omeros de PCL tras la simulaci´on usando agua como disolvente. (b): Captura de la caja de simulaci´on albergando la cadena polim´erica y las mol´eculas de disolvente, en este caso, agua. Para una mejor visualizaci´on, las mol´eculas de agua son transparentes. Tras la simulaci´on en vac´ıo, se realizaron simulaciones utilizando agua como disolvente. Al igual que ocurr´ıa en el caso anterior, el pol´ımero adquiere una configuraci´on totalmente colapsada reduciendo su superficie y las interacciones con las mol´eculas del medio, por lo que podemos concluir que el agua es un mal disolvente para la PCL. Esta conclusi´on se puede confirmar visualmente comparando las configuraciones de la cadena de 30 mon´omeros en agua (Figura 4.3(a)) y en el vac´ıo (Figura 4.2(c)). Adem´as, se aprecian similitudes entre la Figura 4.3(b) y la Figura 2 en el trabajo publicado por Di Pasquale et al. [18], pudi´endose ver en ambas la cadena de PCL totalmente colapsada rodeada por las mol´eculas del disolvente (agua en ambos casos). Por ´ultimo, se analiza la simulaci´on utilizando un conjunto de 70 cadenas de PCL, imitando las condiciones en el fundido. En dicha simulaci´on se puede observar como la cadena ya no est´a colapsada, adquiriendo una configuraci´on totalmente distinta, por lo que podemos asumir que, a simple vista, el fundido no act´ua como un “mal disolvente” ya que la cadena busca en cierto modo interactuar con las mol´eculas del medio (Figura 4.4). 30 Cap´ıtulo 4. Resultados y Discusi´on (a) (b) (c) Figura 4.4: (a): Captura de una cadena polim´erica de 30 mon´omeros seleccionada aleatoriamente entre las 70 presentes en una simulaci´on en el fundido. (b): Captura de la cadena seleccionada dentro de la caja de simulaci´on. (c): Representaci´on de la caja de simulaci´on conteniendo la cadena seleccionada y las dem´as, en color lila, representando el fundido como medio denso. Para una mejor visualizaci´on, las cadenas alrededor de la cadena seleccionada son transparentes. 4.2.2. Cuantificaci´on de las inspecciones visuales A la hora de cuantificar las observaciones visuales, se procedi´o a calcular y definir 2 propiedades estructurales fundamentales en pol´ımeros: Radio de giro, Rg: El radio de giro se define como la distancia al cuadrado promedio entre los mon´omeros en una conformaci´on dada y el centro de masa del pol´ımero [2]. El valor medio est´andar se calcula habitualmente del valor promedio del radio de giro al cuadrado, tal y como aparece en la ecuaci´on 4.1: Rg=qR2 g,(4.1) 4.2. Propiedades estructurales 31 donde: R2 g=P i ||ri−rcm||2mi P i mi ,(4.2) mies la masa de la part´ıcula i,ries su posici´on y rcm es la posici´on del centro de masas del objeto medido. Los corchetes en la ecuaci´on 4.1 representan el promedio sobre el sistema previamente equilibrado. El error se calcula seg´un la teor´ıa de propagaci´on de errores [33] con la ecuaci´on correspondiente para la ra´ız cuadrada 4.3: ¯ Rg=1 2⟨R2 g⟩−0,5¯ R2 g,(4.3) donde ¯ R2 ges la desviaci´on est´andar obtenida. Distancia entre extremos, Re: La distancia entre extremos se define como la medida obtenida de la distancia cuadrado promedio entre ambos extremos de la cadena polim´erica [2]. Su valor al cuadrado se puede calcular gracias a la ecuaci´on 4.4: ⟨R2 e⟩=⟨||re−r0||2⟩,(4.4) donde 〈·〉 representa la media estad´ıstica de todas las cadenas polim´ericas y del tiempo, re representa la posici´on del ´ultimo ´atomo mientras que r0la posici´on del primer ´atomo. En analog´ıa con el radio de giro, el valor promedio y el error se calculan generalmente tal y como se ha mostrado en las ecuaciones 4.1 y 4.3. Relaci´on entre RgyRe: Existe una relaci´on matem´atica entre el radio de giro y la distancia entre extremos [2], tal y como muestra la ecuaci´on 4.5: R2 g=R2 e 6.(4.5) En la ecuaci´on 4.5 se puede comprobar como la relaci´on entre ambas propiedades es directamente proporcional para un valor promedio de un sistema previamente equilibrado (representado por corchetes) [2]. Dicha ecuaci´on es v´alida para una cadena lineal e ideal, entendi´endose por ideal que se trata de una cadena Gaussiana e infinitamente larga [2]. 32 Cap´ıtulo 4. Resultados y Discusi´on Figura 4.5: Representaciones gr´aficas de la evoluci´on del radio de giro, Rgy la distancia entre extremos, Re, con respecto al tiempo en los 3 tipos de simulaci´on. 4.2.3. Evoluci´on de las propiedades estructurales en el tiempo En Figura 4.5 se representan las propiedades estructurales definidas en el apartado anterior con respecto al tiempo para los 3 tipos de simulaciones realizadas con una cadena de 50 mon´omeros. Se observaron tendencias similares en las simulaciones de las cadenas de longitudes restantes. En primer lugar, en las Figuras 4.5(a,b) se puede observar c´omo con la simulaci´on en vac´ıo conseguimos que la mol´ecula colapse r´apidamente adquiriendo esa conformaci´on de “coil” descrita anteriormente. La simulaci´on a 300K (Figura 4.5(a)) presenta un descenso progresivo a lo largo del tiempo, no pudi´endose establecer con claridad el momento en el que la cadena 4.2. Propiedades estructurales 33 colapsa, de hecho, la simulaci´on se detuvo antes de que eso ocurriera. Es por ello que se hizo la simulaci´on a 600K, ya que las mol´eculas adquieren una mayor energ´ıa cin´etica pudi´endose obtener un estado colapsado de manera pr´acticamente instant´anea. Se puede comprobar que subiendo la temperatura en la simulaci´on a 600K (Figura 4.5 (b)) ocurre un dr´astico descenso pr´acticamente al inicio de la simulaci´on, indicando el colapso de la mol´ecula en la caja de simulaci´on. Despu´es de la ca´ıda, el valor de Rgfluct´ua sobre un valor promedio en el tiempo, que indica que se ha alcanzado el estado colapsado de la cadena. El motivo de que exista esta tendencia descendente en ambos casos, se debe a que el sistema en vac´ıo no se encuentra equilibrado. Una vez obtenida la mol´ecula colapsada, se realiz´o la simulaci´on en disoluci´on usando agua como disolvente. Se comprueba como ahora los valores de RgyRes´ı fluct´uan sobre un valor promedio ya que el sistema ha sido equilibrado previamente. En el caso del fundido, el sistema ha sido equilibrado previamente al igual que en caso de las gr´aficas (c) y (d). Sin embargo, se ve c´omo en el fundido la fluctuaci´on tanto en Rgcomo en Rees sustancialmente m´as oscilante que en agua. En las gr´aficas (e) y (f) se pueden apreciar visualmente c´omo las fluctuaciones para RgyRe son m´as oscilantes que en los casos (c) y (d), por lo que se podr´ıa asumir que la cadena polim´erica fluct´ua m´as en el fundido que en disoluci´on. Sin embargo, es posible calcular la oscilaci´on de manera porcentual con la ecuaci´on 4.6: Fluctuaci´on = ∆x x×100,(4.6) donde ∆xrepresenta el error y xel valor promedio de la propiedad en cuesti´on. De esta forma, se puede comprobar que la fluctuaci´on para la simulaci´on en disoluci´on para la cadena de 50 mon´omeros es de 2.7 %, mientras que para la disoluci´on en el fundido la fluctuaci´on es del 2.4 %. Es decir, que son pr´acticamente iguales a pesar de que son visualmente distintas. El c´alculo anterior se realiz´o para todas las longitudes de cadena, pudi´endose as´ı verificar que las fluctuaciones de la cadena en ambos medios difieren de manera insignificante. 4.2.4. Evoluci´on de las propiedades estructurales frente a N Para demostrar la influencia que posee el n´umero de mon´omeros de la cadena, N, en el radio de giro y la distancia entre extremos en diferentes medios, se han recopilado los datos promedios para todos los sistemas. Los resultados obtenidos se muestran en el Cuadro 4.1. En las representaciones gr´aficas obtenidas del radio de giro y distancia entre extremos frente al n´umero de mon´omeros por cadena, se ha decidido prescindir de la simulaci´on en vac´ıo. Esto se debe a que los resultados obtenidos distan mucho de lo esperado y carec´ıan de sentido, el 40 Cap´ıtulo 4. Resultados y Discusi´on Propiedades de sistemas en el fundido Tipo N´umero Densidad τ β de simulaci´on de mon´omeros [kg/m3] [ps] 10 928 ±7 1440 0.77 Fundido 30 928 ±4 3223 0.41 50 928 ±3 12100 0.39 100 918 ±3 318000 0.46 Cuadro 4.2: Resumen de los valores de densidad y los tiempos de relajaci´on para el sistema en fundido. El par´ametro βes el coeficiente de estiramiento. Por ´ultimo, es importante se˜nalar que en vistas a la construcci´on de un modelo computacional ascendente, el c´alculo de propiedades clave tales como los tiempos de relajaci´on, la densidad, RgyRepara diferentes longitudes de cadena es un paso esencial para el procedimiento de paso de par´ametros. Tal y como se mostr´o en el trabajo de Behbahani et al. [9], para obtener el conjunto completo de par´ametros, ser´ıa necesario ampliar el an´alisis para definir y hacer coincidir el empaquetamiento local y tambi´en el movimiento difusivo entre los modelos que operan en diferentes escalas de longitud y tiempo. 4.3. Propiedades din´amicas 41 Figura 4.7: (a): Representaci´on de la funci´on de autocorrelaci´on orientacional del vector entre extremos para la cadena de 30 mon´omeros. (b): Representaci´on de la variaci´on del tiempo de relajaci´on en funci´on del n´umero de mon´omeros de la cadena. 42 Cap´ıtulo 4. Resultados y Discusi´on 5. Conclusiones En el presente trabajo se ha presentado un procedimiento sistem´atico con enfoque ascendente para un pol´ımero biodegradable, concretamente policaprolactona. El m´etodo computacional elegido para estudiar estos sistemas ha sido las simulaciones atom´ısticas con una simulaci´on de din´amica molecular en 3 diferentes medios: el vac´ıo, en disoluci´on y en el fundido. Los sistemas tambi´en han sido simulados con diferentes longitudes de cadena polim´erica. En primer lugar, se desarroll´o su campo de fuerzas bas´andose en tres modelos: OPLS-AA, LOPLS-AA y se us´o informaci´on de un modelo del ´acido polil´actico para un valor de un ´angulo debido a que no exist´ıa bibliograf´ıa que contuviera ese dato para la PCL. Todos los detalles se encuentran adjuntos en el Anexo. Para verificar que el campo de fuerzas creado concuerda con algunos resultados ya publicados, se procedi´o a calcular diferentes propiedades estructurales y din´amicas del pol´ımero. Se han construido los sistemas con las distintas cadenas polim´ericas con las longitudes deseadas y se han equilibrado para poder realizar las simulaciones propuestas. Se han monitorizado las propiedades de estudio para verificar que los sistema han sido equilibrados correctamente, y se han recopilado todos los datos para poder realizar el c´alculo del promedio de las propiedades a estudiar. Se han estudiado dos propiedades estructurales fundamentales en pol´ımeros: el radio de giro (Rg) y la distancia entre extremos (Re). Los datos han sido comparados con trabajos publicados previamente, en concreto se ha tomado como referencia principal los trabajos de Di Pasquale et al. [18, 19]. En base a los resultados obtenidos, se puede concluir que ambas propiedades aumentan su valor a medida que aumenta el n´umero de mon´omeros de la cadena polim´erica, de acuerdo con la ley de potencia cuyo exponente corresponde a 0.5 para pol´ımeros ideales y 1/3 para un disolvente pobre. Estos exponentes, de acuerdo con la precisi´on que se puede obtener en la simulaci´on, concuerdan con lo que se podr´ıa esperar para los resultados obtenidos en disoluci´on usando agua como disolvente, que es un disolvente malo, y en el fundido, que deber´ıa aproximarse m´as a un sistema ideal. Finalmente, se calcularon dos propiedades de los sistemas en el fundido: la densidad y el 44 Cap´ıtulo 5. Conclusiones tiempo de relajaci´on. En los sistemas simulados, se puede observar que los valores de densidad resultan independientes de la longitud de cadena, sin embargo, el tiempo de relajaci´on posee una relaci´on directamente proporcional a la longitud de la cadena, de acuerdo con el comportamiento te´orico esperado. Los datos recopilados de propiedades estructurales y din´amicas resultan de gran utilidad para construir un modelo mesosc´opico, que ser´ıa el siguiente paso del enfoque sistem´atico ascendente. 6. Perspectivas de futuro La investigaci´on desarrollada en este trabajo establece una base s´olida para futuros estudios en el campo de los pol´ımeros sostenibles, con un enfoque espec´ıfico en la policaprolactona (PCL). En las pr´oximas d´ecadas, se espera que los avances en la simulaci´on computacional y las metodolog´ıas de din´amica molecular contribuyan significativamente a mejorar las propiedades de la PCL, haci´endola m´as competitiva frente a los pol´ımeros sint´eticos en t´erminos de resistencia mec´anica y estabilidad t´ermica. Uno de los principales retos a abordar ser´a la optimizaci´on de la PCL para aplicaciones industriales, como la fabricaci´on aditiva y la biomedicina, donde ya muestra un gran potencial debido a su biodegradabilidad y biocompatibilidad. La personalizaci´on de campos de fuerzas y modelos atom´ısticos permitir´a predecir con mayor precisi´on su comportamiento en entornos complejos, lo que facilitar´a la integraci´on de este biopol´ımero en aplicaciones de alto rendimiento. En el presente trabajo, se ha mostrado la eficiencia de un enfoque computacional ascendente capaz de brindar resultados cuantitativos prediciendo diversas propiedades estructurales y din´amicas de los pol´ımeros. Adem´as, se anticipa que el desarrollo de nuevas t´ecnicas de simulaci´on multiescala, que combinen enfoques atom´ısticos y mesosc´opicos capaces de predecir la reolog´ıa del pol´ımero, permitir´a abordar la relaci´on entre estructura molecular y propiedades macrosc´opicas de manera m´as eficiente en diversas aplicaciones, por ejemplo, la impresi´on 3D. 46 Cap´ıtulo 6. Perspectivas de futuro Bibliograf´ıa [1] D. I. Bower, An introduction to polymer physics. Cambridge University Press, 2002. [2] M. Rubinstein and R. H. Colby, Polymer physics. Oxford University Press, 2003. [3] 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. [4] 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. [5] M. Thakur, I. Majid, S. Hussain, and V. Nanda, “Poly (ε-caprolactone): A potential polymer for biodegradable food packaging applications,” Packaging Technology and Science, vol. 34, no. 8, pp. 449–461, 2021. [6] M. P. Allen, Introduction to Molecular Dynamics Simulation, vol. 23 of NIC series, pp. 1–28. J¨ulich: John von Neumann Institute for Computing, 2004. [7] J. Jin, A. J. Pak, A. E. Durumeric, T. D. Loose, and G. A. Voth, “Bottom-up coarsegraining: Principles and perspectives,” Journal of chemical theory and computation, vol. 18, no. 10, pp. 5759–5791, 2022. [8] F. R. Hung, K. E. Gubbins, and S. Franzen, “A graduate course on multi-scale modeling of soft matter,” Chemical Engineering Education, vol. 38, no. 4, pp. 242–249, 2004. [9] A. F. Behbahani, L. Schneider, A. Rissanou, A. Chazirakis, P. Bacova, P. K. Jana, W. Li, M. Doxastakis, P. Polinska, C. Burkhart, et al., “Dynamics and rheology of polymer melts via hierarchical atomistic, coarse-grained, and slip-spring simulations,” Macromolecules, vol. 54, no. 6, pp. 2740–2762, 2021. [10] E. Christofi, P. Baˇcov´a, 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. 48 Bibliograf´ıa [11] Y. Nie, Multiscale investigation on the physical properties of polycaprolactone and its composites by molecular dynamics and coarse grain simulation. PhD thesis, Queensland University of Technology, 2023. [12] C. S. Ezquerro, J. M. G. Aznar, and M. Laspalas, “Prediction of the structure and mechanical properties of polycaprolactone–silica nanocomposites and the interphase region by molecular dynamics simulations: the effect of pegylation,” Soft Matter, vol. 18, no. 14, pp. 2800–2813, 2022. [13] Z. U. Arif, M. Y. Khalid, R. Noroozi, A. Sadeghianmaryan, M. Jalalvand, and M. Hossain, “Recent advances in 3d-printed polylactide and polycaprolactone-based biomaterials for tissue engineering applications,” International Journal of Biological Macromolecules, vol. 218, pp. 930–968, 2022. [14] M. Drenscko and S. M. Loverde, “Molecular dynamics simulations of the interaction of phospholipid bilayers with polycaprolactone,” Molecular Simulation, vol. 45, no. 11, pp. 859– 867, 2019. [15] W. L. Jorgensen, D. S. 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. [16] H. Sun, S. J. Mumby, J. R. Maple, and A. T. Hagler, “An ab initio cff93 all-atom force field for polycarbonates,” Journal of the American Chemical society, vol. 116, no. 7, pp. 2978– 2987, 1994. [17] I. Yungerman, I. Starodumov, A. Fulati, K. Uto, M. Ebara, and Y. Moskovitz, “Fullatomistic optimized potentials for liquid simulations and polymer consistent force field models for biocompatible shape-memory poly (e-caprolactone),” The Journal of Physical Chemistry B, vol. 126, no. 21, pp. 3961–3972, 2022. [18] N. Di Pasquale, D. Marchisio, P. Carbone, and A. Barresi, “Identification of nucleation rate parameters with md and validation of the cfd model for polymer particle precipitation,” Chemical Engineering Research and Design, vol. 91, no. 11, pp. 2275–2290, 2013. [19] N. Di Pasquale, D. L. Marchisio, A. A. Barresi, and P. Carbone, “Solvent structuring and its effect on the polymer structure and processability: The case of water–acetone poly-- caprolactone mixtures,” The Journal of Physical Chemistry B, vol. 118, no. 46, pp. 13258– 13267, 2014. Bibliograf´ıa 49 [20] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, vol. 1 of Computational Science Series. Academic Press, second ed., 2002. [21] S. Nos´e, “A molecular dynamics method for simulations in the canonical ensemble,” Molecular physics, vol. 52, no. 2, pp. 255–268, 1984. [22] V. A. Harmandaris and V. G. Mavrantzas, “Molecular dynamics simulations of polymers,” in Simulation methods for polymers, pp. 178–218, CRC Press, 2004. [23] 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. [24] 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-2, pp. 19–25, 2015. [25] S. W. Siu, K. Pluhackova, and R. A. Bockmann, “Optimization of the OPLS-AA force field for long hydrocarbons,” Journal of Chemical theory and Computation, vol. 8, no. 4, pp. 1459–1470, 2012. [26] B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, “LINCS: A linear constraint solver for molecular simulations,” Journal of Computational Chemistry, vol. 18, no. 12, pp. 1463–1472, 1997. [27] Avogadro Chemistry, “Avogadro.” http://avogadro.cc/. [28] P. Baˇcov´a, E. Glynos, S. H. Anastasiadis, and V. Harmandaris, “Nanostructuring singlemolecule polymeric nanoparticles via macromolecular architecture,” ACS Nano, vol. 13, no. 2, pp. 2439–2449, 2019. [29] H. J. Berendsen, J. P. Postma, W. F. van Gunsteren, and J. Hermans, “Interaction models for water in relation to protein hydration,” in Intermolecular forces: proceedings of the fourteenth Jerusalem symposium on quantum chemistry and biochemistry held in Jerusalem, Israel, April 13–16, 1981, pp. 331–342, Springer, 1981. [30] A. D. Glova, S. G. Falkovich, S. V. Larin, D. A. Mezhenskaia, N. V. Lukasheva, V. M. Nazarychev, D. A. Tolmachev, A. A. Mercurieva, J. M. Kenny, and S. V. Lyulin, “Poly(lactic acid)-based Nanocomposites Filled with Cellulose Nanocrystals with Modified Surface: Allatom Molecular Dynamics Simulations,” Polym. Int., vol. 65, no. 8, pp. 892–898, 2016.