UNIVERSIDAD COMPLUTENSE DE MADRID Facultad de Ciencias Químicas Departamento de Química Física ¡ < UNIVERSIDAD COMPLUTENSE 5310313745 SIMULACION DE LIQUIDOS MOLECULARES FLEXIBLES ~ L ¡ —~ 2’ Noé García Almarza Madrid, 1993 Colecc¡ón Tesis Doctorales. N.0 142/93 © Noé García Almarza Edita e imprime la Editorial de la Universidad Complutense de Madrid. Servicio de Reprograf la. Escuela de Estomatología. Ciudad Universitaria. Madrid, 1993. Ricoh 3700 Depósito Legal: M-12337-1993 La Tesis Doctoral de a. NOR GARCíA ALMARZA SIML’LACTON LE LIQUILCE MOLECULARESTitulada FLEXIBLES! Director Dr.. O. . ..~P.A~I~QI$.ILF~ . fue leída tutía Facultad.de ....cJ1~QIn.caz::As.. de la UNIVERSIDAD COIU>LLITENSE DE MADRID. el dfa .h~t.. da S~PTTEMBRE dc 19 .9k. ante el tribunal constituido por las siguientes Profesores: nrswsrr U. JUAN JOSE FREIRE COMEZ VOCAL U. MANUEL LOMBARL ~ yuca U - JUAN COLMENERO DE LEON VOCAL U. MIGUEL RUEI CAPACETTI src~nAaxa - U~ FAZ PADII~A PI. £{‘-L~&~~QR habiendo recibido la calificación de . • -. . - j. - . Madrid, a lo de (tC £&ú$Éde.19 tj.( EL SECRETARIO bEL TRIBUNAL. UNIVERSIDAD COMPLUTENSE DE MADRID FACLLTAD DE CIENCIAS QUíMICAS SIF&LACICN DE LíQUIDOS MOLECILAHES FLEXIBLES MEMORIA QUE PARA OPTAR AL GRADO DE DOCTOR EN CIENCIAS QUIMICAS PRESENTA NOE SARCIA ALMARZA DIRECTOR D~. EDUARDO ENCISO RODRíGUEZ DEPARTAMENTO DE QUíMICA FíSICA UNIVERSIDAD CONfLUTENSE MADRID. JUNIO DE 1991 Este trabajo ha sido realizado en el Departamento de Química Física 1 de 1. Facultad de Ciencias Químicas de la Universidad Complutense de Madrid, bajo 1. dirección del Dr. O. Eduardo Enciso Rodríguez, a quien deseo expresar mi más sincero agradecimiento por su constante apoyo. Mucha ha sido la gente que, con sus sugerencias y experiencia en diversos campos de la ciencia, ha contribuido al desarrollo de este trabajo; quiero, por ello, citar aquí las interesantes discusiones mantenidas con las siguientes personas: Dr. David Brown, Dr. Marvin Bishop. Dr. Jean—Paul Ryckaert. Dr. Carlos Vega. Dr. Julian MR. Clarke, Dra. Paz Padilla, Dr .Idrn lija Siepmann, Dr. Jesús Alonso. Dr. Javier Bermejo y nra. )~arta Alvarez. Deseo agradecer las facilidades recibidas para la realización de este trabajo por parte de diversos centros de cálculo: Centro de Proceso de Datos de la Universidad Complutense, Sisten Informático dc Ciencias de la Universidad Complutense. Centro de Cálculo del Consejo Superior de Investigaciones Científicas y University of Manchester Computing Centre. Para un buen aprovechamiento de los recursos informáticos de los citados centros tuve la fortuna de contar con la valiosa ayuda de las siguientes personas: Dr. Enrique Lomba, Dr. Claudio Hartin, Dr. J 05é Luis Fernández Abascal, Dr. David Brown y Dr. Sieve Lies; a los cuales deseo expresar mi agradecimiento. También quiero expresar mi agradecimiento al Departamento de Ovimica Física Y por las facIlIdades y el apoyo recibidos para la realización de este trabajo. ~gradezcoa los compafieros del Departamento, a si familia y a mis amigos su constante apoyo, especialmente durante el proceso de escritura de esta memoria. Finalmente, agradezco al Ministerio de Educación y Ciencia la concesión de una beca eorrespondiente al Programa de Formación y Perfeccionamiento de investigadores en la utilización de la radiación de sincrotron y de los haces de neutrones, así como a la Dirección General de Investigación Científica y Técnica por sus ayudas ?BSE—0617—C02—02 y PBSS-OO27—COS—03. INDICE 1. INTRODIJCCION 5 1.1. MECANICA ESTADISTICA CLÁSICA II 1.2. C~RDENADAS lNT~t~AS 21 2. MODELOS DE POTD4CIAL 27 2. 1 MODELO SIMPLIFICAW PARA N—ALCMIOS Si 2.2 MODELO SIHPLIFICAm PARA l.2-DICi.OROETANO 35 2. 3 MODELO REALISTA DE N-BJTANO 39 3. TECNICAS DE SImUCION 45 3.1 ME1tDOS DE MONTE CARLO Sí 3. 1. 1 ALGORITMOS DE TIPO HETROPOLIS 59 3.2 DINAJ4ICA NOLEC~JLAR 73 3.2. 1 SISTEJ4AS CON LIGADURAS 77 3.3 ANALISIS DE CONFIGURACIONES. TECHICAS UMBRELLA 99 3.4 ESTIRACION DE ~RORES 95 4. SIMUI.ACION POR MONTE-CARLO DE NOLECULAS FLEXIBLES 101 4. 1 SIIULACION EFICIENTE DE N—BIJTANO 101 4. 2 METO~S DE R~TACION 133 4.2.1 ANAIISIS DE EFICIENCIA EN REPTACION 139 5. RESULTAmS 149 5. 1 N-BUTANO LIQUI~ (MODELO SIMPLIPICAm) 149 5. 1. 1 EFECTO DE LAS IN1nÁCCI0N~ INTERMOLEC1JLARES 157 5.2 SI1ffJL.ACION DE 1,2-DICIOROETANO 161 5.3 SIWJLACIDN DE N-AICANOS 169 5.3. 1 SMJLACION POR DINM4ICA MOLECULAR 193 5.3.2 REFINM4IENTO DE PARÁMETROS DE POTENCIAL 209 5.4 H—BUTANO LIQUIDO (MODELO DETALLAm) 229 6. CONCLUSIONES 243 2 7. ANEXOS 245 A DE3IVACION DE FU~2AS INTRMCLEOJLARES EN SISTfl4AS CON LIGADURAS RíGIDAS 245 B. CALCULO DE PRESION 251 C. TRAYSFOPflACION DE LA INTECRAL DE CONFICLJRACION 257 0. RANSFOBMACIONES DE COOHDENADAS 269 E. DXNM4ICA DE SISTfl4AS CON LIGADURAS 279 E. MUESTREO DE DR IENTACIONES MOLECULARES 287 II. CALCULO DEL TENSOR NETRICO PARA CADENAS CON LIGAflURAS RíGIDAS 291 E. REFEPENCIAS 299 3 1. INTRODUCCION El estudio de sistemas constituidos por moléculas flexibles constituye un campo de investigaci¿n de creciente desarrollo en los iitiaos ellos. La característica más diferenciadora de estas eoiécuies respecto a otros tipos de compuestos es su posibilidad de sufrir cambios conformacionales estos cambios pueden lugar a coaportaaientos muy variados desde el punto de vista 2 físico—químico La estructura molecular y las propiedades dinámicas de las proteínas son determinantes para su función, tales propiedades dependen de modo fundamental de la conformación molecular, que puede cambiar notablemente con las características del medio. La conformación molecular juega, por tanto, un papel muy importante en procesos bioquímicos (reacciones enzimátlcas. desnaturalización de proteínas, etc). La flexibilidad zolacular detern,ina, así sismo, las propiedades de moléculas orgánicas flexibles, así coso de los sistemas poliméricos. La capacidad que tienen las macromoléculas do adoptar distintas conformaciones da lugar a la diversidad de propiedades físicas y aplicaciones tecnológicas de estos materiales. En el estudio de este tipo de sistemas me emplean diversos modelos para represeniar las moléculas’t Para tratar aspectos fundamentalmente teóricos de la fisica de polimeros se utilizan ,, 5 frecuentemente modelos simplificados . que hacen abstracción del carácter químico de los monómeros que constituyen una determinada macromolécula. Estos modelos permiten estudiar propiedades termodinámicas, estructurales y dinámicas de los sistemas en función del tasaóo del polímero (leyes de escala), de la naturaleza de las interacciones intermoleculares, de las interacciones con un disolvente. En el estudio de moléculas flexibles de Sajo peso molecular es posible utilizar modelos más realistas; en este caso los átomos o los grupos de átomos que definen desde un punto de vista químico la molécula son tenidos en cuenta explícitamente al definir los potenciales de interacción. facilitándose de este modo la comparación de los resultados obtenidos a partir del Modelo con resultados experimentales sensibles a la identidad molecular. Las técnicas de Simulación por ordenador constituyen hoy en día uno de los instrumentos más fructíferos en diversidad de campos de la ciencia. La física de las fases desordenadas de alta densidad, y en particular la física del estado liquido, es uno de los campos de la investigación científica donde la introducción de las técnicas de simulación ha sido fundamental para los desarrollos producidos en los últimos años. Este hecho se debe a las dificultades de tipo matemático que surgen en los tratamientos teóricos de estos sistemas, y que conducen inevitablemente a la necesidad de formular teorías aproximadas sobre modelos que 6 1. IrsoouccloU tienen en cuenta los elementos físicos básicos de la fenomenología que se desea estudiar. Ante esta situación la simulación por ordenador puede aportar resultados muy valiosos, tanto para contrastar la validez de las aproximaciones teóricas, como para obtener resultados con caracter predictivo de propiedades físicas experimentales. Esta segunda aplicación de las técnicas de simulación es cada día mas importante por la constante sejora de prestaciones y precios de los ordenadores. Las técnicas de simulación se utilizan también en la interpretación de resultados experimentales de tipo espectroscópico 8. En el caso de moléculas flexibles los experimentos de dispersión de neutrones permiten obtener información sobre la distribución conformacional9’ ‘~. Sin embargo. a menudo no es fácil extraer conclusiones definitivas a partir de tales experimentos, debido a que las distancias lntrsmoleculares que cambian apreciablemente con la conformación molecular son del sismo orden que algunas de las distancias intermoleculares. En eetos casos, la simulación por ordenador constituye un instrumento indispensable para separar las distintas contribuciones al factor 9 de estructura En este trabajo se ha abordado el estudio de fluidos constituidos por moléculas flexibles mediante métodos de simulación por ordenador. Por una parto se estudió la influencia del medio y del tamaño de la cadena en la estructura de moléculas de n—alcanos, por otra parte se investigó el efecto de las interacciones intermoleculares de tipo electrostático sobre el equilibrio conformacional utilizándose en este caso como ejemplo la molécula de 1.2 dicloroetano. Finalmente se abordó la comparación de los resultados de simulación, con resultados estructurales obtenidos mediante técnicas de dispersión de neutrones. Una de las particularidades del tipo de moléculas estudiadas en este trabajo es la existencia de diversas conformaciones estables separadas por barreras de alta energía. Este hecho tiene consecuencias importantes desde el punto de vista físico—químico (situaciones de no equilibrio, fases metaestables, fenómenos de histéresis, largos tiempos de relajación de la estructura intra— e intermolecular, etc). Por otra parte, ia existencia de las barreras energéticas da lugar a dificultades desde el punto de vista computacional, haciéndose necesario en multitud de casos el 11.42 desarrollo de nuevas técnicas de simulación para obtener con la precisión requerida las propiedades de los sistemas simulados. En este trabajo se han desarrollado algunos métodos para la simulación eficiente de sistemas constituidos por moléculas flexibles; estos métodos han permItido establecer interesantes conclusiones sobre los distintos factores que influyen en el equilibrio conformacional de las moléculas flexibles. En las secciones 1.1 y 1.2 de este capitulo introductorio se presentan los conceptos de la Mecánica Estadística utilizados en este trabajo, así como una breve descripción de las coordenadas e U ItflSOOtJCCIoN internas empleadas en las simulaciones. El capitulo 2 está dedicado al problema de los potenciales de interacción, recogiéndose los modelos utilizados en este trabajo. En el capitulo 3 se presenta la descripción general de las técnicas de Simulación más frecuentemente empleadas -en el estudio de fluidos moleculares, prestándose especial atención a las aplicaciones a modelos que tienen en cuenta grados de libertad intraaoleculares. En el capitulo 4 se describen en detalle los procedimientos desarrollados en este trabajo para la simulación de moléculas flexibles. En el capitulo 5 se recogen los resultados obtenidos, análizandose las conclusiones que se pueden extraer de ellos. Este análisis hace hincapié en la comparación con resultados experimentales, en los efectos físicos observados para los distintos sistemas y en la eficiencia de las distintas técnicas de 5iaulac 1 ón. En el capitulo 6 se resumen las conclusiones generales. Finalmente se incluyen vatios anexos, donde se tratan diversas cuestiones técnicas del trabajo realizado. io 1. iWTROOIXCIO ti MECANICA ESTADíSTICA O..ASICA En esta sección se recogen exclusivamente algunas de las ecuaciones y definiciones que se usarán más adelante. 33 La función de partición en el colectivo microcanónico Q,w~ (número de partículas, N. volumen. V, y energía total, t. constantes) se define como: 314 a~• 1 dq’ Mdp~” fi[fl(X .1> —6 1 (1.1.1) Nt 3M 3M q representa las coordenadas de las partículas. p los momentos conjugados, 3< es la función energía total, función de coordenadas y momentos. h es la constante de Planck y A es la función delt, de Dirac. La función de partición en el colectivo canónlcoS ~ (número de partículas, volumen y temperatura absoluta. 7. constantes) es: t s3W Qwvt — 1 dfldp” exp (—$fl(E .p >1 (1.1.21 NI > 13W > Siendo $—i’kT, donde k es la constante de Boltzaann. En sIstemas conservativosl~íS. en los que la energía total — Uno depende explícitamente del tiempo, la energía total fl(R .p 11 puede separarme, trabajando en coordenadas cartesianas, en energía potencial, ti. función de las posiciones y energía cinética, 3<, función de Los momentos: 311 35 31< 3Mfl(X .p, 1 ) a ¶i(X 1 • X(p~ 1 (1.1.3) 35 2 — (1.1.4) ‘—a Donde m~ representa la asta asociada al momento p¡. Integrando la ecuación (1.1.2) respecto de los momentos se obtiene para el caso en que todas las masas son iguales: a .......faexp i—pinx 3Mn (1. 1.5) Id! A Siendo A: A * (h2fl~aT>’12 (1.1.6> La función de partición configuracional se define, trabajando en coordenadas cartesianas como, ZSVT.J dx” np I—$ii(X’)1 (1.1.7) 12 1. IrTRoBlJccloN La probabilidad de una cierta configuración en el subespacio de las posiciones será: 314 314 3*’ 38p(X ) dI — np (—0¶I(X )1 dX Ci. 1.8) Cuando se utilizan coordenadas generalizadas la energía cinética puede depender de coordenadas y momentos, RCX<4 3*tpN>. La función de partición conf iguracional se puede expresar en este •tU 16—caso como Ci. 1.9) El termino aparece al integrar generalizados (véase anexo C). El valor medio propiedad, A. que depende exclusivam.nt, de las obtendría” como: 111 i’2 311 A(q ) 35 35 >/2 los momentos de una cierta coordenadas se (1. 1.10) La función de partición canónica se relacione con la energía libre de 4 mediante la ecuación: ACM, Y, Ii a kT la ~ 11.2.11> 13 A partir de las relaciones de la teraodinásica se pueden evaluar diferentes propiedades. Por ejemplo el valor medio de la presión en el colectivo canónico viene dado por: p — (ii. i2> Utilizando las ecuaciones (í.I.Ií) y (i.i.12) pueden determinarse diferentes maneras de evaluar la presión del sistema mediante simulación en el colectivo canónico ~y (1.1. 13) Teniendo en cuenta la relsción entre energía cinética y temperatura~ 3 3< - NkT (1.1. i4) Se obtiene: ti Nk • (BU/OT) ~2] (1.1. 19) La estructura de los fluidos suele analizarse mediante la función de distribución radial1’~~, g(r). La definición de esta función puede plantearse utilizando funciones de probabilidad para las posiciones de varias partículas. De este modo P (r>.r 2) define la probabilidad de que las partículas 1 y 2 estén situadas, respectivamente en las posiciones r, y r2 independientemente de las posiciones del reato de las partículas del sistema: P’ 2> (r>,r 2)dr,dr2 — d.,dr¡ fdrs.. fa.-. g.-’> (1.1.201 1$ Donde p(r”> representa la función densidad de probabilidad para las configuraciones de un sistema de <4 partículas. Análogamente para la probabilidad de que la partícula 1. esté situada en la posición r,, independientemente del resto de partículas se define la función Nr, Y; P(r,)dr, — dr, fdr 2...fdra p(r”) (1. 1.21> Para un sistema homogéneo: P(r, )dr, — dr, Jdra.. .fdrs pCi.N) (1. 1.22) P/V (1.1. 23) Integrando la ecuación (1. 1.20) sobre las posicione, de la partícula 2 se obtiene: fdt2 P(r,,r2) a P(r,) (1.1.24> La probabilidad de que una partícula cualquiera ocupe un •Iemento de volumen dr, alrededor de r,. n(r, > será: n(r,)dr, — Id P(r,)dr, (I.i.25) nf,, )dr, — (N/V>dr, — n dr, (1.1.26> 16 i. turuobúccto« Siendo n la densidad de partículas en el sistema. La probabilidad de que ¿os partículas cualesquiera ocupan las posiones r, y r 2. n>Z>(r, r2>, será: (2> (2>n (r,,t2)dr,dr3 a N(Id—i) P (r1,r2) (1. 1.27) a> (r,.r2)a (r~a-í> * alt,) Mr2> (aí a n tr,,r2) 1 n De acuerdo con la ecuación (1.1.22) se cumplirá: nf g’flr,.ra> dr2 — N—l La función de correiaclón de orden 2. d.fhilda entre átoaos para sistemas homogéneos e isótropo, dependerá únicameñte de la distancia entre partículas. ¡r:—r>j: (a>g (r>.r2) g(1r2—r,[) Donde g es la función de dlstribucidn radial. g(r) verifica: (1. 1.28) (1.i.29> (1. 1.301 17 (1.1.31) Estadisticamente g(r) puede plantearse coso: < ,.v ~ 85 (1.1.30) En este trabajo se han incluido otras funciones de distribución, por ejemplo para considerar la estructura intrm,iolecular’ de moléculas flexibles se utilizaron funciones 5(r) definidas sobre grupos pertenecientes a una misma molécula: 5(r) — < Z Z~~~’j I—r> > (1.1.31) Ea esta ecuación el doble sumatorio se realiza sobre todos los pares de grupos pertenecientes a una misma Solénla. Integrando se obtiene: IS(r) dr • 14114•—l) /2 (1.1.32) Donde N~ es el número de átoaos por molécula. Ami sismo se definieron funciones de distribución para ángulos de rotación interna (sección 1.2): 19 5. lNTROOJcCIOk Donde •> representa un cierto ángulo de torsión de una determinada ao)écula. De modo análogo se utilizaron otras funciones para los distintos grados de libertad intraaoleculares. 19 20 1. INIXCIO« 1.2 COORDENADAS INTERNAS Existen diversas posibilidades de describir la posición de los átomos que constituyen una cierta molécula. Una opción es la utilización de las coordenadas cartesianas atómicas. Desde el punto de vista práctico resulta frecuentemente más útil utilizar coordenadas generalizadas’ ~ 15~ Las fuerzas relacionadas con el enlace químico condicionan de manera importante el movimiento de los átomos, por ello la manera habitual de elegir un conjunto de coordenadas generalizadas consiste en separar coordenadas internas y coordenadas externas. Las primeras dan cuenta de las distancias entre átomos de una misma molécula (distancias intrimoleculares), las coordenadas externas indican la posición y orientación de la molécula et~ el espacio. Dadas las coordenadas atómicas cartesianas es posible definir unas ciertas relaciones que permitan construir un sistema de ejes cuya orientación respecto al sistema de referencia del laboratorio permite definir la orientación ‘4 molecular La elección habitual de coordenadas internas para cadenas linealesl 4>? con >4 centros consiste en Id—l distancias de enlace, <4—2 ángulos de enlace, QYZ y N—3 ángulos de torsión (o de rotación interna), •»3. Estas coordenadas internas están esquesatizada. en la figura 1.2.1. La formulación matemática de 21 las coordenadas utilizadas a largo de este trabajo se recoge en el anexo 0. Desde el punto físico las distintos tipos de coordenadas internas utilizados en la descripción de la estructura tr~traisolecwlar de cadenas lineales son de muy diferente naturaleza’. Las variaciones de distancias y ángulos de enlace respecto de sus valores sedios en el equilibrio son poco importantes, debido al efecto de las energías de los enlaces moleculares. Los ángulos de torsión, sin embargo pueden variar considerablemente, dando lugar a caeMos notables en la geometría molecular; por este motIvo gran cantidad de modeles utilizados en el estudio de estos sistemas consideran distancias de enlace (y a menudo también ángulos de enlace) fijos. Las coordenadas de estos sistemas podrán clasificarse’ • ‘< por tanto en rígidas (que variarán poco respecto de sus valores medios) y coordenadas flexibles. Desde un punto de vista formal en este planteamiento se considera que no existen acoplamientos entre ambos tipos de coordenadas. La energía potencial. ti, de un cierto sistema podrá escribirse como: • ¶1r() • ¶i,Uq,.}) (l.2.i> Donde representa el conjunto de coordenadas flexibles del sistema y representa las coordefladeas rígidas. La contribución ¶1< podrá escribirse como función de las coordenadas 22 5. l,rrooucclo,. flexibles con una dependencia parámetrica de los valores medios de equilibrio de las coordenadas rígidas. a V~((qfl > (1.2.2) Desde el punto de vista dinámico pueden nistir dos maneras de fijar los valores de algunas de las coordenadas internas~’ m. La primera de ellas dará lugar a lo que en este trabajo se denominará modelo de ligaduras flexibles Desde un punto de vista formal la manera de fijar el valor de cada una de las coordenadas rígidas, q,> puede plantearse, por ejemplo, introduciendo un cierto potencial armónico alrededor de su posición de equIlibrio. q%: 1 (1.2.3)2 ><~~ En el limite cuando las constantes de fuerza. K, 1 tienden a infinito, la distribución de probabilidad de la coordenada q,1 depende únicamente de la contribución energética por tanto el valor de la coordenada q~1 será q~,. o lía p(q51> — 5(q,,—q,<) (1.2.4> 1<... 23 La probabilidad de una cierta configuración del sistema. referida al conjunto de coordenadas flexibles, será por tanto, en el colectivo canónico: o ,,~znp ~—~¶l~(q,.q~i1~G(q..qr3~ pÍ. La matriz G esta relacionada con el caabio de coordenadas cartesianas por coordenadas generalizadas en la descripción del sistema (ver anexo C). El segundo procedimiento de fijar los valores de las coordenadas rígidas consiste en introducir en el sistema ligaduras que mantienen los valores de las coordenadas oen los valores de equilibrio. q,. Este desarrollo conduce a los modelos que se denominarán modelos de ligaduras rígidas. Desde el punto de vista físico en este método se reduce el número de grados de libertad del sistema. Planteando la resolución de la. ecuaciones del movimiento del sistema en coordenadas ‘4generalizadas se tendría: fl(p.q) e t¿(q> • R(p.q> (1.2.6) 24 1 ,.rTmotm,ccloN o Los valores q~ son los valores de equilibrio de las coordenadas rígidas, y p~ son sus correspondientes momentos conjugados, que para el sistema con ligaduras rígidas son iguales 4,24.26a cero . Las ecuaciones del movimiento serian: j, ( 83< ] (1.2.8) • 0 (1.2.9) La eliminación de los momentos conjugados a las coordenadas ligadas provoca, cuando coexisten coordenadas internas rígidas y flexibles (ver anexo C). una variación en la probabilidad de las configuraciones, para el modelo de lIgaduras rígidas se obtiene: p4q,) £ s/~• (1.2.10>f dq, expI%8ti,(q~.4’3I la ~q,,q4I 4.16 Siendo O (4<.q~) una matriz cuadrada de dimensiones N~xM~. En el anexo C se describe un posible procedimIento para evaluar el valor del determinante de esta matriz. 25 992 FiCURA 1.2.1 COORDENADAS GENERAlIZADAS DE UNA CADENA LINEAL 4 ~ y ;.., REPRESENtaN LAS COORDENADAS DE ROTAC~N se han utilizado para la incorporación de los efectos cuánticos a algunos de los grados de libertad de un cierto sistema, para los que no es adecuada la descripción cíásicatZ. Por otra parte se han desarrollado procedimientos de simulación (flétodo de Car y ParrineIio 2SO¿) en los que no se introduce un potencial de interacción intermolecular. En estos métodos se consideran los sistemas como conjuntos de núcleos atómicos y electrones, el comportamiento de los nucleos se describe desde la mecánica clásica (mediante dinánica molecular), los electrones son descritos mediante la mecánica cuántica (utilizando la aproximación Born-Oppenheimer). La solución del problema cuántico se aborda mediante la teoría de funcionales de densidad33 (dentro del formalismo desarrollado por Hohenberg y Kohnt y Kohn y Sham3m> utilizando la llamada aproxiamcidn de densidad local’3 (LIlA). Dentro de este formalismo la resolución del problema cuántico es válida sólamente para calcular las propiedades del estado electrónico fw,daaental, sin embargo el volumen de cálculo requerido es considerable—ente inferior sí de procedimientos de tipo Itartree-Fock-Roothant ya que se siapilfica notablemente el cálculo de las interacciones de repulsión e intercambio 28 2. >tOOELo5 DE PoTEMcIAL electrónicas, consiguiéndose pese a ello resultados de una calidad a, comparable La clave del método Car—Parrlneiio consiste en el acoplamiento del proceso de solución de los problemas clásico >. Por otra parte la cantidad de cálculos requeridos limita considerablemente la integración durante largos tiempos de las ecuaciones del movimiento de los nucleos, así como el número de partículas a utilizar. En este trabajo se han utilizado potenciales de tipo tradicional, así coma el formalismo de la mecánica estadística 29 clásica para el estudio de diversos sistemas constituidos por moléculas flexibles, en las siguientes secciones de este capitulo se describen los modelos de potencial utilizados para los distintos sistemas estudiados. 30 2. 14000435 05 poincláL. 2.1 MODELO SINeLFICAflO DE NALCANOS En este trabajo se ha utilizado el modelo propuesto por 24.15 Ryckaert y Bellemans para el estudio por simulación de n—alcanos. En este modelo no se consideran explícitamente todos los átomos, empleándose en su lugar una serie de centros de interacción, cada uno de los cuales representa de manera global a los átomos pertenecientes a un grupo metilo (—CH,) o metileno (-CH,-l. La descripción de la estructura geométrica de una de estas moléculas puede realizarse mediante la implesentación de coordenadas internas. De este modo, para una cadena lineal con grupos, la estructura intranolecular quedará perfectamente definida con 3n, -6 coordenadas, agrupadas del siguiente modo n,-l distancias de enlace: Ib,. b, bh .,>. n,—2 ángulos entre enlaces contIguos <61,0 O,.~> y n,—3 ángulos de rotación interna o torsión: 14,.#,....4~,> (figura 1.2.1). En el modelo de Ryckaert y Bellemans se consideran distancias de enlace y ángulos de enlace fijos. b. y 6~ respectivamente. Las interacciones Intratoleculares pueden separarse en dos grupos: interacciones de torsión, utr que son función de los ángulos de torsión • e interacciones de van der Waais, las interacciones de torsión se modelan como: 31 n. —3 (2.1.1) a.’ 5 V~ (4) — 5~ cos (•> (2.1.2) a —o Los coeficientes a> se presentan en la tabla 2.1.1: Tabla 2.1.1 i 0 1 2 3 4 5 a,/k (Kl 1116 —1462 —1578 368 3156 3788 a 1/Ild aol”) 9.279 —12.156 —13.120 3.060 26.241 31.496 t 1< a Constante de Boltzmann Las interacciones de van der Unís se consideran para los pares de grupos de la misma molécula separados por más de tres enlaces. Para modelar esta interacción se utiliza el potencial de Lennard—Jones truncado a una cierta distancia R~. Las interacciones inter,aoleculares se representan de manera análoga sobre los pares en los que cada grupo pertenece a distinta cadena. Para un sistema constituido por ?¿ cadenas estas interacciones serán: 32 8 mt.., n, 7 7 7¾(Ir -r 1 L L L. •I ai jfl kl.4 N~1 • o, 05 U~jrtmZ Z ZE VisiR) m 4c [e.)í2 (0.)’ VisiR) =0 En estas ecuaciones los vectores del tipo posición del grupo j de la cadena a. Los parámetros utilizados se recogen en la tabla 2.1.2: simbolizan la Tabla 2.1.2 b 0/n. 80/gr.dos s/tj .~0> (r/kT)/X o’ /t¿z fi0 lo’ 0.153= 109.47 599.6 72.0 0.3923 2.50 En la figura 2.1.1 se contribuciones intramoleculares para el caso del o—pentano. muestra la influencia de las en la distribución conformacional 33 2. 140051.05 05 PO~fltcI&L (2.1.31 y (Ir -r 1) LI ‘1 ,RSR~ 1 (2.1.4) (2.1.5) R > MC¡ 2~4 ~S4 O £0 ~ •— ~4Q — O £0 lOO lOO 140 4 ‘a ~ RO ‘Os 3 245 Sc “O FiGURA 2.1.1 PROSASIIIDAD CONroR.jCIoÑAi. PM~A EL N.PENTAflO (7.720<> ‘a 34 2. IOOfl.O5 ot ponOOcIA,. 22 MODELO SINCLIflOADO DE 1.2-DiCLOROE’TANO El modelo utilizado para el 1.2—dicloroetano ~tr + JaLar (2.2.1> >j U 99 utOr= Vteú(t, (2.2.2> >=1 3 tor*V (•l a tal cos (4) <2.2.3) —o 35 M——l 8. a imtn V V r st. — L L L LLS (Ir —r ~flt~ 1k así boa, 5 1.1 k—I SI.—’ NS n~ 5. t ntr rI¼QSLLLit afl b..i J.I k..i y (Ir —t 1. 00 1 bk 0 (R.qj.q,) — Y (8. — o (2.2.5) 8< R~ ~ (2.2.6) .Rt 85 3 V¶R, q 1. q~) ~í q1 q~ fi (2.2.7) N• representa el número de moléculas del sistema y n5 el número de grupos por molécula (des grupos Cl y dos grupos —CH2—). En las tablas 2.2.1. y 2.2.2 se recogen los parámetros utilizados: Tabla 2.2.1 1 0 1 2 3 a1/(kj mol’) 8.021 —12.067 1.393 21.481 bc £~d/nm ben ~CJ¡”~ O~j~g¡~~/grados 1530 0,1785 109.47 36 (2.2.4) 2. $005105 05 POtENCIAl. Tabla 2.2.2 t 1 t 3 ¡ CcN,c$ Lele, C 052c1 dcS2cN2 0cíe¡ acs~c’ R 5/nm qc~ ¿e q~1/e 0.4779 1.1785 0.7404 0.3982 0.3635 0.3805 1.018 0.25 —0.25 1 < Id aol’ 1 t(nm) e u —carga del electrón 37 38 2. tofoS Ot POTENCIAL 2.3 MOCELO REALISTA DE NEUTANO SS Este modelo está basado en el propuesto por lillo y Yip , en este cato no se fija ninguna de las coordenadas Internas, incorporándose los hidrógenos de manera explícita. La energía intramolecuiar se define para cada molécula cono: izara ten Lles tem—fies ten Lles—ter. vdsu —u tu tu tu tu <2. 3. 1) La contribución de tensiones de enlace, UtF£. se representa como: Ut~0= ~{~[E K 11 (R,1—R~1 í~ 3 <13> <2.3.2) Donde el susatorio se realisa (tres enlaces C—C y diez C—H>. El adopta la forna: sobre pares de átonos enlazados término debido a las flexiones 1 ~ g¡lk tít 411k> (2.3.3) En este caso el asatorio se realiza sobre tríos de átomos unidos sediante enlaces L 1 1 • cos 1 3 •cecc) 1 (2.3.4) ccc Siendo *ccc el ángulo de torsión que definen los cuatro carbonos presentes en la molécula. La energía de acopiamiento tenslón-f lesión. ~ se representa como: <2.3.51 Los sumatorios se reali:an sobre los siseos términos que en la contribución U”~’. Los términos de acopiamiento flexión torsión, u””~. Vienen dados por: U •o.- (aík, 1 <2.3.6> Donde el sumatorio se realiza sobre grupos de cuatro átomos en los que se dan secuencias de enlaces representables coso i—j—k—l. Para este tipo de interacción habrá una contribución C—C—C-C. 10 C—C—C—H y 16 H—C-C—M. La interacción se define entre átomos separados por tres o más enlaces, utilizando un potencial Lennard—iones 9—6. 40 2 $005105 fil ~oflflc1AL 7 t”’( I~í—~íI o’ . £ ) <2.3.7) L ti ti R, (2.3.8) Las interacciones intereoleculares se establecen entre pares de átomos pertenecientes a distintas moléculas mediante un potencial Lennard—Jones 9—6. utilizándose los mismos parámetros, c y o’, que para la interacción lntraaoiecular de van der Unís. Los parámetros utilizados se recogen en las tabla 2.3.1. 41 Tabla 2.3.1 Tensiones K 11/(kJ aol”A 2 ) R%/ A C-C 923 1.540 C—H 2397 1.099 Flexiones 1< / C—C-C—C 4.858 C—C—C—H 0.000 0.000 Tensión— X.~í/(kJ mol’A’’) H~¡ A K’ /. en las que el tamaño limitado del sisteaa simulado impone una serie de limites en las fluctuaciones de las propiedades del sisten 3’0 de tal manera que tales condiciones no son accesibles de un modo directo mediante los métodos de simulación habituales, aunque en algunos casos pueden desarrollarse métodos alternativos para solventar estas dificultades’~41. En la simulación de sisteeas homogéneos es preciso ajustar las condiciones de contorno de manera qu.e todos los puntos en el sistema sean equivalentes, para ello se recurre habitualmente a la 4inclusión de condiciones periódicas de contorno De este modo se evitan las inhomogeneidades que surgirían por la presencia de una superficie que limitase las posiciones que pudieran adoptar las partículas. Algunos autores han utilizado, sin embargo, otros métodos para preservar la homogeneidad del sistema, simulando sistemas tridimensionales sobre la hipersuperflcie de una hlperesfera tetradimensionaí’2 41• La peridiocidad del sistema se introduce de manera análoga a como se hace en la cristograíoíia~~45. La periodicidad de la red se expresa matématlcamente como: n pueden construirse a partir de la posición de una partícula las coordenadas de todas sus Imágenes o réplicas en el sistema periódico. A partir de los parámetros de la red puede construirse la llamada celdIlla prImitIva de .4 Ulgner-Seitz , que contendrá una Imagen de cada una de las partículas del síatesa. Las condIcIones periódicas de contorno más frecuentemente utilizadas en sinulación son las que se corresponden con redes de .4 fravais cúbicas . y particularmente la red cúbica simple. Para ésta la celdilla primitiva’de Wigner—Seitz es un cubo, lo que facilita cíe manera considerable los cálculos de distancias y de posiciones de imágenes periódicas, por ejemplo las componentes del vector que une una cierta partícula 1 con la réplica de la partícula j más cercana podría calcularse como: — x, - d anint< t representa las coordenadas absolutas de las particulas. d es la longitud de la celda y anint es una función de que cálcula el numero entero más próximo a un cierto argumento real. 47 Otras condiciones de contorno utilizadas en ciertas ocasiones son las correspondientes a las redes de firavais cúbica centrada en las caras y cúbica centrada en eí interior, que dan lugar respectivamente a celdillas primitivas con geometría de dodecaedro rómbico y octaedro truncado’. Existen ciertos procedimientos de simulación en los que la celda unidad primitiva es de tipo triclínico, con parámetros geométricos que varian a lo largo de la simulación (método de Parrinello y Rahm,sn’ 6). su aplicación fundamental se encuentra en la obtención de estructuras de sistemas en estado sólido. tos métodos de simulación que se aplican en el estudio de fluidos moléculares son de dom tipos, por una parte los métodos de Monte Carlo, que generan una serie de puntos en el espacio fásico del sistema R3. Este hecho no afecta a la estructura del fluido, que para las distancias 3, usadas habitualmente es poco importante. Las propiedades termodináalcas pueden corregirse adecuadaaente. asumiendo un cierto co,portaaiento para la función de distribución radial, g(r). para distancias intermoleculares mayores que R~. De este modo se reduce notablemente el volumen de cálculo a efectuar en la simulación. Sin embargo, para potenciales de largo alcance (interacciones entre cargas. carga—dipolo. dipolo—dipolo. etc 1 a veces no es adecuado utilizar potenciales truncados, en estos casos suele ser conveniente calcular las interacciones con todas las imágenes periódicas (sumas de toíaldt o utilizar algún método alternativo para minimizar el efecto del truncamiento sobre las propiedades que dependen de la estructura a distancias largas (por ejemplo. métodos del campo de reacción medlo’~ 1. Las simulaciones por ordenador pueden realizarme en ‘o diferentes colectivos . Habitualmente las sisulaclones por Monte Carlo suelen llevarme a cabo en el colectivo canónico . integrable analíticamente. La ecuación (3.1.1> pued. reescribirse coaot %‘(1) — “ (3.1.3) Jx¿S dx La estimación del valor medio de A podrá realizarse mediante la ecuación: Sp E A(x±)p — hm (3.1.4> Np E p(x,l/p0(s~l 5—’ donde los suntorio. se realizan sobre N~ puntos pertenecientes al dominio Y elegidos aleatoriamente con probabilidad dada por la función Po En medios densos el procedimiento anterior resulta poco útiiS La presencia de las interacciones entre pare. de centros 52 3. 75c51cA5 ~ 51U54C10S restringe considerablemente la zona del espacio de tases representativa, sólo una pequeflisima parte del bipervolumen Y contribuy, de anera relevante a las propiedades de equilibrio. Es por ello que resulta mucho más conveniente recurrir a procedimientos dinámicos. En éstos mediante un cierto algoritmo se genere una secuencie o trayectoria de confIguraciones de tal manera que en el limite de trayectorias suficientemente largas las configuraciones aparecen con una cierta probabilidad~, PO~ A diferencia de lo que sucedía en el caso anterior la función PO no está sometida a restrlociones especiales. por otra parte la. conf iguraciones generadas sucesivamente no Son independientes, lo que ha de ser tenido en cuenta en el cálculo de errores 4. Las técnicas de Monte Carla dinámico se basan en el formalismo de las Cadenas de MarRov’~’5. Las cadenas de Markov describen la evolución de la distribución de probabilidad de los estados de un cierto sistema mediante un procedimiento de tiempos discretos5m. Supongan un cierto sistema que puede adoptar una serie de estados y considérese que en un cierto paso de tiempo t~ la probabilidad de los estados viene dada por un. cierta función p(t, > normalizad.: ~:. P2. <3.1.5> Ph} Pi ~ 0. Vk 53 1< —~ ~. 1.6) ini es la probabilidad del estado A en el pato de tiempo 1, It es el número de estados. Consideremos un mecanismo de evolución tal que la distribución p en ci paso de tiempo t 1,, dependa exclusivamente la función p en el paso de tiempo inmediatamente anterior, y no de las distribuciones previas, y en el que la evolución viene descrita por las ecuaciones: PUlZlTxí p~ (3.1.7) I—1 La condiciones impuestas en las ecuaciones <3.1.51 y <3.1.63 restringen los valores de los coeficientes ll~, debiendo estos cumplir: E ~ l .1—1.2 a 1 (3,18) k O J 1.. ecuación (3.>.?) puede expresarse en forma matricial: 54 1. troncas 05 SIMILACIOS Ifl Pi 5.’ 02 ~1 A ‘4’ • . .n.I lite ..Ala. nl 1 2121 ITt. ITt. 1~ ~Sm ‘‘5 1 p -hp Aplicando varias veces la matriz II: “1 1o —fl P Pl P 93 Dm (3. 1.9) (3.1.10) <3.1.11) Las matrices IT. que verifican las condiciones expresadas en la expresión (3.1.8) se denominan matrices eetocásticas~ 5~~. El teorema de Perron—Frobenius54~~ establece que los valores propios. X, de estas matrices verifican: xii SI 3 j. Aj — 1 (3.1.12) (3.1.13) rs posible desarrollar una cierta función Olineal de los vectores propios de la o tal., vectores propios, y sea # distribución de probebilidad. expresable cOSO: o como combinación matriz II. Sean a cierta función 55 t—1 Los vectores propios habrán de verificar: Zic&s • ½ Z’ckí (3.1.15) ja, Jal Donde 5kJ representa la componente ,I del autovector k. Teniendo en cuenta la ecuación (3. 1.8) se habrá de cumplir: * ] —> — 0 (3.1.161 1~l Por tanto para que la distribución inicial y las sucesivas verifiquen la condición de normalización se cumplirá que la suma de los coeficientes de los autovectores con autovalor unidad es igual a 1. Aplicando la condición de evolución podrán obtenerse las 1 2sucesivas funciones de probabilidad o . • etc. > 2~ ~ c~ Xj ít~ <5. i. Vn sai Ial 56 ,. ?ttmlc.s br s,íeAsc!c. Al cebo de sucesivas aplicaciones de la matriz U se obtendría: a 4 4 4 (3. i. 18) ‘a, En limite de infinitas aplicaciones sucesivas de la aatriz II sólasente los vectores propios cuyo valar propio cuapla: ¡A~1aI contribuirán a la función p . Se dice que una cadena de Markov es 55irreducible cuando para cada par de estados 1<. 1 existe un cierto 1. tal que Il~~ y ~hk son mayores que cero, lina cadena de Markov irreducible que presenta un único valor propio Ab. que 55verifica ¡Ah Ial se dice que es ergódica . El valor propio correspondiente es, lógicamente Akal. De acuerdo con la ecuación <3.1.18) el limite de la distribución p vendrá definido por: l~ Do - PO a 4 (3.1.19) La rapides con la que se alcance la distribución PO dependerá del resto de autovalores. así como de la distribución inicial. 5? SS 3. Tscrnc.s DÉ sIIlutác,c* 3.1.1 ALGcAITMOS DE TIPO t el conjunto de coordenadas cartesianas de las partículas de un cierto sistema tridimensional con condiciones periódicas de contorno en el paso de tiempo t. Para generar la configuración del sistema en el paso de tiempo t.l se construye una configuración de prueba mediante la variación de las coordenadas de una de las moléculas en un cierto intervalo alrededor de su valor en X~. De este modo se genera la configuración de prueba X~. Por ejemplo, para un sistema tridimensional, moviendo una cierta partícula 1. elegida al azar, se obtiene la configuración de prueba: 59 41 — rtj . vi e i (3.1.1.1)XLI — X~j ¿x } Yt, — nr ¿y a z.~+ Iz Cosi los desplszaaiemtos ¿x. ¿y, ¿a elegidos aiestoriaente, con distribución uniforme en el intervalo (—AAAI. En función del procedimiento utilizado para generar la configuración de prueba se puede definir una función de probabilidad a priori. «(XL4XL) que representa la probabilidad de elegir X~ coao configuración de prueba a partir de X~. El método de )4onte Carlo establece además un criterio para aceptar la configuración de prueba XL como nueva configuración. XL,, • de la cadena de Markov. La configuración será aceptada con una cierta probabilidad” 555. A(Xt 2(3, tal que en el limite de infinitos pasos de simulación cada configuración aparezca con su probabilidad p(X). La probabilidad de un tránsito desde la configuración X. a la configuración 1b’ vendrá dada. en función de las probabilidades ft y A como: a I1 a Q <3.1.1.4) La ecuación (3.1.1.4> se cuaplirá si todos los sumandos del primer alembro son nulos: Utilizando la ecuación (3.1.1.2): A(Xh4X,1 ~ <3.1.i.61 AtXa4Xb) -P~ • p(X~) t(Xa..X,) ] (3.1.1.9) En los procedimientos más habituales de simulación, se utilizan probabilidades a priori, a. que a ~<~b.~e) (3.1.1.10) Este es el caso del movimiento descrito en las ecuación (1.1.1.1). La relación de probabilidades de aceptación vendrá dada como: A(X,4Xb) P(Xb) A(Xb4X.l p(X,) (3.l.l.iil La ecuación (3.1.1.5) representa un criterio suficiente, pero no necesario, para conseguir que el muestreo efectuado en la simulación conduzca a la distribución de probabilidad de equilibrio. p. siempre y cuando la cadena de Markov conseguida sea ergódica. es sin duda una condición más estricta que la expuesta en la ecuación (3.i.i.3l; sin embargo resulta muy práctica su 62 2. TUSIC*5 rs SnstAclou utilización ya que permite calcular de rodo sencillo los criterios de aceptación de configuraciones de prueba. La citada elección suele denominarse coso condición de reversibilidad microscópica o de balance detaiiado’”. La eficiencia de un cierto .lgoritao de simulación por Monte Carlo depende de tse modo funasental de 1. elección de las funciones a. Para sistemas moleculares flexibles resulta, por ejemplo, más adecuado utilizar coordenadas generalizadas la eficiencia del método de simulación 4dependerá del valor atulmo del de.plazaailento . A. El valor óptimo de A será función del estado termodináauico. en particular de la densidad del sistema. Puede decirte que la eficiencia de un muestreo por Monte Carlo dependerá de la cantidad de información, que sobre el espacio de fases configuracional. se introduaca en la construcción de las funciones a. Como ejemplos extremos consideres. primeramente el caso en el que la configuración de prueba se elige al azar sobre el espacio de fases del sistema. Si se parte de un estado de energía baja, y dado que sólo un parte muy pequefia del volumen físico será representativa del tistes en las condiciones 65 de temperatura y densidad propias del estado liquido, prácticaaente todas las configuraciones de prueba que se generen serán rechazadas, ya que todos los puntos del espacio físico tendrán la misma probabilidad de ser elegidos como configuración de prueba. Partiendo de un estado de alta energía la evolución hacia configuraciones de baja energía será, por las mismas razones, conmiderablemente más lenta a medida que el sistema desciende hacia los estados enérgeticamente favorecidos de acuerdo con las condiciones termodinámicas. En el otro extremo estaría un hipotético algoritmo de simulación tal que la función a estuviese definida como: opLa — a al azar con probabilidad uniforme en un intervalo alrededor de los valores originales <#.O.*>: • a •+¿~ ,5~ • I~-At.*+A41 1 O — ja • s.¿o .38 e Iw-AS,6.&91 a $•¿~ e <~‘¿*,*~A~I J Donde J es un núaero entero tal que el valor de O permanece acotado en el intervalo Lot>. Para un ,aoviaiento puramente oriantacional desde la configuración X a la configuración de prueba X por variación de los ángulos de Euler de una de las moléculas seglJn el anterior algoritmo, el criterio de aceptación, teniendo en cuenta el jacobiano de la transformación de coordenadas cartesianas atómicas 66 ~. TÉOICÁ5 Di SINIA.At104 a coordenadas generelizadas y la slmetrla de las funciones a. adoptaría la forma (véase anexo fl): A(X41) p(X~) exp(—$1i(X)1 serde) ________ a a <3.1.1.16> A cuando 840. se recurre a sustituir la forma del muestreo, remplazando la variable 0. y utilizando en su lugar cosO como coordenada .éoleculart En tal caso la forma de originar configuraciones orientacionales de prueba seria, representando cose con el símbolo c y coso como c c — 2j.c.¿c ,¿c e (—Ac,ÁcI • — •.A• 3~ e <—t•.A•1 1 <3.1.1.17) • — *‘~* ,¿~ e I—A@,~1 j Donde j es un número entero tal que c permanece acotado en el intervalo 1—1.11. Las probabilIdades de aceptación serán en este caso: A(X.eX) P< 1b> exp¡-tU(X)I _______ — — (3.i.1.is) 67 Este último método propuesto para el muestreo de las orientaciones toleculares puede plantearse formalmente de una manera diferente, aunque equivalente en la práctica, que consiste en seleccionar el valor o en un intervalo de dimensiones variables con probabilidad no uniforme, las ecuaciones correspondientes son: sen(O 1 a(O-.O 1 — (3.1.1.19) eche1 sen<,~) dii 9t0 Las variaciones má%isas , te’ y 1.6 se definen sediante: OCAS oc — f sen <~) dii.. cose — cosfe.A6) (3.1.1.20) e o te af sen d~s cos(B—06) - cos y velocidades IR,> en función del tiempo. t. La integración del anterior sistema de 3M ecuaciones diferenciales de orden 2 requiere la introducción de 614 constantes, que definen el estado del sistema en la situación inicial. La resolución de <3.2.1> requiere, salvo para casos excepcionalmente simples, la utilización de métodos numéricos. Los métodos usados habitualmente 73 (3.2.2) físicas (2.2.3) se basan en la aproximación de diferencias finitas~. Desarrollando en serie de Taylor el desplazamiento de una cierta partícula en función del tiempo: a R, at • ~ a, it) st. <3.2.5) 74 3. ttcWiC&s DÉ 51QJ1.AcIOm II, It—St) = R, It) — y. (Zótí + eat 3 í Ile modo análogo: y, It) y, v, se obtiene uno las ecuaciones del 1 <3.2.11) La expresión (3.2.11) representa el algoritmo de Verlet de integración de las ecuaciones del movimiento en la llamada versión (3,2.6) (3.2.7) (3.2.8) <3.2.9) 13.2.10) 75 leap frog . Este algoritmo permite una buena conservación de las 4 •constantes del movimiento (energía total y momento lineal) . es reversible, y fácil de adaptar a sItuaciones en las que existen 17,71 ligaduras rígidas en el sistema . Las simulaciones por Dinámica Molecular llevadas a cabo en este trabajo fueron realizadas en el colectivo microcanónico (Número de partículas. vcluaen y energía total constantes> utilizando el algoritmo leap frog. 76 1. tÉCE~.5 DÉ SIXuLACION 3.2.1 SISTEMAS CU’J L>GADLSAS Los modelos moleculares empleados en cate trabajo contienen ligaduras holonómicas’”’ 2’ que restringen el movimiento de ciertas coordenadas internas de las moléculas. El tratamiento de estas ligaduras en las técnicas de simulación mediante Dinámica Molecular se realiza habitualmente considerándolas ligaduras rígidas’5’21. Como se discutió en secciones previas, la utilización de ligaduras rígidas o ligaduras flexibles puede dar lugar a distribuciones de equilibrio diferentes en las propiedades configuracionales del sistema cuando en las moléculas coexisten grados de libertad Internos ligados y flexibles. En este trabajo se han simulado mediante Dinámica Molecular cadenas de hidrocarburos considerando distancias de enlace y ángulos de flexión rígidos. Estas ligaduras pueden describirse matemáticamente cOmo: ji«ifia «~~ 1 d~ .0 .k1 8 (3.2.1.1> R nR -R (3.2.1.21 Pi a1 Siendo R el vector que une las particulaa enlazadas mediante la ligadura ¡e, y 14 el número total de ligaduras. 77 Derivando respecto del tiempo, una ecuación del tipo de las anteriores se transforma en: R -R —o <:3.2.1.31 «kOt «tA& •R« 8+R ~ Rak$k O (3.2.1.4> £1 significado físico de estas ecuaciones se puede expresar diciendo que la velocidad relativa de dos centros conectados mediante una ligadura es perpendicular al vector que los une, y por otra parte que la aceleración relativa tI~j es tal que compensa la aceleración centrífuga generada por la rotación del enlace entre los centros 1 y J. Las ecuaciones del movimiento para un sistema con ligaduras rígidas pueden expresarse de acuerdo con las ecuaciones anteriores, y según el algoritmo de integración de las ecuaciones del movimiento anteriormente definido como: R,tt+Atí — Ijt) • tt Rjt+Lt/2> (3.2.1.5) a a, < f, + G, 1 (3.2.1.61 (3.2.1.71 78 ‘. rrcaíc.s DÉ s¡wjtaclow Donde el termino O. representa las fuersas debidas a la ligadura, es decir las fuerzas a añadir pata que la trayectoria del sistema mantenga las restricciones expresadas en las ecuacIones <3.1.2.3) y (3.1.2.4). Por otra parte las fuerzas Chan de ser normales a los desplazamientos del sistema al objeto de que no introduzcan variaciones en la energía total del sistema. Por tanto las fuerzas de ligadura actuaran en la dirección de los enlaces R, ~. Lean los vectores entre los «3~ t 2B2 centros conectados mediante ligaduras rígidas, tiendo M el numero de ligaduras holonómicas del sistema. Las fuerzas asociadas serán del tipo (véase anexo El: a,~, «~$~ Los parámetros X~ serán tales que me verifique la ecuación 12.2.1.4). Las aceleraciones É pueden expresarse en función de ej ~ las fuerzas presentes en el sistema 1< -1 fl ~a F —m F + 1 A1fl Cjj 13.2.2,91 si «~ Donde los coeficientes C1¡ vienen dados por: ej ~¡ej 5a 1aj }C1,= { n~’ [fi —5 1 — Ya’ [5 1 (3.2.1W) 79 Las ecuaciones 13.2.1.4) y <3.2.1.9) permiten construir un sistema de II ecuaciones acopladas con ti incognítas. cuya resolución permitiría obtener los valores de los coeficientes A¡. 6 —, —lmr -m r-Ii 1 }«í~i ajflj $j aj a~ 1=1 <3.2.1.11) Se han desarrollado algunos algoritmos de s1mulación~” para sistemas con ligaduras que implementan la resolución del anterior sistema de ecuaciones para obtener los valores A¡ . Sin embargo a menudo presentan inconvenientes de tipo numérico, ya que tales algoritmos garantizan el cumplimiento de la ecuación (3.2.1.4). sin embargo la acumulación de errores de redondeo puede dar lugar a desviaciones en la ecuación <3.2.1.3) que ocasionarían al cabo de un cierto tiempo de integración de las ecuaciones del movimiento desviaciones apreciables de las distancias intramoleculares sometidas a ligadura de sus valores En este trabajo se ha utilizado para mantener las ligaduras del sistema el procedimiento denominado SHAXE””’’ 2, que evita las desviaciones comentadas anteriormente. La aplicación del citado método depende del algoritmo utilizado para resolver las ecuaciones del movimiento, el procedimiento consiste en calcular las fuerzas de ligadura de manera que para cada paso de tiempo en so 2. yzcflcas os slmA...cloN el que se especifican las posiciones atómicas se verifiquen las relaciones expresadas en la ecuación (3.2.1.1>. Para el algoritmo de integración ieap frog se llevaría a cabo el siguiente desarrollo: fi, (t.At) fi, <3.2.1.15)a. Siendo R~ + AtamlZ A,R(t) la —5. I.$í ~.«, 1=’ 3.2.1.17> Donde los términos a son deltas de Kronecker. vectores entre átomos ligados serán: Ji ji ~ Los nuevos (3.2.1.183 I=1 Operando se obtiene: ji tt+ttiR (t+ótl = jiO ~ ~«ial YE +átau:$(t.ÁtíZA,n ~(ti C>~ * I—1 +óO ~ 7j A~ A1 C~6 C~í R«60~t) It (tI Ial <3.2.1.19> Las ecuaciones para resolver los coeficientes A serán las (3.2.1.17> y <3.2.1.19). El método SRA/CE se basa en una solución 82 3. ltcJflc*5 Os slm>Laclo» num4rica de tipo iterativo. El procedimiento cm el siguiente: a partir de una cierta solución Inicial [6 —S 1 ‘1$, l.«j (3.2.1.20) De este modo se obtiene: 16 ji 1 La condición <3.2.1.19> puede expresarse como: It (t+At)R (t+At> = U’ (t+At>R’ se construye un nuevo conjunto. 1 1 1{A,,A A.,>, en el que se variará uno de los coeficientes, por ejemplo A,. 1 0 A, A¡ , 1 5 1 <3.2.1.24) Sustituyendo estos valores en la ecuación <3.2.1.22) asociada con la ligadura 1. y considerando la. definición de los coeficientes Cj R (t+At)R (t+Atl = fi’ (t+tt).Ru R asociada con la ligadura establecer el nuevo valor A 1: 1 0 1 A, — A, + f -‘ ‘1 a 4. la, a,> ¡ d~ — fi <5 ).R 1 0(t+ht> 3 A partir del nuevo valor A, puede obtenerse una nueva serie de posiciones aproximadas de las partículas: 1 0 2 EstasR.4}. , 1posiciones serán idénticas a las (R~, U2,..., fi.,>, excepto para los dos átomos implicados en la ligadura ‘1’. para los cuales se ha cambiado el valor A, apretinado. Las nuevas coordenadas serán: = ~‘ (t+At) — At 2m¾)4—A~) ji (t) a, - a, «1 U5 + AY.2’ ‘—4) U‘aa, anterior, la 1’ permite <3.2. 1.26) SS partículas «~ y/o 3,. cuyas posiciones han sido alteradas. A partir del nuevo conjunto de valores, fi 2 se puede obtener un «101 nuevo conjunto de valores A: <4,4,... 4>, que difiere del conjunto <4.4,4> en el valor del coeficiente A,. El valor mejorado de éste, se obtendrá de modo análogo a como se obtuvo el valor sejorado 4: 4 — U’ A 2 — A2 + [ 1 <3.2.1.29> 1. - 2 ¡2 —I — ti (GR (tv-ót)At a 4. ~ a,0, «,3, «2 A partir de este valor se podrán calcular los vectores U, y 1 3ji • con los cuales podrá estimarse un valor mejorado A,, etc. 5’ . U, Este proceso se repite ciclicameate sobre todas las ligaduras hasta que se verifica el criterio de convergencia (3.2.1.23) para todas y cada una de las ligaduras. La velocidad de convergencia del método depende de diversos factores’: paso de tiempo de la simulación, solución inicial, grado de acoplamiento entre las distintas ligaduras, etc. En los sistemas estudiados en este trabajo se utilizó se utilizó como solución inicial: oA1 — O . lal,2 14 (3.21.30) £6 3. rEcHicas 05 SlwtJtAC1ON El valor de la tojerancÍa, o, utilizado en el criterio de convergencia < ecuación (3.2.1.23) 1 fue de lOa, En los casos estudiados en este trabajo no se observaron problemas de lentitud de convergencia, ello pudo ser debido a que los pasos de tiempo utilIzados fueron relativamente cortos < del orden de 2.201% En aplicaciones del método SI/ARE o similares” en la simulación de 74moléculas rígidas donde se utilizan pasos de integración del orden de 10’ ‘a. es conveniente optimizar la elección de la solución inicial de loa coeficientes A con el fin de que el tiempo de cálculo requerido en el reajuste de las posiciones de las partículas no redunde en una disminución de la eficiencia obtenida 75al utilizar pasos de integración relativamente grandes 87 £8 a. TÉcnicas DÉ sI~u.acYoN 3.3 ANALISIS DE CONFUGIRACIONES Según se comentó en la sección 3, los procedimientos de simulación permiten obtener un conjunto de configuraciones en un cierto espacio fásico de tal manera que en el limite de gran nunero de configuraciones éstas aparecen de acuerdo con una cierta probabilidad p(p,q). La utili&d d~. .~ .——“ ~i’n,,1ar¶ón no se reduce, sin embargo, al cálculo de las propiedades en unas ciertas condiciones físicas . En otros casos surgen problemas de ergodicldad al emplear una cierta función de probabilidad p o para el número de configuraciones (Netropolis Monte Cario> de la simulación. Resulta útil en estas situaciones reducir esos impedimentos mediante la simulación en unas condiciones más favorables, utilizando una cierta función densidad-de probabilidad PO en lugar de la probabilidad p. Las funciones p y Po pueden diferir en condiciones físicas — <3.3.1) U Transformando la ecuación anterior se obtiene: 90 3. rEcalcas LE s,nuLAc2on J dR PO<~> [p(R>/p, a 1 3.3.2) dR po s (3.3.3>~ Explícitamente pueden escribirse los valores medios obtenidos para los sistemas p y p0 a partir de una simulación del sistema p~ como: a A(k) <3.3.4) kC PO 2 A(k> kE(ji> PO — (3.3. Sl p(k > Z Pu representa un conjunto de Po configuraciones obtenidas en la simulación par. las condiciones p~. siendo N, el número total de tales configuraciones. Obviamente la ecuación <3.3.5) se reduce a la (3.3.4) cuando p y PO son la misma función. Definiendo los pesos estadisticos 5k como: 91 p(k) U W = <3.3.6) L p, MC) Po Puede reescribirse la ecuación (3.3.5) de manera mas compacta: paZ A(kl w~ (3.3.7) kdR> PO Ya que se verifica: Z 4k — ¶ (3.3.2> he p La ecuación (3.3.7) expresa el sentido matemático de la transformación entre sistemas. Al pasar de unas ciertas condiciones a otras varia el peso estadístico de las conliguracicnea. Es por tanto necesario, para realizar una evaluación adecuada de las propiedades, disponer de un amplio conjunto de configuraciones del sistema simulado para q,le aparezcan, aunque sea con pequeña probabilidad, las .ronfigursclnnen más representativas del segundo sisteaa. Al objeto de analizar la calidad del método de ,muestreo para condiciones distintas a las de la simulación se utilizará un parámetro de eficiencia estadística ~: 92 3. TSC>IIC*5 DE 51MU~sCOoN La — <3.3.9) 2 2 2 > N~< W~ > Para el sistema de referencia se verificará: 1/Id, V k, ca 2, 2 <3.3. 20) SaI (3.3.11) Si sólo una de las configuraciones de partida tiene importancia en el sistema p: 1. ~ ..sw,.= 0 (3.3.12) £ a 1/1<, (3.2.13) —l El valor de S puede interpretarse, por tanto, como el número de configuraciones que para el sistema p dan una información análoga a cada configuración en el sistema p 0. En este trabajo se analizó, mediante el formalismo anterioraente desarrollado, la influencia de los parámetros del potencial intersolecular para n—slcanos en algunas propiedades termodinámicas. Los resultados obtenidos para la presión y la 93 energía intermoiccuiar con el modelo de potencial utilizado no concordaban con los valores experimentales. A fin de estudiar la dependencia de las propiedades termodinámicas con los parámetros de potencial intermolecular se llevaron a cabo una serie de cálculos que permitieron establecer algunas interesantes conclusiones sobre la forma que habría de tener un modelo de potencial transferible para fluidos constituidos por moléculas de n-alcanos. los detalles del método de cálculo, así como los resultados obtenidos se presentan en la sección 5.3.2. 94 2. TÉcNICAS OC SiMutACiofl 3.4 EST¡MACiON DE ERRORES Mediante métodos de simulación es posible obtener una serie de configuraciones representativas de un cierto sistema. A partir de tales configuraciones se puede obtener información termodinámica, estructural y dinámica. Según se indicó previamente los métodos de Monte Carlo y Dinámica Molecular dan lugar a una cierta trayectoria en ía que las configuraciones sucesivas no son estadisticamente independientes, y por tanto tampoco lo son los valores de las propiedades de interés. La existencia de tales correlaciones entre configuraciones debe, por tanto, ser seriamente tenida en cuenta al evaluar la precisión de los valores medios obtenidos para las distintas propiedades Supóngase que se han obtenido en una cierta trayectoria de simulación Id valores de una cierta propiedad ‘2<’, habiéndose extraído tales valores de configuraciones o conjuntos de configuraciones igualmente espaciados: Y~.) (3.4.1> La varianza de los datos vendrá dada como, 2 2 <3.4.2)e, — <2<> SS Agrupando valores sucesivos de la propiedad 2< podremos definir: <~I~ Y Y.— —.~.— ~ (3.4.3> (3.4.41 Obviamente se cumplirá: (3.4.5) Con los valores y, pueden definirme desviaciones análogas a de orden a: 2o 5 — < y >2 El valor medio de será: N/a m —k—~[ E E E Xm4tai)•k 4<4a014J ) ‘al k1 jal <3.4.6) (3.4.7> (3.4.8) 96 ,. ,Ec,’ycas os slws.ÁcloA La ecuación previa puede reescribirse en forma de promedios: m—I a ..1y4 a < E 2> + 2 (m—k) c 2<, ~ > ] <3.4.9) a.’ Teniendo en cuenta la igualdad <3.4.5): a-t a [~<~> + a <1—2—) F(k) 1 (2. 4.10> tal La función de autocorrelación E se define como: c 2<, X,,~ > — < F(kl a <3. 4.11) 2 2 c 2<, ~ — < De modo semejante a como se define la varianza de la media en situaciones de medidas Independjentes~ pueden definirse los 2 parámetros a~«2<>) en función de m, que estarán relacionados con las varianzas mediantes ae o~(X) a Mediante las ecuaciones (3.4.11) y <3.4.12) se obtiene: 97 — ——o — 2 ~ (l—=~) <3.4.13) a.(.CX>) — ~«2<~ 1E(O) + L m <‘.0 Que puede expresarse de modo aproximado introduciendo una integral como: m — e~(CX>) [2 [F(k) ( — .±.....> ~jc] (3.4.14) 6 En el limite de valores de a altos, y suponiendo que para un cierto valor de 1< 1 k << e 1 la función de autocorrelación E es próxima a O. podrá sustituirse la ecuación <3.4.14) por: — o~(cX>) [21 F(k> dic ] (3.4.15) Si la función de autocorrelación r se comporta bien, lo que sucederá normalmente si se trabaja con trayectorias suficientemente largas se podrá utilizar como estimación adecuada de la desviación estándar de la media: fl ¿«2<>> — 4() [2 F) 41<2<>> .t <3.4.17> Siendo ~ el téraino entre corchetes de la ecuación (3.4161. que se denominará longitud de £neficacie estadlstica o tiempo de rejajacjón. En función de la varianza 4txí se podrá escribir: 2 a N En la ecuación anterior queda patente el sjgnuficado fisico de t El cociente MIt representa el número de configuraciones independientes equivalentes (para la propiedad Xl a la trayectoria de longitud Id obtenida en la situlación. Las barras de error. ±AX. que se dan para los resultados obtenidos en este trabajo se obtuvieron como: £1 a 2 2 <3.4.19>«2<>] 99 4. SIM(LACION POR MONTE CARLO DE MOLECU..AS FLEXIBLES 4.1 SIM4,LACION EFICIENTE DE N-BUTANO La molécula de ro—butano es la pricera de la serie de los n—alcanos que se caracteriza por la existencia de ángulos de rotación Interna. Desde el punto de vista del anflis conforisacional presenta tres confórmeros estables: t, g y g (ver sección 2.1). El estudio de la influencia del medio en el equilibrio conforesacional del n—butano ha atraído la atención de numerosas investigaciones tanto desde el punto de vista teóricos’ e’ como de síaulaciór,1 21&10 ~ Sin embargo, pese a la aparente simplicidad de la molécula, los estudios mediante simulación han ‘eencontrado, hasta fechas • diversas dificultades para obtener conclusiones definitivas. Por una parte, el ya comentado (véase seccion 1.2) problema de la diferencia entre las distintas vías de fijar los grados de libertad internos poco relevantes ha provocado no poca confusión a la hora de analizar los resultados. Por otra parte, la precisión de los resultados obtenidos mediante simulación por Dinámica Molecular está considerabícaente condicionada por los largos tiempos de la relajación de los movimientos de rotación interna • a ello contribuye, no sólo la presencia de importantes barreras de ial energía en el potencial torsional, sino también las interacciones intermoleculares que generan impedimentos estéricos que dificultan la transición entre los distintos confórmeros. Desde el punto de vista de la simulación por métodos estándar de tipo Metropolis Monte Carlo al tratar con modelos moleculares que combinan coordenadas internas rígidas y flexibles se hace preciso el trabajo con coordenadas moleculares (ver sección 3.1.1). Frecuentemente se ha tratado la coordenada de rotación interna de modo análogo a como se hace con las coordenadas cartesianas en simulación de sistemas atómicos” ~ es decir. generando una nuevo ángulo ~ dc prueba uniformemente distribuido alrededor del ángulo ~ de partida. Este tipo d.c procedimientos adolece de las miagas limitaciones que se exponían anteriormente para el caso de la Dinámica Molecular: la presencia de las barreras Serma de manera notoria la eficiencia del método. lográndose pocas transiciones entre los diferentes confórmeros. Para evitar las dificultades Inherentes a la existencia de la barrera algunos autores han recurrido a movimientos internos en 90 los que se generan grandes desplazamientos o incluso el ángulo de prueba e se genera al azar, uniformemente distribuido en todo SIsu intervalo de definicion . e e 10. 2.]. En estos casos gran parte de los intentos de cambio en la conformación son rechazados debido a que la probabilidad de los ángulos de torsión está fuertemente concentrada en las zonas de los mínimos del potencial torsional. 102 4 KOI.ECVISS fl.flICLZS Finalmente, en este breve recorrido por las diferentes maneras de abordar las dificultades que conlleve ia obtención de precisión en el cálculo del equilibrio conformaclonal del n—butano y otras moléculas flexibles de pequefio tamaño, es preciso hacer •0, U • 60, El referencia a las técnicas de tipo umbrella sa,npLing Estas técnicas pueden aplicarse tanto a simulaciones por Monte Cario como por Dinámica Molecular, en ellas en lugar de utilizar el potencial de rotación intern, del modelo Utora en la simulación ter. se emplea otro potencial U 0 Este nuevo potencial permite una mayor facilidad en las transiciones entre conlórnoros. Las tar. conforeaciones obtenidas en la slaulaciM, utIlizando U0 pueden utilizarse. mediante el conveniente repesado (ver sección 32), para calcular las propiedades en el sistem, con energía de torsión - Al utilizar este tipo de métodos es preciso mantener ciertas precauciones. Una variación e1 U 1’ 0 0It O 5 k 1 Donde RU, g~. dR’ y di? engloban las coordenadas dc translación, orientación, así como los correspondientes términos que aparecen en el jacobiano de la transformación de coordenadas atómicas a coordenadas generalizadas (véase anexo Dl. La ecuación (4.1.4) se transforma, mediante el cambio de variables anteriormente comentado, teniendo en cuenta (4.1.3) en: z—f a’ J ~s? .. .{¾Ú~. exp ¡~olP(RM.oN.,’l se sisplifica. obteniéndose. 7 = ~ j dR” dO’ {dui.. .Jd~ exp [—$tt’(R’,It U(léN)fl It’ flM 0 0 (4.1.10) A partir de este desarrollo, puede caracterizarse la técnica empleada para el muestreo de los grados de libertad torsionales en las sisulaciónes de n—butano en fase líquida. El procedialento consiste en generar, para una cierta molécula. la coordenada de prueba o> con probabilidad uniforme en un cierto intervalo alrededor de la coordenada de partida ~. La amplItud de este 107 intervalo es constante para todos los valores de u~. Denotando como tu el máximo desplazamIento en la coordenada u, y teniendo en cuenta el carácter periódico de los ángulos de torsión interna, se pueden definir los intervalos, para la coordenada j, en los que se genera el desplazamiento a partir de úa~ como: De acuerdo con este criterio en los desplazamientos la probabilidad de elección de elegir u 5 como coordenada interna de prueba a partir de u, es igual a la de elegir u. cono coordenada de prueba a partir de — a(u1...a,1 (4.1.12) En la segunda vía para explicar el método de simulación se trabaja directamente con las coordenadas internas ~, en este caso se genera un nuevo ángulo • (de prueba), a partir de un cierto ángulo de partida 4~, en un intervalo alrededor de éste con probabilidad no uniforme. Los limites del intervalo dependen del potancial de referencia V,: • c(~> exp I-0V4•)l. Y4 • I (4.1.13)• O ~•4 los 4. ItoIrctL.s ttflIBLtS En esta ecuación c(4) es simplemente una constante de normalización. La relación entre los dos formaliseos aparece con la definición de los li,,ites de los intervalos. ya que los parámetros t( y tó’ dependen de 4: •0 Aa> = ~ J exp 1—0V~(O)! d4 (4.1.14> 4-t •0 ~AO tu — 4~ J ex~ (-$V4•)3 dO <4.1.15) •0 Mediante estas ecuaciones puede deterainarse el valor de las constantes c(41. Utilizando el criterio de normalización: r «<4o~l d$ — 1 (4.1. i61 se obtiene: 4 .t• c(4) f exp E—$V,(~)1 d• • 1 (4. i. 17) 4 c(4) — 1 / ( 2 ñu 4,.) (4.1.18) 109 La característica más importante de este resultado es que no depende del ángulo t~. sino que es una constante para todo O en el intervalo 10,2.1. Teniendo en cuenta lo ya expuesto en la sección (3.1.1). la condición de reversibilidad microscópica permite calcular la relación que habrán de verificar las probabilidades de aceptación de las transiciones directa e inversa entre dos configuraciones (1 y Ji que se diferencian en el cambio del ángulo de torsión de una cierta solécula. Siguiendo la nomenclatura desarrollada en secciones previas: A,j Pi «ji (4.1.19> 1’ P 1 «±j Para el primer desarrollo (donde se utilizan las coordenadas o» se verifica de acuerdo con las ecuaciones (4.1.10), <4.1.111 y (4.1.12). = = 1 (4.1.20> Pj e%» ~jLJ ~~vtr ~ +SVÁ*(a>ífl 1 (4. 1.21) Pi C%P t¾n~term<41 1> •0V,.(OGa,l> 1 Para el segundo desarrollo las correspondientes relaciones se obtienen haciéndo uso de las ecuaciones (4.1.5). (4.1.13> y (4.1.18): no 4. MOtEOLA5 ,tntecss te,. p 5 exp E -atPuí — pV (4) 1 eXp E —OU’ 4(i) — 0~’~ (o,) (4.1.23) Mediante las ecuaciones anteriores se comprueba que ambos formalismos conducen al mismo criterio de aceptación: A 11 exp ¡ —0J’~(il ~v~” (tjí •PV,Átjl 1 (4.1.24> ( p)JL>(fl ~ tsr. exp y (~,l +0V44,) 1 Para demostrar le equivalencia de ambos formalismos sólo resta comprobar la equivalencia en las funciones «~. Considerando implícitamente la periodicidad de 4 para no complicar el análisis, partiendo de las expresiones (4.1.11) y considerando o1 como una constante, la probabilidad de escoger una cierto valor de prueba de con las 4 (dentro del priser formalismo> o el cálculo de los factores t4 y A* ~en el segundo formalismo) no se adaptan de isodo coneeniante a la eficiencia que se requiere en los cálculos mediante simulación por ordenador. Ello es debido al hecho de que la función ~ (4) en el modelo de potencial utilizado viene dada en función de un desarrollo en potencias del coseno del ángulo Ita 4. sottcuuás FLfllsi.ts ter. torsional y por lo tanto la función expl—$V (4)] no es integrable analíticamente”. Para evitar la complicación de efectuar cálculos númericos cada vez que en la simulación se genera un movimiento torsional se puede elegir entre distintas alternativas; sustituir la función potencial de referencia Y 0. empleando en lugar de usa expresión cuya exponencial de Boltzeann sea integrable analíticamente (en tal caso habría que utilizar como criterio de aceptación la ecuación (4.1.24) en lugar de la más sencIlla (4.1.30)); tabular previamente a efectuar los ciclos de simulación de manera adecuada la función 4=4dw) e interpolar a partir de la tabulación los valores requeridos durante la simulación. En este trabajo se eligió, sin embargo, un método más sencillo. La idea básica es la de sustituir el ángulo de torsión 4 continuo, por un conjunto de ángulos discretos representativos del sistema. De esta manera se facilitan notablemente los cálculos a realizar dentro del ciclo de sisulación. Considerando un número suficiente de puntos la descripción del grado de libertad interno mediante el conjunto de conformeros discretos ea totalmente análoga a la del sistema continuo. Del mismo modo que se hizo anteriormente la discretización puede plantearse desde los do. formalismos. Considerando las variables u, el cabio consiste en sustituir la variación continua de u por una serie de valores discretos uniformemente espaciados. 113 Para un número de confórmeros discretos m representación de esta discretización es: 04 1 0. 1 1 ‘ «.>~. W • i—i/2 . i.’l.2 a } (4.1.31> Cada valor u 1 estará relacionado 4j según la ecuación: 4, i—i/2 — ~ cm [—$V4~fl d~ con un cierto valor discreto (4. 1.32) Considerando implícitamente el carácter variables a>. se puede escribir la probabilidad a 4 1 1 = o .vj.ij—ii ‘Am cíclico de las priori « como: (4.1.33) Desde el punto de vista del segundo formalismo. La discretización puede plantearse como un proceso en el cual se sustituye el dominio de variación continua del ángulo 4 por una serie de intervalos definidos mediante los ángulos auxiliares 7<.: 114 4. uottC’nsm ÉtUIBIn exp [—~V4g4] de— ~;‘ ¡ (4.1.34) 1< • 0.1 e El intervalo k se define como 1~t.,’ nl A partir da estos intervalo se definen los ángulos k, representativos de cada Intervalo COSO: It <4>k = (4.1.35) J;:$XP V$V,(f)1 d~ De este modo se puede establecer una relación entre cada ángulo en el continuo con un ángulo discreto que lo representará en la sisulación. Escribiéndo la transformación realizada como: O = 0(4) (4.1.36)} 0(4) • 4 e E a~., n 1 Donde la función 4 adopta los valores discretos <4><. definidos anteriormente, en función del intervalo (rkl. Tkl en el que está contenido 4 (vémnee figuras 4.1.1—4.1.3). Este segundo formalismo permite entender mejor el significado físico de la aproximación que se realiza al utilizar la 115 representación discreta. Si 4,.4 4 repreSentan los ángulos de torsión de un sistema constituido por N moléculas de m-butano. la aproximación realizada es: 41 — Lt’(t,t,$( 4~) 4<4)) (4.1.37) La bondad de esta aproximación dependerá, obviamente, del número de puntos discretos utilizados. al aumentar éste mejora la concordancia (4(4) • 4). Por otra parte, debido al procedimiento utilizado en la definición de los intervalos, la anchura de éstos es menor (y por tanto la concordancia entre 4 y O es mejor) en las zonas de mayor probabilidad del ángulo torsional para la molécula aislada — 74 j dR j dO j d4’ podrá calcularme como: = z;1 J dt J ~gU frío’ expi—0&’(R’, O”.*”(t’h11 ,. s 5 2 (2 4 U k1 El valor de que se incluye en las dos ecuaciones anteriores se define como: 7~<... J dR’ J di? J d4” exp(~MtJ(RU,Oit*U(•~ 1)3 5 U UA O 4 It k1 En general resulta más sencillo calcular los valores medios de las propiedades estáticas a partir de la ecuación (4.1.39). ya que en la simulación se trabaja directamente con los ángulos 4,. Para profundizar más en el significado físico de la 217 aproximación discreta consideremos dos posibles valores del ángulo torsional: ~, y w2. Si tales ángulos pertenecen ml mismo intervalo de discretizaclón se verificará: (4.1.411 La relación de probabilidades entre ambos ángulos se puede obtener Integrando la función de probabilidad de la aproximación discreta sobre las variables externas y sobre ifrí ángulos de torsión: ~(flUgU,U~t45) = exp 5—, —I ter. Integrando sobre todas las coordenadas externas, y sobre las coordenadas internas de 8—1 moléculas se obtiene: Donde ~;(*(4,)) es: z;t*<.i.i>—f .in’f río’J #» ‘exp(-0.F k1 Como la función depende de 4~ a través de 4, para dos ángulos torsionales. y>, y wz pertenecientes al mismo intervalo la relación de probabilidades será: p(~> > 5,9 1~toru (y>~ >> 1 Ñy>2) trcex? (—~V (y»)) (4.1.45) ~dy>,) = Este resultado permite establecer otra consecuencia de la discretización, y es que se considera que la influencia de las fuerzas intermolecuiares sobre la distribución conforTeacional es la misma para todos los ángulos pertenecientes a un miseo intervalo. De este modo la expresión máa adecuada para calcular la energía intramolecular es; t=1 tora ~ter.) ~y> (4.1.46) dy> 119 Coso ya se cosentó anteriormente desde ci punto rís vista práctico resulta sAs directo el cálculo de propiedades inureoieculares (energía intermolecular. presión) utilizando la ecuación (4.1.39) en lugar de (4.1.38), sin embargo mediante la expresión (4.1.45) se pueden relacionar ambas vías de promediado. Tal refinamiento no tiene importancia en eí contexto de este trabajo, ya que el número de confórmeros discretos empleados en el cálculo no afecta al tiempo de cálculo, pero puede resultar adecuado en otro tipo de procedimientos de simulación que aproximan situaciones continuas mediante la utilización de confortaciones moleculares discretas. Todos los desarrollos anteriores se simplifican en el caso particular en el que el ángulo de prueba 4. se genera con una cierta función densidad de probabilidad P,(4) sobre todo el intervalo de deflniclón~ esta situación puede ser descrite como: exp 1—$V,44)) = (4.1.47) ir AL> • 1/2 (4.1. 48) Oc 4 0+A4 40tA4 f PÁ*> d4 • ¡‘44) d4 =5 P,(4> d4 a 1 (4.1.49) Oc 120 4. 5oLScULs flU!BtES Teniendo en cuenta la condición de normalización utilizada para P,(4) y la periodicidad de 4 se cumplirá: A4—(41~t4) = ñO’. A4 = 2. (4.1,50) Esta elección elimina, por tanto, el problema del cálculo de los valores At y A4. Este planteamiento se ha utilizado en algunas de las simulacIones llevadas a cabo en este trabajo (secciones 5.2, 5.3 y 5.4). En la sección 4.2 aparece la discusión de la misma idea aplicada en una situación ligerasente diferente. 121 E ~~1 >1 FIGURA 4.1.1 POtENciAL TORSJ~AL PaR EL Ñ—BLiTANO 122 4. UnLEcuLAs rLrximl.ts 0 o o rICURA 4.1-2 DiSCRCTIZADGR4 DEL ANGULO DE IORSION 722 ¡ ¡ 020 1 2< 0.40 - 0. >0 a-> 000 OCUPA A QT~2SV6 II) 0 60 120 180 240 40 / grados <3 CISnIBUCION DE 4NGULDS 3CRSJON.LES Pk~ EL N—BUTAÑO ON DE LAS CONFOR<.4ACICNES a-rS Y IRANí SOORE LOS EJES kAoLEC uLARES 130 4. SoLicitas FtEXIM.55 e. ¿ 4> o A ‘9 e r e o b ¿3 ‘e e ~CAUCHE + .4.44 CAUCHE — 7ICURA 4.1.5 REP~ESENtACION DE LOS CONFORMEROS O Y O SOBRE LOS LIES MOLECUlARES 2.CCD~C GAUCHE •...~. GAUCHE + FIGURA t~6 cAMBIO DE LA ORIEÑTACIoN o>OtECULA~ EN UN MOYII*lEmO IFJTE~NO O — O b O ¶31 132 4 MolÉcuLAs tt.tXlBtfl 4.2 NCTOOOS DE REPTACI~ El término reptaclón” ~““ indica un cierto tipo de procedimientos utilizados en simulaciones por Monte Carlo para construir configuraciones de prueba. Su aplicación fundamental se da en el campo de simulación de cadenas lineales, bien sea en modelos de polímeras en red 3 o sistemas de cadena, en el continuo. Tienen especial importancia ya que, a menudo, en situaciones de gran empaquetamiento permiten el muestreo conformacional de moléculas flexibles. En situaciones de alta densidad otros 3 procedimientos pierden eficiencia de manera notable La aplicación del algoritao para una cadena lineal se puede explicar con los siguientes pasos: 1. Elección al azar de uno de los grupos extreaos de la cadena. 2. Eliminación del citado grupo. 3. Adición de un nuevo grupo en el otro extremo de la cadena. teniendo en cuenta las características particulares del sistema concreto3 (posibles ligaduras, sistemas en red, etc) 4. Aplicación de los teste adecuados (según cada caso) de aceptación o rechazo de la nueva configuración. En este trabajo so han realizado diversas simulaciones mediante métodos de Monte Carlo aplicando técnicas de reptación en 133 la generación de configuraciones de prueba. Los sistemas estudiados fueron moléculas de n—alcanos. Con el objeto de estudiar la influencia de la densidad en la estructura intramolecular se efectuaron simulaciones para sistemas con varias cadenas con condiciones periódicas de contorno (representando la fase líquida) así como simulaciones para una molécula aislada (situación de gas ideal). El modelo de potencial utilizado fue el propuesto por Ryckaert y Bellemans (ver sección 2.1); En este modelo cc :¿:cIJ:rzr distancias y ác’—.’”~ A. ~n!ace fIjos. sin esbargo los ángulos de torsión pueden variar de forma continua (ver sección 1.2). Al analizar el problema en coordenadas internas se demuestra’ (anexo D) que para el modelo de ligaduras flexibles el jacobiano de la transformación no depende de los ángulos de torsión. Existen diversas posibilidades para la construcción de configuraciones de prueba atendiendo al modelo utilizado en este trabajo; un posible esquema es el de generar el nuevo ángulo torsional con probabilidad uniforme en el intervalo [O.Zxl. sin embargo la existencia de la parte torsional del potencial intrasolecular condiciona de manera notable la efectividad de tal método, ya que la distribución de ángulos de torsión en el equilibrio esta considerablemente localizada alrededor de las posiciones de los miniaos de energía del potencial de torsión. Es por ello, que gran cantidad de configuraciones de prueba generadas mediante una elección uniforme del angulo torsional de prueba serio rechazadas por efecto del citado potencial de torsión. De 134 4. IO.ECILS rtnlmLn acuerdo con las ideas desarrolladas en la sección (3.1.)) resultaría interesante introducir de algún modo información adicional en el método de muestreo p.ra mejorar la eficiencia. lIna manera sencilla de hacerlo es recurrir a una elección del ángulo de torsión de prueba con probabilidad no uniforme, «(4). que tenga en cuenta los rasgos más importantes de la distribución de ángulos de torsión en el equilibrio. En este trabajo se utilizó: «(4> — c exp [ —0V~” (4)] . •s [0.2w] (4.2.1) Donde c es simplemente una constante de nornílzación. Esta elección de la función a es no simétrica (ver sección 3.1.1); la probabilidad de elegir una cierta configuración j (caracterizada por la inclusión de un nuevo ángulo torsional en uno de los extremos de la cadena y la eliminación de otro cierto ángulo 4~ en el otro extremo) como configuración de prueba a partir de la Configuración i es diferente a la probabilidad de elegir ¡ como configuración de prueba a partir de la configuración J. Este hecho, según se comentó en la sección 3.1. ha de tenerse en cuenta al construir el criterio de aceptación de movimientos. El criterio de reversibilIdad microscópica establece la relación: A 11 Pi «it 135 La relación entre probabilidades de las configuraciones j e puede escribirme como; — exp { -0[ lIL>(j)~tH(i) + lItera(j)—lI~~ í í} (4.2.3) — exp { —~ A(J’I -$ AU~~ } <4.2.4) La relación entre probabilidades a priorI, teniendo en cuenta la ecuación (4.2.1) y las características del modelo de potencial utilizado, resulta ser: U Cxp { ~0AlIt;r } (4.2.5) La relación de aceptaciones será por tanto: = exp { fllIt } (4.2.61 A diferencia de lo que sucede para el n—butano (en el modelo de Ayckaert y Belleaans) para los sistemas constituidos por moléculas de n—alcanos de mayor longitud de cadena la energía de interacción intrapolecular no viene exclusivamente determinada por los términos de lIt~ (sección 2.1). La estructura intramolecular en el caso del gas ideal depende también de los términos de 136 e. ,~tgctD.As rLtxlai.n interacción intratoleculares de tipo Lennard—Jones. Estos términos contribuyen a una mayor estabilidad de ángulos de rotación interna de tipo trans. este hecho se debe, entre otros motivos a que las Interacciones Li difIcultan considerabIes.nte la presencia de secuencias cZ G en las cadenas. Este efecto provoca, por otra parte, distintas distribuciones dé probabilidad para los ángulos de torsión en función de sus posiciones dentro de la cadena. Podrían plantearme cambios en las funciones ~ para tener en cuenta este efecto, sin embargo los cambios que aparecen en las poblaciones conforaecionales al tener en cuenta las interacciones intramoleculares de tipo van dcc Waals, en las condiciones de temperatura estudiadas en este trabajo, no son excesivos, por lo que pareció adecuada la elección expresada en la ecuación (4.2.1) para las funciones a. Por otra parte, tal elección simplifica de manera notable la interpretación del método: La influencia de las interacciones de torsión en la probabilidad de una cierta configuración se tiene en cuenta en el método de generación de configuraciones de prueba, mientras que los efectos de los términos Li se consideran según el procedimiento estándar en algoritmos Hetropolis=Nonte—Carlo. A menudo en simulación mediante técnicas de reptación se incluyen movimientos adicionales. Es frecuente encontrar en la literatura esquemas que combinan la reptacíón con movimientos moleculares globales de translación y e, rotación • la necesidad de estos covimientos es obvía para 137 modelos que trabajan en el continuo desde el punto de vista de las translaciones y orientaciones moleculares, pero que modelan los cambios conformacionales mediante la utilización de un conjunto 1,96discreto de conformaciones internas representativas El procedimiento utilizado para la simulación de n—alcanos en estado liquido consistió en la generación de una serie de ciclos de simulación. En cada ciclo se generó una configuración de reptación de prueba para cada molécula, esto se llevo a cabo de manera secuencial (sin elegir aleatoriamente la cadena a mover). 138 e. n.rc,a.ás rLn¡au 4.2.1 ANAUSIS LE EFICDEIA EN ~.ETWOSDE REPTACICt~ Los métodos de reptación empleado. para simular sistemas constituidos por cadenas moleculares difieren sustancialmente de otros tipos de técnicas aplicadas en la simulación de fluidos moleculares mediant, técnicas de Monte Carlo. En los experimentos de sísulación mediante reptación llevados a cabo en este trabajo no se incluyeron movimientos moleculares globales de translación o rotación. Aunque en principio el muestreo de tales grados de libertad se realiza implícitamente esa el proceso de reptación (cuando no se utiliza un conjunto de conformaciones discretas), los términos de energía intratolecular del potencial de interacción restringen de manera notable los posibles valores que las coordenadas internas (ángulos de torsión) pueden adoptar. Este hecho puede incidir en la eficiencia del muestreo dei espacio fásico configuracional. Es por ello necesario analizar detalladamente la calidad del muestreo que, para los grados de Libertad externos, lleva a cabo esta técnica de simulación. Por otra parte, las situaciones de alta densidad para líquidos constituidos por cadenas flexibles pueden dificultar de modo notable, coao ya se comeató en secciones precedente., el muestreo de los grados de libertad internos. LA generación de 139 configuraciones de prueba se basa en la modificación de los extremos de las cadenas, manteniéndose la estructura de las lonas sUs Internas de la molécula, por tanto el proceso de simulación ha de permitir la difusión hacia el interior de la cadena de los nuevos ángulos conformacionales generados en los extremos, a fin de que el muestreo de la estructura intraaolecular sea eficiente. Para llevar a cabo el análisis de la eficiencia en el muestreo de los distintos tipos de grados de libertad de estos sistemas se utilizaron diversas funciones de autocorrelación sobre las trayectorias generadas en las simulaciones. Estas funciones son en algunos casos análogas a las utilizadas en el cálculo de propiedades dinámicas en métodos de simulación por Dinámica Molécular. Lot procesos de translación solécular fueron analizados mediante el cálculo del desplazamiento cuadrático medio molecular. Z~. como función del número de ciclos de simulación; It(k) es el vector posición del centro de masas de la molecula a en el ciclo de simulación ‘Is’. A los efectos del cálculo el promediado anterior se realiza sobre todas las moléculas del sistema utilizando un conjunto de M’configuraciones generadas en la simulación separadas por un cierto número de 140 4. ZCI*45 FLUIW ciclos At. Li desplazamiento cuadrático medio se cálcula por tanto como: a 121 t(jht) — fr4¡zjy ~ fl.ÁlAt) — R.((iilAt) .3 = 0. 1 *4 j En esta ecuación ti es el número de moléculas, 14 el número de configuraelones utilizadas y ¿t el Intervalo de separación sisboliza la etiqueta del grupo que ocupa la primera posición de la cadena a en el ciclo de simulación ‘it De modo análogo a como se hizo para el desplazamiento cuadrático medio. el proceso de cálculo de 7, se expresa más adecuadamente mediante la ecuación; 1 42 4. IWLECW.A5 n2xImLzs II 51 F,(jtt) ¡4(14—ji 2 Zcos{~freiát>£.~íííAti 1 j(4.2.!.6)•fl I=l j0. 1.14—1 Esta definición cumple con los requisitos exigibles a una función de auto—correiacidn normalizada’: F%(0)=1. y F,(t)—0 cuando no existe correlación entre las etiquetas del extremo de una cierta cadena para configuraciones separadas por t ciclos de Simulación. La definición (4.2.1.5) puede parecer un tanto arbitraria. Teniendo en cuenta la periodicidad del proceso de reptación se puede plantear un isomorfismo entre el cambio de atiqueta del primer grupo de la cadena y el movimiento estocástico de una partícula en una red periódica monodimensional con N, nudos en la celdilla funda.ental. Las coordenadas de los nudos en una celdilla fundamental de longitud unidad serán: lc—i xk a . k • 1. ?4• (4.2.1.7) Supóngase que en un instante Inicial r 0 la partícula se encuentra localizada en la posición xÚ al cabo de un cierto tiempo t la probabilidad de que la partícula esté situada en una cierto nudo de la red podrá expresarme sediante una cierta función 143 periódica en el espacio P(Xkt)• esta probabilidad podrá desarroliarse en serie de Fourier. p(x~,t) a m 0 + (a100s (2iwx~> + b1sen(2in~)I (4.2.1.8) Donde los coeficientes del desarrollo dependen del tiempo. Suponiendo que la probabilidad de que la partícula se mueva hacia a la derecha es igual a la de que se mueva hacia la izquierda la tuncion p será par idada la periodicidad del mistemal. con lo que los coeficientes b1 serán nulos, la función p podrá desarrollarse como: 5e1 p(xk,t) — m,(t) cos (2XiXk) (4.2.1.9) 1 .0 El número de coeficientes del desarrollo, N~, será como máximo igual al número de puntos independientes: 14, • (N.+2)/2 . 14, • 2. 4. <4.2.1.10]1 ¡4, — (N,+1l/2 , 14, — 1. 3,.., Los coeficientes a 1 se obtendrían’ 1 como: ‘e’ R~i a 1(t) — p(x~.t) cos (2wix~) [~jcos2(2.ixk)] (4.2.1.11) lc.O k.O 144 4. Kotrcu,..s rInISLES Cuando la probabilidad p está uniformemente todos los puntos de la red: distribuida sobre <4.2.1.12) Se verfica: Nc —‘ Z Pj — 5—1 ‘=1 ¿4.2.1.17)j=O .1 14—1 Esta función tiene el valor 75=1 para jO y el valor F~aO cuando las direcciones de los vectores extremo—extremo de las moleculas en dos configuraciones separadas por un cierto número de ciclos no están correlacionadas. La elección de la función F~ puede justificarse de modo análogo a como se hizo para ~ En este caso Las funciones para desarrollar adecuadamente la probabilidad de un cierto ángulo entre loa vectores de una molécula para dram configuraciones separadas por un cierto número de ciclos de simulación son los 97.99 polinomios de Legendre del coseno del ángulo. Eligiéndose la función de mutocorrelación de orden 2 por la simetría de las cadenas. 147 5. RESUSADOS 5.1 ~-BUTAJ’JOLIQUIDO ( MODELO S»VPLIFIOADO En este trabajo se han llevado a cabo diversas simulaclones para fluidos constituidos por moléculas de n—butano utilizando el sodelo de potencial de Ryclcaert y Bellemasas (sección 2.1). Se realizaron simulacIones con los dos modelos de ligadura de coordenadas internas duras (ligaduras flexibles y ligaduras rixldas}, estudiándose la eficiencia del procedimiento de simulación desarrollado, ocasiona en la estructura En la tabla 5.1.1 simulados. Tabla 52.1 así como los efectos que intramolecular. se recogen los estados el medio denso termodinásicos a Estado T/lO kT/c d,lkg a~~) pc A 200.2 2.78 669.8 0.419 3 292.6 4.05 583.5 0.365 C 432.0 6.00 583.5 0,365 La elección de estos estados se realizó al objeto de comparar con resultados de Dinámica Moleoular~ 3 existentes en la biblIografía. Los puntos normales de fusión. T~, y ebullición, T,. 149 dci n—butano son ~<. 134.8 K y Tj 272.7 1< , considerando masas iguales para todos los grupos y teniendo en cuenta los valores de las ligaduras para el modelo de Ryclcaert y DeHesas (sección 2,1), adopta la forma: 2 3 5¿(. En la tabla 5.1.2 se presentan algunos detalles técnicos de las simulaciones llevadas a cabo para ambos modelos de ligaduras. En la tabla 5.1.3. se presentan las poblaciones conformacionales obtenidas para ambos modelos, en los tres estados estudiados. En la tabla 51.4 se comparan las fracciones de confórmeros trans del líquido (obtenidas mediante simulación) con las correspondientes al gas ideal para el mismo modelo de energía intrasolecular. La la tabla 5.1.5 se presentan los resultados para las propiedades teraodináaicas: 152 5. aU<AUS Ea la tabla 5.1.6 se comparan los resultados para el modelo con ligaduras rígidas con resultados bibliográficos de simulación por Dinámica Xolécular: Tabla 5.1.2 ¾ Br Cr A~ 3, 4 N.conflg/iO m 2.0 2.025 2.0 2.0 2.0 2.0 - %Trasl~ 30 30 30 30 30 30 t %Rotac. 30 25 30 30 25 30 t %Intra. 40 45 40 40 45 40 tAXmoa/nm 0.03 0.04 0.05 0.03 0.04 0.05 *A¿~ 1/grados 36.0 45.0 45.0 36.0 54.0 54.0 %Attca~ 43 46 42 43 46 41 § 40 47 52 41 41 46 XAíet A 57 57 56 54 55 55 ± ~1 %A(g n) 84.1 87.0 86.0 67.2 88.4 86.3 86.0 87.0 84.1 85.1 86.5 23.2 1 %Mt4g) 11.8 27.0 34.0 11.7 27.1 34.3 ¡1.2 25.7 31.9 10.5 26.0 32.0 35.2 52.3 56.3 36.0 51.9 56.4 N~(t4g)/10 4 16.8 55.1 65.6 17.4 57.1 66.6 Nt(g4t)/l0 a i&.s 55.1 65.6 17.4 57.1 66.6 N~(g±4gfl/l0 34 8.1 29.2 44.0 10.1 34.1 51.4 * Tanto por cierto de movimientos de cada tipo * Desplazaaientos máximos en un movimiento 4 Porcentajes de aceptación de los distintos tipos de movimiento 1 Porcentajes de aceptación para los tránsitos conformacionales * Número total de tránsitos entre confórmeros 153 Tabla 5.1.3 Tabla 5.1.4 Tabla 5.1.5 pA4Pa A, —39.40<4) 2.73(4> 11±4 8, —32.67(43 4.21(4) 1±4 -30.77(41 6.06(4) 95±3 A, —39. 37(41 2.95(5) 10±4 8, —32.70(4> 4.45(5> 2±3 —30.75<4) 6.39(4) 95±5 t Energías en unidades reducidas U aU/c. caSgS.6 J 154 5. RÉSUI.T5005 Tabla 5.1.6 Los resultados anteriores muestran la eficiencia del método de simulación desarrollado en este trabajo. El muestreo de la rotación interna es altamente eficaz, como se puede ver en el número de transiciones que tuvieron lugar entre los distintos confórseros, este número de tránsitos es considerablemente mayor que el obtenido en otras simulaciones por Monte Carlo utilizando el—ms procedimientos de tipo estándar . La precisión obtenida para las poblaciones cor,formacionales en el estado 8 es similar a la 7.obtenida por Brown y Clarice mediante simulación por Dinámica UieLar/C • 10±4 0.722(8) 1’DM —39.39 3.050 17.2 0.713 B, MC -32.70(4) 4.45(5) 2±3 0.613(61 —32.52 4.613 5.9 0.606(15> DMt —32.71(1) 4.47(3) 1.7(5) 0.609(7) C. XC —30.75(4) 6.39(4) 95±5 0.511(5) DM 12 (e/R)’) + c • RS 2 u (5.1 1.1) .0 ,n>a”%’ J Los sistemas con potencial intermolecular repulsivo se estudiaron utilizando el modelo de ligaduras flexibles. Los detalles técnicos de estas simulaciones son esencialmente idénticos a los remeDiados en la sección 5.1. En la tabla 5.1.1.1 se recogen los estados simulados. Tabla 5.1.1.1 Estado T/X ka/e d/(kg a’) pu 4” 200.2 2.78 583.5 0.365 r•1 9< 291.6 4.05 583.5 0.365 432.0 6.00 593.5 0.365 157 En la tabla 5.1.1.2 se presentan los resultados obtenidos para distintas propiedades, así como la comparación con resultados teórico. obtenidos a partir de la ainimizacibn de la energía libre de Helsholtz respecto a la utilizando para ello una ecuación de propuesta por Soublik’ 01. distribución conformacional. estado para cuerpos convexos Tabla 5.1.1.2 ACP 8L” 3.50(3> 4.77(4) 6.48(4> 2.B9(5) 4.38(3) 6.19(3> $p/p 13.9(1> 12.3(1) 10.7(1) (0p/Phroea.t 13.8 12.2 10.9 (X~,5,, >555 0.771 0.658 0.565 (Xt,...>i.5 0.729(61 0.611<6> 0.528(5] 1~ (K,.,,. ~ 0.687 0.603 0.515 e — 595.6 .1 mol’ t Referencia 80 Comparando los resultados obtenidos para los sistemas con potencial intermolecular completo (sección 5.1> con los sistemas con el potencial repulsivo VCA se observa un comportamiento diferente de las poblaciones conformacionales. £1 efecto de las fuerzas repulsiva, favorece el aumento de la fracción de 158 5. REsfiAs=s confórseros gauche en la fase líquida respecto a las poblaciones en fase gaseosa. Lete desplazamiento está parcialmente compensado por el efecto de las contribuciones atractivas al potencial intermolecular. Nte sismo efecto ha sido observado mediante un desarrollo teórico de tipo perturbativo para un modelo del ji—butano diferente al aquí considerado. Resultados obtenidos por 3. Brown’o’ mediante Dinámica Molecular para un modelo de ji—butano semejante al utilizado en este trabajo, en el que en lugar de utilizar ángulos de enlace rígidos se incluyó un potencial armónico, confirman el distinto papel que en el equilibrIo conforaaclonai juegan las contribuciones repulsiva y atractiva del potencial intermolecular. A diferencia de lo que sucede para la estructura Intermolecular de los fluidos atómicos densos 7 la estructura intramolecular de las moléculas flexibles no se describe adecuadamente en términos de las fuerzas repulsivas. Este hecho ha de ser tenido en cuenta en posibles desarrollos teóricos mediante teorías de perturbaciones. 159 160 5. REStLflXaO5 5.2 SM..LACION DE 1.2-DICLOPOETAN2 Se llevaron a cabo simulaciones de 1—2 dicloroetano utilizando el modelo propuesto por .)orgensen y colaboradores (véase sección 2.2). La técnica de simulación es análoga a la utilizada ea las simulaciones de ra—butano, aunque en este caso no se hizo uso de la aproximación de confórmeros discretos (ver sección 4.1), generándose los ángulos torsionales de prueba sobre el intervalo completo [D,2m1 de acuerdo con una cierta función de probabilidad P,(~>. No se incluyó el cálculo de las interacciones electróstaticas para el sistema periódico, en su lugar se procedió al truncamiento de las interacciones electrostáticas Según SC muestra en la sección 2.2, este hecho puede afectar a los resultados de propiedades que dependen de la estructura a grandes distancias (constante dieléctrica, presión, energía interna, etc). aunque tiene aenor influencia en propiedades que dependen principalmente del empaquetasiento de una molécula con sus moléculas vecinas 4. Las condiciones del sistema estudiado se encuentran recogidas en la tabla 5.2.1, 16! Tabla 5.2.1 TAO d/(lcg a~’1 H,,, R 5/nm 298.2 1246. 64 1.018 - Se realizaron tres simulaciones potenciales de referencia. u,(4); a exp [ —0U,4~) 1 3 114•) — Zat ~‘<~ =0 En la cada caso. utilizando distintos (5.2.1) (5.2.2) tabla 5.2.2 se recogen bm coeficientes utilizados en Tabla 5.2.2 1’ 1’ 1’ tCaso a1 a3 ter.A(U,aY ) 8.022 -12.067 1.393 21.451 B(U,.0 > 0.000 0.000 0.000 0.000 C 5.858 —17.473 0.000 24.421 • en Xi mol’ Las simulaciones se realizaron moviendo las distintas moléculas de manera cíclica. En cada ciclo de simulación se construyeron 3M configuracIones de prueba. (3 para cada molécula, una por variación de la posición del centro geoséti-ico molecular. 162 5. ~É5Ut?AI~S otra variando la orientación molecular y una tercera mediante la generación de un nuevo ángulo de torsión con la probabilidad descrita por las funciones P4qO). La forma de realizar los distintos movimientos fue análoga a la realizada para las simulaciones recogida. en las secciones 5.1 y 5.1.1. Ea la tabla 5.2.3 se muestran los detalles técnicos más interesantes de estas simulaciones. Tabla 5.2.2 N.ciclos(eq)1 N.ciclos(proa)* AX~ 54 A~,51 2.~l0~ 0.04 36 ¶ Número de ciclos de equilibrado * Número de ciclos de promediado 1 Demplmrmmientos de prueba máximo ( translación ). en na 1 Angulo de giro máximo ( en el muestreo de orientaciones, eai grados En la tabla 5.2.4 se incluyen los resultados de aceptación de los distintos movlsientos, así como el número de tránsitos conformacjons les. Ea la tabla 5.2.3. se presentan las poblaciones conforsaclonales obtenidas en cada simulación, así como las correspondientes fracciones para el gas ideal. Ea las figuras 5.2.1 y 5.2.2 se presentan respectivamente las las preb.bilidadem P,(•) asociadas con cada una de las 163 simulaciones y las las distribucion conformacionales obtenidas para el liquido. En la tabla 5.2.6 se presentan los resultados para las propiedades termodinámicas: Tabla 5.2.4 A 2 C ttres 23.5 23.4 23.4 ~ret 32.6 32.4 32.5 24 nt,. 22.0 11.4 23.6 • 1 XA(g 4g> 62.0 33.0 70.7 65.4 31.2 66.4 %A(g4t >1 + ZA(g.g 3.6 CA 0-0 1.3 1.2 1.2 1.2 0.0 0.0 Njt,g7 Mt(g..t) t N~(g±4gflt 3570 3578 33 3565 3419 1561 3417 69 82 § Porcentajes de aceptación ¶ Porcentajes de aceptación 0.117 164 5. ICSuLTAWS Tabla 5.2.6 Para el caso del 1.2 dicloroetano la eficiencia obtenida en la descripción del equilibrio conforaacional mediante el algoritmo de simulación descrito en la sección 4.1 es considerablemente menor que la obtenida en los cálculos llevados a cabo para el n—butano. Aunque existe una ligera mejora en el número de tránsitos conformacionales tránsitos respecto a los resultados 3., obtenidos por Jorgensen mediante técnicas de tipo Monte Carlo estándar, la variación nra es demasiado importante. La calidad del muestreo conforeactonal es prácticamente insensible a la función P,(W~ utilizada; no obteniéndose mayor precisión en los resultado, al utilizar una función de probabilidad a priori P,(#) próxima a la función de distribución para el ángulo conformacional. S(w>. en el liquido 4.82(16> 7.2(4) 4 4.0 2.0 Ii Y 1 1 j ~1 ‘\ ji It ¡ .1 2,’’ 1 0.0 0.00 0.25 0.50 075 ‘.00 rICURA i.z.~ PROBABWIDAD DE cLfcc,ON A 0R10R,, ?.(~). DE UN ÁNGULO DE ¶OPSL0N ~aRA uaS 4 3.0 2.0 -o 0.0 O - DO rICURA 5.2.2 DISTRÉSUOIONES SIM LI tAO lO NE5 DE ÁNGULOS DE TORSION EN EL UOULOO EN LAS TPEO 5. NESOLTADOS 40/2 ir 167 165 5. RÉstTaiaos 5.3 SINtLACION E N-ALCANDS En este trabajo se llevaron a cabo simulaciones de ra—alcanos en estado liquido utilizando técnicas de reptación. La formulación del método de simulación se recoge en la sección 4.2. Los sistemas estudiados fueron: ji—pentano. n—he,cano. n—heptano, n—octano. ra—nonanra y ji—decano. Las simulaciones se realizaron en el colectivo canónico (MVT), utilizando el potencial de Ryclcaert y Bellemana (sección 2.1>. En todos los casos la temperatura de los sistemas fue 298.1 K. Se utilizaron las ¡00 densidades ezperimentales correspondientes a la citada temperatura y presión de 1 atmósfera. En la tabla 5.3.1 se recogen las densidades y otros parámetros de la simulación: 169 TABLA 5.2.1 Se utiiizsz-on condiciones de contorno periódicas, con geometria cúbica para la caja de simeotación. Las coofiguraciones iniciales se construyeron situando los primeros grupos de cada cadena sobre una red de tipo cúbico consistente con el núsero de moléculas de cada sistema (cúbica simple para N=64, centrada en el las caras para N~32 y centrada en el interior para N—54L La posición del segundo grupo se fijó mediante la generación de un punto aleatorio con probabilidad uniforme sobre la superficie de una esfera de radio b (distancia de enlace) centrada en el primer grupo. De modo análogo se fijó la posición dci tercer grupo. mediante la generación de un punto al azar sobre la circunferencia definida por las posiciones de los dos priacros grupos, la distancia y el ángulo de enlace ib y e). Dadas las restricciones L/r Lime C 5Ha2 64 220 626 1.58 5.873 2.20 C611,4 54 224 659 1.67 5.789 2~7 C71$~ 54 172 624 1.74 &011 2.26 C511,5 32 256 703 1.79 5.230 2.05 %H20 32 282 712 1.83 5.299 2.12 CaoHaí 32 220 730 1.87 5.552 2.12 N~Núaero de moléculas, N,=Nú¡sero de centros de Interacción dDensldad, p,r’flenaldad de rentroslreducidal. Lbongitud de la caja de Simulación 170 5. RLSULTAOOS geósietricas del modelo de potencial, el resto de los grupos se aDiadieron sucesivasente generando ángulos de torsión • con probabilidad: exp (= 0ytO~’ (4)? d4 pt•)d• a (5. 3.1) 21! Las configuraciones así generadas resultaron ser de alta energía debido a los solapamientos de pares de grupos con interacciones Lennard—Jranes. Los sistemas se equilibraron, entonces, mediante movimientos de reptación según el algoritmo expuesto en la sección 4.2. En la tabla 5.3.2 se recoge el número de ciclos de equilibrado ( M~~) para cada sistenia. Un ciclo equivale en este caso a N intentos de reptación (uno para cada cadena>. TABLA 5.3.2 Fi número de cielos necesarios para la fase de equilibrado se detereinó analizando la estabilización de los valores de energía intermolecular, energía intramolecular. presión, poblaciones 171 confrarmaciona)es. desplazamiento cuadrático medio y funciones de autocorrelación de reptación y de orientación (sección 4,2.1). Tras descartar las configuraciones generadas en la fase de equilibrado, me prosiguió el proceso de simu)ación generando una nueva serie de configuraciones sobre tas que se evaluaron distintas propiedades: energias. presión, estructura ínter e intrasolecular. Por otra parte se grabaron en cinta nagnótica configuraciones cada cierto número de ciclos para utilizarlas posteriormente en el cálculo de errores y en la evaluación de las funciones dinámicas utilizadas para analizar la eficiencia del algoritmo. En la tabla 5.1.2 se muestra el número de ciciioa pr generados para cada sistema en la fase de promediado, >4 . y ex número de configuraciones equivalente TABLA 5.2.2 Sistema ~sH¡~ C¿H,• C,}l,a C 5li15 09>120 010822 >4Pr/íoa 0.112 0)98 0.400 0.528 1.280 1.500 N Nconí/lO 7.169 10.69= 21.60 16.896 40.96 48.00 Al obAeto de estudiar el efecto del estado liquido en la estructura intramolecular se llevaron a cabo simulaciones para sistemas con una sola molécula (representando la situación de gas ideal) a la misma temperatura y con idéntico modelo de potencial. En este caso no son precisas las condiciones ‘periódicas de 172 5. Rt5ttTAOO5 contorno, siendo el número de ciclos de simulación equivalente al número de configuraciones generadas. El procedimiento para generar la confIguración inicial fue análogo al utilizado para las moléculas en la fase líquida. Se realizaron para todos loa 6 6 sistemas 2.510 ciclos de equilibrado y 4.520 cIclos de promediado. En la tabla 5.2.4 se comparan los porcentajes de aceptación de movimientos obtenidos en las simulaciones en fase liquida y gas ideal /.A.~ y %Aaa! respectivamente. TABLA 5.3.4 Sisiemaa ges ZA 091a2 C 6H16 C~H16 C5fl~ 09Hac C,0H22 90.9 88.8 87.2 87.0 8~6 86.2 ZA 9.1 5.2 8.0 1.8 7.9 7.3 En la tabla 5,2,5 se comparan la las fsses líquIda y gas Ideal (x~ ~‘ fracción ‘a’ Xt de ángulos trajis en TABLA 5.3.5 e.’Xt U5 1! 12 C¿B¡4 ~ 09E18 CoTila CíoHzz 0.715 0.714 0.714 0.713 0.711 0.709 4 0.721(6) 0,725(22) 0.738(7> 0.726(20) 0.744(20> 0.754(12) )73 En la tabla 5.3.6 se desgiosa la distribución de probabilidad de los ángulos torsionales según la posición de cada ángulo en la cadena. Las barras de error para los resultados en estado liquido se estimaron en torno a 0.02. TABLA 5.2.6 Angulo(s) 24 e..24 C 58,2 •,. •~ 0.727 0.715 0.715 ~ 4~, ~, 0.699 0.782 0.692 0.686 0.771 %H~~ ,. 42. s 0.711 0.772 0.702 0,768 0.684 0.745 CaHas 4,, 42. 4~ 44 0.694 0768 0.246 0.705 0.771 0.684 0.742 0.712 09H35 •fl 4¡~, 4~. 46 4~ 4, 0.599 0.769 0.759 0,717 0.769 0.752 0.682 0.742 0.709 CioBaz 4,, 02. 4,, 4~ 46 0s 0.722 0.781 0.765 0.765 0.714 0,782 0.751 0.681 0.740 0.707 0.704 En las tablas 5.2.7—5.3.10 se presentan conforsacionales obtenidas en fase liqulda ji—pentano. ra—hexano y n-heptano. las poblaciones para las cadenas de 174 s. ~a.Ta~ rAtA 5.3,7. PoblacIones conforsacionales: ni-pentano Tipo M.C.E 1 X(Ll4uido> ~C (Gas ideal! T Gt 4 51.0(9) 51.4(1) T T 1 46.6(91 45.9(1) 0±0± 2 2.4(3> 2.7(1) 0±a 2 0.0 0.0 Número de conforsaciones equivaientes TABLA 5.3.8. PoblacIones conforsacionales: n—taexano Tipo M.C.E % (Liquido) %(Gas ideai) 7 7 0± 4 35.4±2.0 36.6(1) 7 7 T 1 32.2±2.0 29.9(1) 7 0±T 2 28.0±1.6 28.2<1) 6±7 0~ 2 5.5±2.0 5.1(1) 0±7 6± 2 5.1±0.9 5.5<1) T 0±0± 4 3.5±0.8 4.2(1) 0±0±0± 2 0.3±0.1 0.3 175 TABLA 5.3.9. Poblaciones conformacionales: n—heptano L En la tabla 5.2.10 se muestran los resultados obtenidos en el cálculo de la contribución torsional a la energía intraaolecular por cada grado de libertad interno (el número de grados de libertad internos es, en el modelo utilizado. igual al número de ángulos trarsionales, n 5) TABLA 5.3.10 C6H,2 09H,¿ C7H,6 4H18 Ceibo C,o}hz 2 3 4 5 6 7 teja ter. U115 /(tn33.93(4) 3.92(7) 2.83(8) 3.85(8) 3.82(9) 3.77(11) e a 599 Kj sol’ Tipo M.C.E %(Liquido> %(Gas ideal> 7 T 7 0: 4 25.8±1.4 24.0(1) T T 6±7 4 25.2±1.4 26.0(1> T 7 7 T 1 22.0±i.3 19.3(1) T 0±7 o; 4 7.1±0.8 7.4(1> 7 0±7 0± 4 6.8±0.8 8.0(1) 0±7 1 0± 2 3.4±0.6 3.9(1) 0±7 7 0 2 2.4±0.6 4.1(1) 7 7 6±0± 4 2. 6±0.5 2.6(1> 7 0±6±7 2 1. 1±0.3 1.8(1) 0±7 0±0± 4 0.7±0.3 1.1(11 0±7 6~ 0$ 4 0. 6±0.2 0.9(1) 7 0±0±0± 4 0.2±0.1 0.5(1) 176 5. O.LTAS En la tabla 5.3.11 se muestran los resultados obtenidos para las contribuciones Li a la energía interna, así como la presión obtenida para los sistemas en estado liquido. La presión fue calculada utilizando el formalismo molécular (ver anexo 81. Los resultados de energías se dan por grupo. TABLA 5.3.11 >gfltta¡(N) LL,tgr 3 U... 1 ¡(He> p« ¡e Gas Liquido CsI4í, “0.157 —0.158(1) 8.70(2) 1.1(3> C4{,, —0.338 —0.339(2) —9.09(4) —2.3(4> 09H,6 —0.492 -0.487(3) —9.39(2) —3.5(3) C5H,5 “0.620 —0.611(5) 9.59(3) ~4.1(6) 4H20 -0.729 —0.706(6) -9.fl(S> —4.7(4) CaoH2z —0.821 —0.787(8) —9.98(5) —5.8(3) —a 3c=5991ao1 e/o’ a16.47>4Pm162.Satm A pmrtir de los resultados de energías internas de las fases líquida y gas ideal se puede estimar la entalpia de vaporización, mediante las ecuaciones: A>!, = MI + pAV (5.3.2> 177 AH, = T4~S — + p ( V 555—V1~4) (5. 3.3) Despreciando el volumen del liquido respecto del gas. y aplicando la ecuación de estado de los gases perfectos se obtiene: AH.. a ?15,, - 1Ii,~ + RT (5.3.41 Este resultado se puede comparar con la entalpia de vaporización estándar, que se puede obtener experimentalmente. Dicha comparación se muestra en la tabla 5.3.12. TABLA 5.3.12 En las figuras 5.2.l-S.3.lE se presentan los resultados obtenidos para las distintas funciones de distribución. En las figuras 5.3.19—5.3.21 se recogen las funciones definidas en la 09H,~ CeNa, CiNte CeHae 091120 C,0}122 5 6 ‘7 8 9 10 8.93 8.78 8.72 8.67 862 8.58 AH, ¡(en.)t 9.6 9.8 10.0 ¡0.! 10.2 10.4 e = 599 .1 mol” 1 Resallados experimentales ( Referencia 104 * Resultado obtenido a partir de la simulación, utilizando la ecuación (5.3.41 178 s. mnttra~ sección 4.2.1 para analizar la caiidad del muestreo mediante la técnica de reptación. A partir de los resultados obtenidos se comprobó la capacidad de las técnicas de reptación para realizar un m&iestreo adecuado del espacio de fases de sistemas constituidos por cadenas lineales de corta longitud. Sin embargo, la eficiencia del método se reduce considerablemente para los n—aicanos de mayor longitud. Ente efecto se observa también en la aplicación de estas técnicas para otros tipos de modelos~ ~. En el campo de la simulación de modelos de polimeros en redes 3 se suelen combinar los movimientos de reptación con otros tipos de movimientos localms que modifican las posiciones de ura número limitado de centros de una cierta cadena de polímero, lográndose de este modo un tela Jamiento mis rápIdo de la estructura conformacional en el interior de las cadenas. Recientemente Xuaar. Vacatello y Yoon” han propuesto un tipo de movimiento de los grupos Internos análogo para modelos de cadenas lineales de polímero en el continuo, sin embargo no se ha demostrado que el citado procedielento moestree el espacío de aZ,C,,loS fases correctamente. Siepeanra ha desarrollado una técnica de Simulación por Monte Carlo generalizable para la simulación de sistemas de cadenas en el continuo, el procedimlmnto está basado en el algoritmo del paseo aleatorio autoevitante (Seif—avoidlng rindo. waJk) desarrollado por Rosenblutb y Rosenbluth’06, Este procedimiento, aunque costoso compwtacionalvaente en la simulación 179 de sistemas constituidos por cadenas de corta longitud, parece ser una vía adecuada para la simulación de cadenas de mayor longitud en medios de alta densidad, situación en la que los métodos de reptación pierden gran parte de su eficiencia’ 5. Como se puede comprobar en las tablas de resultados el equilibrio conformacional en el estado liquido está desplazado hacia un incremento de la población de ángulos torsionales ttans. Esta tendencia se incrementa con la longitud de la cadena, como consecuencia de ello el estado liquido favorece conformaciones moleculares alargadas. Este resultado puede ser debido a efectos de empaquetamiento. La concordancia con los resultados experimentales de entalpia de vaporización y presión no es muy adecuada, este hecho se atribuye al modelo de potencial utilizado. y en particular a un valor excesivo de los parámetros e del potencial de Lennard—lones para alpinas de las interacciones intermoleeulares. En la sección 5.3.2 se analiza con más detalle esta cuestión. mo 60,0 5. RwtTa~ 50.0 - 40.0 .~.3o.0 - 20.0 - 10.0’ 0.0 0.00 0.50 R/a 3.00 rICURA 5.3.1 ruNcloN DE DIS’TRIBuCION DE DI5TMIC~S IUTR*~4OLECUURE5. S PMA EL 5—HEXANO LI0LJIno o—PENTANO H n—HEXANO AIAiJ 181 60.0 50,0 40.0 - a.30.0 Eh 20.0 ~0.0’ 0.0 OLIO ~IOUP4£32 5% PaR. EL N—É4EPTÁ!40 &20.0 Ji 20.0 10.0 0.50 4.00 >50 200 ~ 7 a n—HEPTAIJO 1 ESO 3.00 rICURA 5.3.4 5(R) PARA Ea. N—oCTANo 60.0 50.0 40.0 182 60.0 5- .TAOOS 1.00 1.50 2.00 R/a FICURA 5.3.5 S(~> PARA EL N—NONAJ%i0 1.00 .50 R/a 2.00 2.50 3.00 FiGURA 5.3.6 5(R) PRA EL N—OECM40 n—NONANO50.0 - 40.0 - ~20.0- Cf> 20.0 - 10.0 0.0 0.00 o.k 2 3.00 60.0 50.0 400 Eh 20.0 ~0.0 n—DECAJ4O i0.0 0.yO .50 i83 1.40 1.20 - 1.00 o 80 ~‘O6O 040’ 0.20 0.00 oso 2.00 2,bO FÉCURA 5.3.7 FUNCION DE DISTRIBUCIOI4 RADIAL PARA lAS cONWíBUcIONES IN1’cRi.~OLECULAPES. g~(R), CN EL 5—PENtANO LLOUIDO FICURA 5.5.0 g,.RA N—HD(A~O 0011100 ti—PENTANO 1 140 1.20 l.00 0.80 ~‘0.60 OtO 0.20 184 1.40 1.20’ .00 0” 0.60 0.40’ o 20 - 0.00 0.50 1.00 1.50 rICuRA 5.3.9 ~.,n(R>PSRA N—I’IEPrM.J0 UCUIDO 140 1.20 >00 0.80 00060 0.40 0.20 000 FiGURA 5.3.10 a~ro ~A 5—OCTANO tIOllIDO 5. U.TAWS ti — H EPTANO tas 1,40 1.20 1. DO o so 020.80 0.40 0,20 0,00 0.. rcURA 5.3.11 ~enl~) PARA N—NONANO LIQUIDO 1~40 1.20 1.00 .0.80 00050 0.40 0.20 0.00 o. FIGURA 5.3.~2 ~jS~ PASA 4—DIDANO UOUíZO R a 186 n—PENTANO oM íM i.Éo R /u FiGURA 5.3.13 flItCCN DE DISTRIBtJCION RADIAL TOTAL íWTRA E INTEPI.AOLEOUL.ARES) P.C N—PENTANO L~uIDO 2.00 1,75 1.50 ~-ioo o’ 0.75 0,50 0.25 0.00 o- FIGuRA 5.3.14 9(R) PARA N.44C~40 uOIJOO g 32 730.0 1.87 298.1 *n—decano(b.d) 32 91.25 0.23375 298.1 * Número de moleculas CsH,2/11 292(2) 292(5) 298(2) 289(4> 4.06(3) 3.99(20) 4.16(5) 3.90(23) 29212) 297(15) 299(4) 231(16) —4. 0094 4. 7051 ‘5, 0526 4.3646 • 1 - u.u’j~, ~.D00S 0.0034 0.0002 Ueíac,1c5 4 3.23(2) 3.22(6) 2.67(21 2.59(4) • 1 0137 0.110 0.146 0. 155 UTRA —0155(3) —0.159(9) —0.835(5) —0.958(6) 406(1) 4,05(13) 4.23(2) 4.46(5) po-3/c —1.3(1) 0.21(1) —5.5(2) 0.134(7) p/plcT —10(1) 1.52(6) ‘7.2(3) 1.4313) * Temperatura reducidas ( 7 a kl/e,) 1 Energías reducidas por grupo 1 ¶1 •ti/fln,e X~(t) > — < X~(0) > e )C~(t) > <1 4 2X~(0> X~(0) 3’— e X~(0) (5.3.11) 199 Donde X~(t) representa la fracción de ángulos trmns de la molécula i en el tiempo t. En las figuras 5.3.1.5—S.3.i.8 se presenta el resultado obtenido para esta función en las distintas simulaciones efectuadas. De modo análogo puede estudiarse la relajación de cada ángulo torsional de manera individual. La función de autocorrelación adecuada para ello es — < X¿ 1(0) X~(t1 > — < ¶4kG) 3’ Xht ) rl < )4h0) )4’tD) > — ~4¼oí0 (5.3.1.2) En este caso x4¾t>es una función que sólo puede adoptar dos valores: x41(t)=í, si el ángulo de torsión j de la solécula 1 en el tiempo t pertenece al intervalo <2s/3.4;/31 y xV (t)aO en otro caso. Los promedios se realizaron sobre todos los ángulos de torsión de las moléculas. De este modo se puede visualizar el tiempo de relajación promedio de las rotaciones internas consideradas individualnente. En las figuras 5.3.1.9-5.3.1.12 se representan estas funciones de autocorrelación, A partir de estas funciones de autocorrelación se calcularon tiempos de relajacIón de la fracción 4 (4H~) 2.5 -10’ 80 ±0.014 MC (C~oH32) 4.8 -10’ 94 ±0.012 OS <4o~~,z) 2.3 -10’ 84 ±0.009 O Numero de configuraciones generadas (MCI o NSIDAD) 21 tiempo ¡ Ps 203 0.80 0 L0 O y 0.70 0.65 0.60 U ‘ño 200 • ióo tiempo / ;s 500 FICURA 5.1.1.1 EvOLt cLON Df LA FflcCION OS ANOLLOS DE RojAcIoN IÑOERNA PARA LA SIMULAZON POR 0;NA&AIÓA 14OLEOULAR LE N—IEOANO UOU!DO 0.73 0.70 0 c o OES >6 0.60 0.55 N—DGCANC (8 .D) 1 rICURA 5.3,1.4 EVOlUcION DE LA rRAccoI4 DE ANOLJLOS DE ROTÁCaON INItRNA PASA LA SIMuIACION POR DINAMICA ..doLEOuLAR DE N—DEcANO (FASE DE BAJA OCNSIDAQ) N — DECANO 204 N—PENTANO 1 1 ~1 44 ~1 1 1500 50 100 t/ps FIGuRA 5,3.1.5 FUNOLON ~L N—PENTAJ0 1.00 0.75 - 0.50 ‘a u- 0.15 0.00 ‘0.25 DE AuToCORPEtACION DE TORSIONES INDMOuALES PARA 50 t/ps 160 1 50 FIGuRA 5.3.1.5 ruNCION DE AutOCORÑtlACION DE TORSIONES INOMDIJALCS PeRA Si. 51575W 05 N—PEN1fl40 EN CONDICIONES CE BAJA OEI.ISIDAQ 205 1.00 0.75 0.50 tfZ 0.25 0.00 5. R~ULTAOOS N—R~NTANO (B.D) 1.00 0.75 - 0.50 tr 0.25 0.00 - —0.25- rICURA 5.3.1.7 nINcCN DE 0.. N—DEcANO .00 0.75 5 0,50 0.25 0.00 —0.25 50 t/ps 1¿0 150 AuTOOORREUCION DE tORSIONES INDIVIDUALES PaRA ricuRA 5.3.1.8 FUNCION Df AISTOCORPEl.ACION Of TORSIONES INOMDtIALES PARA EL SISTa Of N—OEOaiJo f*. cONDIcIONES OC SA.» OCNSIOAO Y N—DECANO t/ps 206 1.00 5. W~flU.TAO05 0.75 0.50 0.25 0.00 —0.25 flGVAA 5.3.1 .C ruNcloN DE AuTOCo~tlJCION DE LA Fft&CCION DE AMGULOS DE TO#~SON TRA>É5 POR NOLEEIJIJA PARA NPDJTM4O uQuáflO 1.00 0.75 - 050 0.25 - 0.00 t/ps 160 O ricisá s.3A .10 FUm004 Dc AJfl0CO4~AtJU4 DE LA FRACCION DÉ DÉ TORSm IRASIS POR CUVIAA Pa fl-PDflMIO EN CONOfOCNE5 DE SAlA 0068W t/ps N—PENYANO (8.0) 2 207 1.00 0.75 Li. 0.50 O .25 0.00 —0.25 rICURA 5.3.1.11 rtJNCI0t~ DE AUToCORREUC¡0N DE LA FRACCLON DE M~4GUL0S DE ¶ORSION NANS POR iAOLECuLA PARA N—DECANO LIOUIDO 1.00 0.75 0.50 ~1< u- 0.25 0,00 - —0.25 1 50 t/ps do rICuRA 5.3.1.12 rUNCION DE AuTOCORRELACION DE LA rRACCION DE AHGULOS DE TOR9ON X~Ab.S POR RtAOLECutA PaRA N-DECANO E$ C01401C04t5 OE BAJA DEI4SIDW N—DECANO (9.0) 2014 Rt. arsuitAnOs 5.3.2 REFINAMIENTO DE PARANETROS DE POTENCIAL En este trabajo se analizó. mediante el foraslismo de repesado de configuraciones. la influencia de los parámatros de interacción Li del potencial propuesto por Ryckaert y Hellemans en las propiedades termodinámicas de líquidos constituidos por moléculas de n—alcanos. La razón fundamental de llevar a cabo tal análisis fue el hecho de encontrar considerables discordancias entre los resultados e%perimentaies y de simulación para propiedades como la presión y la energía intermolecular. Estas deficiencias aumentan con la longitud de la cadena de ra—alcano. de tal manera que los parámetros adecuados para el cálculo de las propiedades del o—butano dan lugar a maloa resultados para n—alcanos de mayor longitud. Diferentes investlgadores5IOl haza comprobado la imposibilidad de construir un modelo general para estos sistemas lalediante la utilización en la descripción de las interacciones intersloleculares e intramoleculares de tipo van der Waals de los mismos parámetros Lennard-Jones para los grupos metilo (—C1l~ 1 y metileno 1—Lila—) que constituyen cate tipo de moléculas. Era el caso de utilIzar los mismos parámetros Li para ambos tipos de unidades se ha visto que los parámetros óptimos varian en función de la longitud de la cadena. Se han propuesto diferentes 209 alternativas ante esta situación. El procedimiento etAs habitual para construir un potencial transferible para estos sistemas consiste en la utilización de diferentes valores de los parámetros o y o en función de la naturaleza quimba del grupo. Parece bastante razonable la distinción entre grupos metilo (—CH 0) y metileno (—ElY2—) y es en es esta línea donde se han desarrollado algunas de las mejoras Otras mejoras propuestas consiten en emplear diferentes ángulos de enlace o distancias de enlace en función de la longitud de la cadena de hidrocarburo~O?. En otras modificaciones se desplazan los centros de interacción desde tos nucleos de los átonos de 107 carbono a nuevas posiciones que tienen en cuenta la presencia de hidrógeno en los grupos metilo o metileno. En este trabajo se ha analizado la dependencia de las propiedades termodinámicas con respecto a los parámetros de potencial que definen las interacciones de van der ¿asís entre los distintos grupos. tratando de buscar un conjunto de parámetros que de cuenta de las propiedades en fase líquida de n—alcanos con distintas longitudes de cadena. Por otra parte es posible encontrar en la bibliografía diferentes funciones para la descripción del potencial de rotación internaS?, La dificultad en la obtención de la población a 05 conformacional etediante técnicas experisentales así coso la complejidad de realizar cálculos mecano—cuánticos relacionables directamente con un modelo de potencial que considera distancias y 210 5. eSSUtTáDOS angulos de enlace rígidos son las razones que llevan a esta indefinición. Sin embargo, la influencia de los potenciales de rotación interna en las propiedades tertodinámicas puede cor,slderarse como senos transcendental, ya que las distintas funciones propuestas se diferencian básicamente en las alturas de las barreras entre los distintos sinimos de energía (conformaclones & , g y tI, lo que tiene sin duda importancia desde el punto de vista dinámico, sin embargo la población de conformeros no suele variar notablemente entre los distintos modelospropuestos. El problem. estudiado puede, por tanto, plantearse en los siguientes teratitos: El modelo de energíapotencial tiene la forma: u • ~ U~ (5.2.2.1) itt ¡E + ¡E 15.2.2.2) e 1 k~’ 1*1 La forma fui~cional y los parámetrosempleadospara la energía de torsión son los del modelo de Ryckaert y Bellemana LTAOOS parámetros propuestos por Ryckaert y Bellemans dan lugar a predicciones bastante alejadas de los resultados experImentales para los valores de presión y entalpia de vaporización, estás han sido las propiedades que ha,, sido utilizadas en el proceso de optiwización. La entalpia de vaporización se expresa de ando aproximado (ver sección 5.3) como: AB,5lIg,g~¶Jjj~ + RT (5.3.2.8) Separando las contrIbucIones Intra e Intermoleculares puede escribirse: Alta < 9494<~94;~ra> + ~ — Int.r (5.3.2.91~Iíq Suponiendo que la diferencia entre energías intramoleculares entre las fases liqulda y gaseosano depende considerablemaentede los parámetros LI del potencial podrá escribirse: 09 letra Aiim AH, — ~T 5 (fllDtt~) ~ 1 (5.3.2.11) flq a ínter (11 ) a —AH, • lot,. (5.3.2.12) a ~(~9es~tlIQ )«~ 213 En esta ecuación el subíndice a representa un cierto conjunto de parámetros íí y {ij}~ podrá escribirse: tkr = ~ 4c,í[ ( Ú,i/R,il’=— ( o,1/R~l 6 ] (5.3.2. 14) 1~1Itt 214 5. Rt5tA.tAOOS Separando las contribuciones repulsiva y atractiva <~~r Y UCC U~e • (5.3.2.15) Ws 4 ~ o,, (5.3.2.161E podrá relacionarse la contribución U~ a la energía interaolecular de una cierta conf ignraciofi Ic’ calculada para un cierto conjunto de parámetros a con la obtenida para los parámetros «o mediante la ecuación: 12 12 IUrc(k)] [Uh(kfl ( É1, «~Q ~0 Oo (5.3.2. 18) «o De modo análogo se obtiene: • e 6 6 lUcc(k)l« ~Wc(kl] «o ( e,, e1,’ e0 o0) (5.3.2.19) Expresiones semejantes se obtendrían para el resto de contribuciones, de tal manera que ci cálculo de la energía de una 215 cierta configuración utilizando parámetros a se obtendría sediante adición de las distintas contribuciones en el estado de referencia, «o. multiplicadas por un cierto factor de escala: IuIatrocn = c 110’U Oc)] CnCfl (14, Oc)] + « «~ 4=)~ 0t~ % 01 ~ [U~,(k)] + j 12> a 6> a O ~¶ 201( C220lZ) fu~(k)l + ______ e0 oc o co 4) « De modo completamente análogo se puede proceder con la presión, separando las distintas contribuciones y escalando adecuadamente. Utilizando configuraciones generadasen las simulaciones de a—alcanosse llevaron a cabo diferentes cálculos de optimizaclón. El criterio para optimizar consistió en obtener el parámetro o los parámetros que minimizaban una cierta función objetivo, O (a), donde a representa el conjunto de parámetros a optimizar. La función Ola) utilizada fue: O(a) = ap E pj(a)—4 )~ + a, ~¡ A114(a)—A14 ]~ <5.3.2.21) a a Donde los suisatorlos se llevan a cabo sobre uno o varios omistes.., dependiendo de los casos, p4 y AI4 son los valores experl.entales. Los coeficientes a, y a, están relacionados con 216 5. esslJttAoOS las barras de error obtenidasmediante la simulación. “‘1. —1a, ~4—E ~ (9>) (5.3.2.23) 1—l •36• aj’ ~7~’ E SJ(’Ai{,>) (5.3.2.24) 1=~ Siendo N• 1• el número de sistemas considerados (6 en el presente caso). Los valores utilizados fueron a,—6.5 Ifa y man 21 4 (aol de -1 0 0grupos) Los valores experimentales MI1 y p utilizados se recogen en la tabla 5.3.2.1. Tabla 5.3.2.1 n-C,111~ n—C68,4 >rC7H34 n—CeB,s n-CeHzo n—C,0822 p 0fl4Pa 0.1012 0.1013 0.1013 0 1013 0.1013 0.1013 1’ o —, A11 5/(Uaol 1 4.85 4.85 4.87 4.88 4.89 4.89 Expresada por aol de grupos Se llevaron a cabo dos tipos de optíaizaciones: optíaizaciones para sIstemas individuales y optimizaciones globales. Incluyendo todos los sistemas. En ambos casos se 217 calcularon los valores de las propiedades utilizadas para el ajuste para todos los sistemas, Por otro lado se evaluaron el valor del parámetro de efIciencia para cada sistema y las funciones O individuales y globales. En primer lugar se abordó la optimización del valor del parámetro t de Lennard—jones. considerando e igual para los grupos metilo y metileno, y utilizando el valor de o de). potencial de Ryckaert y Bellemaxis. En este caso se llevaron a cabo optíaizaciones para cada sistema independiente (6 optíaizaciones diferentes). Los resultados obtenidos se reflejan en la tabla 3.3.2.2 219 a. v’flAflDos Tabla 5.3.2.2 Posteriormente se llevaron a cabo optimizaciones para cada sistema independlente, considerandodos parámetros a ajustar, e y o’. Coso en el codelo de Ryckaert y Bellemans se consideró: Ocx,ocn~’.C y Cc’4,t~’4 2*C. Los resultados se presentanea la tabla 5.3.2.3. Obviamente se obtienen valores muy buenos de las funciones objetivo individuales, ya que se ajustan, en cada caso, dos propiedades (presión y entalpia de vaporización) mediante la OPTIMIZACIOI¿E3 INDEP»¿DIflI1~ PARM4EflOS A OPTIMIZAR: e l~ESTRICCI0NES: t • — 53. 0jS03 — 0. 3923 n~ WIIflD bE 0P’TIM!ZACIONES • 6 sístna ca~ pAlPa 611 ~ o~ n—C 611,2 67.S 0.717 —4.5 4.85 0.6 n~CóI 4í, 63.0 0.359 —10.7 4.84 3.0 1’10’ n—C,W. 63.2 0.168 —21.6 4.26 II. 7.102 n—C,H,, 62.2 0.240 —25.4 4.87 16. 7.102 n—C,ff 20 61.6 0.092 —24.7 4.88 15. 9’10 n—C,0H23 60.2 0.068 —36.4 4.87 32. 210’ • Xélvin * Eficiencia del repesadosobre el misten optíaizado 4 LI! < ¡sol de grupos) O Función objetivo individual del sistema ajustado 1 Función objetivo global: Suma sobre todos los sistemas 219 introducción de dos parámetros. Tabla 5.3.2.3 Como se deduce de los resultados anteriores la utilización de los mismos parámetros Lennara-Jones para los grupos metilo y metileno no conduce a una buena concordancia de resultados experimentales y de simulación, los parámetros que ajustan las propiedadesde un cierto n—alcano en estado liquido fallan de modo obstensible al aplicarlos a otro ti—elcano. Se han propuesto~ OPTIMIZACIOHES INOEPDCI ENTES PARA>4ETROS A OPTIMIZAR: e. o HESTRICCI0N~: cma, mc 1, o = o1 = o’2 MUflO DE OPTIMIZACIONES = 6 SISTfl4A t1 0 o’/¡jjj p/MPa Al!) o n—C 5H,2 67.5 0.3934 0.89 0.1 4.85 0.01 310~ n—C6H,, 63.6 0.3951 0.77 0.3 4.85 0.01 líO 2 n—C,B, 6 60.1 0.3985 0.13 0.4 4.87 0.07 7.102 n’~CH,6 59.8 0.3975 0.24 0.8 4.88 0.07 7.102 ~—%>~~e 58.6 0.3982 0.32 —0.2 4.89 0.10 l’l0’ n—C,01122 56.6 0.4000 0.14 0.2 4.89 0.01 t Kelvin * Eficiencia del repesadosobre el sistema optimizado Idi ( aol de grupos) * Función objetivo individual del sistema ajustado 1 Función objetivo global: Suma sobre todos los sistemas 220 5. R15t5.7A005 modelos de potencial para estos sistemas en los cuales se establece une distinción de parámetros en función de la naturaleza química del grupo, siguiendo esta idea se realizaron ajustes globales de parámetrosutilizando todos los sistemas. El primero de estos ajustes considera el .ais.o diámetro para los grupos metilo y metileno, pero diferentes valores de e. Los resultados obtenidos se recogenen la tabla 5.3.2.4. 221 Tabla 5.3.2.4 Finalmente se llevó a cabo eí ajuste de cuatro parámetros, considerando o y e diferentes para los grupos metilo y metileno. Los resultados obtenidos se muestran en la tabla 5.3.2.5. OPTIHIZACION GLOBAL PARM4EROS A OPTIMIZAR: e,, e 2, o’ RESTRICCl0N~: o = o’, o’2 PARMCTROS OPTIMOS r,ik 5 Cc5~ /k 89.0K = ccc 1k 52.1K o’ = «CH3 o’052 0.3950 ni. 1 O, — 37. * SISTD4A pil4Pa AV. 4 0.06 19. 4.85 8.8 n—C’ b 5 RZ5¶LtAI~S lm e u e O “5 ‘e — — — — — — •2 — ‘0 ,e. EFICIENCIA: cc,,./k=SOK. cc,s/kSOK FIGURA 5.3.2.3 EFicIENCIA MEDLA DEL MUESTPEO MEDIANtE REPESAOO o’cc3.SOA, ~c 4s4.OOA 18 80 84 U 22 90 1 <»1 “o EFICIENCLk 20 8* 10 72 a- ~3-B~ 1 <3 U’” 5,5 3-” 3S 7* 72 e 54 SI .4 —~ •~ ‘~ ~CE4,•0 •0 40 U *5 S ,oe FiGURA 5.3.2,4 rncIcNct MEDIA Da MUESTREO MEOWTE REPESADO’ 226 Log 0~: Ecs~Vkm90K, , 4/kn5OK 3,5 3-22 SS 3-tí 5.03 5,5 S.4~ — 4S1 4.05 4S ‘o. ‘o, 4,” 1 Ql 3,5 311 FiGURA 532.5 LOG FIGURA 55.2.6 LOG 0, Sal 3-42 Se 5,5 5. RWA,TAVOS Loq 0~: ~ 3.2DA. ~.r4.O0A lO U U 35 18 lO U 80 48 80 80 80 ‘3 ‘4 1 £3: •0 SS tCH, ‘e ~ ‘U 227 228 5. RtSttTAflO5 5.4 N’BLJTAHZ UOUI~ < MODELO DETALLADO 1 Se han llevado a cabo simulaciones para el n—bvtano mediante métodos de Monte Carlo espleando .1 modelo propuesto por Vilo y Yip (sección 2.3). En la tabla 5.4.1 ser recogen las condiciones de los sistemas simulados. Tabla 5. 4. 1 SISTD4A T/X d/m. <«IP. r,, son los valores inicial y perturbado del ángulo de torsión definido por los carbonos, Análogamenteb15, y b.~ representanlas distanciasde enlace. Para elegir la función VÁr1 ) se realizó en primer lugar una sisulación relativamentecorta para el misten ideal en la que V~fy~) se consideró constante, a partir de los resultados de la función de probabilidad S(11 1 (sección 1.1) obtenida y teniendo en cuenta la temperatura se obtuvo la función VÁi1í utilizada en el resto de las simulaciones: (<7,i 5$t lo (5.4.8) So ~n 1 Donde < representa el valor sáximo de la función de probabilidad obtenida. 5o~ La función se representó como un desarrollo en serle sobre el cosenodel ángulo: 235 2 Vt(y) Z a~ tao <5.4.9) En la tabla 5.4.3 se recogen los coeficientes de este desarrollo, así como las variaciones máximas permitidas en el resto de los sovialentos: Tabla 5.4.3 Para generar la configuración de partida en la simulación de la fase líquida se situaron las moléculas sobre una red cúbica, con orientaciones al azar, como distancias de enlace iniciales se tomaron los valores de equilibrio de las contribucionesde tensión a la energía intraseolecuiar . Los Translación molecular: AX — Ay — Az a 0.40A htaci6n molécular: — 45 grados Movimientos intraaoleculares: Ab — o.oíÁ Aca 0.01 A 7 — 0.02 rad 5 1.15 grados Coordenaday,. coeficientes del desarrollo de 1< 0 1 2 3 —la8/(Kj aol ) 8.369 —14.226 0.000 22.594 236 5. RWL¶A5 demás ángulos y fueron escogidos de maner, que los ár>guios de torsión definidos obre lineas de enice valiesen 60, 190 o 300 grados. En la simulación del gas Se procedih de modo análogo. En 1* tabla 5.4.4 se recogen algunos detalles técnicos de estas simulaciones: Tabla 5,4.4 GAS LIQUIDO Fase de equilibrado: t 5 Número de ciclos 10 510 Número de configuraciones 10~ 0.9610~ Fasede promediado: t .7 5Número de cielos 10 1.6’1O Núaero de configuraciones 10’ 3.m2~loe II Aceptación: Translaciónmolecular — 22. 5 Rotaciónmolecular — 21,8 Intramolecular 45.5 20. 9 t Un ciclo equivale en estecaso a haber generado una configuración de prueba por cada tipo de movimiento para cada una de las mo]átul.s del sistema. 1 1 ciclo < > 3 al configuraciones, en el liquido 1 ciclo < > 1 configuración, en el gas ideal, donde los moviaiemtosde translacióny rotación no son necesarios 1 237 En la tabla 5.4,5 se recogen los resultados obtenidos para las distinta, contribucionesa la energíapotencial (véasesección 2.1). Tabla 5.4.5 GAS LIQUIDO “u U±nt.r/(kJaol’) — —24.6(1) vn, -t U 1,~~10/(kJ sol ) 21.5(1) 21.5(5) tse —, U lOa •ol 1 12.1(1) 11.1±1.0 tj 4~”/(IcJ aol’) 15.2(1) 14. 3*1.1 tOra —l U ¡(Id mol ) 0.71(1) 0.49(1) (Id soY1) 9.1(1) 7.3(9) líq 88~ 1•(a) 1.555<2> 0.031 0.031 1.563<1> 1.556<2) 0.032 0.032 9 — 3 ángulos en grados En la tabla 5.4.7 se presentanlos resultados las distribuciones de ángulos torsionales. Tabla 5.4.7 239 obtenidos para En la figura 5.4.2 se comparan los resultados de simulación para la función de distribución radial global. g 8(r), construida pesandolas distintas contribuciones parciales de acuerdo con las longitudes de dispersión de neutrones de los núcleos para el 9sistema n—butano deuterado, C.~~0 con resultados experimentales para el mismo estadotersodináslco. Como se puede observar la concordancia es siasplesente cualitativa, la función obtenida mediante simulación presenta picos más altos que la experimental. Esta discordancia puede ser debida al proceso de inversión de los datos escperiaentales para obtener la información en el espacio real. al modelo de potencial empleado o a la utilización de una descripción puramente clásica en la obtención de los resultados mediante simulación. 2>40 a. TA03 5.0 4.0 3.0 o, 2.0 ‘o 0.0 FIGURA 5.4.2 COMPMACION X ¿A PUNc¿on DE D~WaIC>ON&AD& DE NOSTsiU4tS . En ciertas ocasiones existirán contribuciones a la energía potencial que vendránexpresadascomo funciones de las coordenadas internas no ligadas. En estos casos podrán aparecer a..bigtiedades al calcular las fuerzas originadas por tales términos potenciales sobre los centros de interaccion. Tal problema no apareceríaen el caso da trabajar con coordenadasgeneralizadasadecuadas, sin embargo resulta notablemente más sencillo en simulaciones por dinámica molecular operar con coordenadascartesianas e incluir fuerzas adicionales que mantengan las coordenadasinternas ligadas Co! sus posiciones de equilibrio. Los términos de potencial intrasiolecular pueden venir dados en funcion de coordenadas internas tales que al expresarlas en coordenadas cartesianas incluyen coordenadas internas ligadas. En tal caso pueden utilizarse dos caminos para proceder a la derivación de las tuerzas. La primera posibilidad consiste en tratar tales coordenadascoso constantesal calcular las fuerzas intramoicculares. alternativamentepuede no tenerse en cuenta la Z45 A—os consideracion anterior. Los resultados difieren según el procedimiento utilizado. Estas diferencias no conllevan un comportamiento dinámico diferente, puesto que al aplicar el correspondiente método de ajuste de posiciones para mantener las ligaduras del sistema las fuerzas resultantes coineiden. Para apreciar de manera más clara lo expuesto anteriormente se incluyen en este anexo dos ejemplos. 1. POTENCIAL DE FLEXION Sea un potencial de flexión dado por una cierta función U (cose). Siendo e es el angulo entre dos enlaces sucesivos de unae cadena lineal. Si A A y fi son las posicionesde tres centros 1 2 5 sucesivos de la cadena. cos O vendra definido por: fi .R cosi o 2! 23 ~ 2 (A.13)acos• L b sine 0 2 aco0 1. b 2 sen2O 0 au 7< • r __________ ¡ ChIS) Bcos$ L b 2sen2e .> o Bu (cos•l ~ft~ 3 1 0 ___ ,~ acos• b 2senO o BU (B. 10)—E-. 2S2 Transformando la ecuaciónanterior se obtiene: < —L 1 BU) > LB.]]) a (~VJ3N OL 3L e — • < -L [.4j—)~,,> iB. 12) En el primer término del segundo miembro de la ecuación anterior aparece la contribución cinética, el segundo miembro incluye el efecto de las fuerzas lntermolec’alares. En el desarrollo previo se ha definido la presión medianteun proceso de escalado de posiciones atómicas, para sistemas en los que existen ligaduras que mantienen constantes ciertas distancias intramoleculares tal procedimiento no es el más adecuado~ resultando sás conveniente un proceso de escalado de posiciones moleculares en el cual se escalan las posiciones de los centros de masas soleculares (o cualquier otro punto de referencia molecularí. En este caso es conveniente clasificar las coordenadas del s~atema en des categerias: por una parte apareceran 3M coordenadas (siendo 4 el número de moléculas del sistema) que indicarán las coordenadascartesianas de los puntos de referencia de las moléculas (simbolizadas como q~’<) y por otra parte el resto de coordenadas del sistema (rotación y coordenadas moleculares Internas), que serán denotadasmediante un vector de dimension SU.353M—3M: q, 2 . Operando de modo análogo a como se hizo en el caso 253 atómico pero realizando el canblo de variable (E. 12) (2.13) <8. ±4) 3±4 2±4 — L «t 2±4—3M 3M.3< Se obtiene una ecuación análoga a la (BtU: II 1 (ay 2V <—L (~L ia~ q t rl Oc nodo semejante la presión, p. pt±ededefinirse a través de 13la función de partición canónica ~V9) e. T 1. 0 - L I dx, dxínl { dp, ~. dp 31 exy, Realizando el siguiente cambio de variable: — L a, Pi i,/ L (8. ±5) (8.18) 254 se obtiene: 1 1 0 0 0= ~“4d.”dasn j..’~ dy±”’’dfl 1 exp ~j—0M(L.sv)j (9.19) Siguiendo la misma ruta que en los casosanteriores: = —~ < L 1 afl(L.a.,flpv ~ BL > (9.20> ‘rs, y Según la transformación de coordenadas previa la energía cinética podra escribirme coso: 1±4 2 <2.21) 1~1 Desarrollando la ecuación (2.20> se obtiene: «.7 2 1 (BU ~L taL 1> <2.23) Por otra parte, desarrollando la derivada de ti respecto de L se obtiene: 253 AMEMOS a. JLL 8x¿J BLJ« (2.241BL ¶5~ Mediante las ecuaciones de cambio de coordenadas (ecuaciones <2.17) y (B.l8l)~ la expresión anterior se transforma en: 3M Bx <2.251(BLJL( Y la ecuación de presión podrá reescribirse, utilizando notaclón vectorial, como: m pV x ~ (B. 26) —l O en función de las fuerzas que actúan sobre cada centro de interacción: ±4 pv = < ~ + < E, > (3.22) 1—l Aparentemente las ecuacIones (3.11) y (3.27) son totalmente equivalentes, dada la relación existente entre energía cinética y temperatura, sin embargo esta última es más adecuadaal trabajar con dinámica molecular. =56 APETOcE C TRANSFORMACI~4 DE LA NTEGRAL DE CfrFIGtEACD~ La función de partición canónica clásica expresada en coordenadascartesianaspresentala forma: Q e J.. 4 exp [ -~ ~ dp 2’t es la energía total del sistema. En sistemas conservativos, en los que no existe intercambio de energía con el exterior puede expresarse 1< como: < 3±431% = y + R< es la energía cinética, que sólo depende de los momentos. La energía cinética puede expresatsecomo: 2 3% 2 3ct 4 es una matriz diagonal en la que aparecen las masas de asociadascon cada grado de libertad, y el símbolo representa la transposición de un vector o una matriz. Para no complicar la notación se han suprimido los superindices que indicaban las dimensiones de los vectores. Utilizando la ecuación — CC. 29) La obtención de Z~(q) puede lograrse mediante un método diferente, en lugar de llevar a cabo el desarrollo sobre variables generalizadas y proceder posteriormente a la integración sobre los 263 AMEMOS momentos generalizados, podría procederse en primer lugar a la integración sobre los momentos cartesianos, y a continuación realizar el necesario cambio de coordenadas de posición. Utilizando coordenadas cartesianas el valor medio de una cierta propiedad Mx) viene dado por J dr exp[.-st¡ CC. 30) f dx exp[-sii~ El cambio de variables z a las nuevas variables q origina una trmnsforn±acionde las integrales para dar lugar a: f dq exP[—PtLC~)j iCq)j (C. 31) J dq eXP[—$tI(q)] J(ql¡ Siendo J(q)l el jacohiano de la transformación (ecuación (Cg)). Obviamente habrá de verificarse la equivalencia de las ecuaciones (C.29> y (C.311. Recordando la definición del teMor métrico G en funcion de las matriz .7 y 71 (ecuación (Cl?)) podrá escribIrme: 1/2 — ng’” jJ¡<3 (C.32) 264 Como 71 no depende de las coordenadas obviamente se verifica la equivalencia de las dos rutas en el cálculo de Z,(q). A menudo no resulta sencilla la determinación de Z, (q) mediante .1 cálculo del determinantedei tensorSirico covsriaot. métrico covariante O o del .lacobiano de la transformación.Existe un procedimiento alternativo 16 que bate uso del llamado temor métrico contr.variante, E. E es una matriz cuyos terminos se definen como: ‘8~k~ (8qI~ ~kI Za; ti~j ~-~—j CC. 33) 1—1 De acuerdo con las definiciones de It y O, e> producto de ambas matrices, (GE). será: * (GH)kla Z Gth II,, CC. 34) Ii. 1 Bx Bx Bq~ Bq 1 donde es ‘ma swbmatriz cuadrada de G con dimensiones (3*44,) x (314—71,1. siendo 71, el nísero total de iigadurn rígida. del sistema. El cálculo directo de G’(411 resulta notm*le!tt 266 bastante complicado, afortunadamente existe un procedimiento alternativo”’’ 5 que facilita nótablesentela obtención de Z,,Uq). Definiendo la matriz RrCq). de disensiones ~ x 74. que es una subntriz de 71. en la que se incluyen los elementosHu±.dondelos indices Sc y 1 se refieren a las coordenadas ligadas se verifica: 1 G’~ 1 G ¡ CC. 41) Para demostrar la anterior relación se define una matriz auxiliar It’. de dimensiones SN x 371, definida como: CC. 42) 5 18.., 1 4 { q > (C.43) Donde simboliza el conjunto de coordenadasgeneralizadas rígidas, y {q<> el reato de coordenadas. Esquemáticamentepuede representarsela matriz it como: — 1 : ? } donde por convenienciase han ordenado las coordenadas,según sean rígidas o flexibles. Según las definiciones previas it’ es una subaatriz de E de dimension (314—71,) ±4 74,, 0 simboliza una subeatriz de dimensiones N~ x (371—74,) cuyos elementos son todos 267 sinos nulos y, finalmente 1 es una matriz unidad de dimensiones (3N-4~,) (3~I—M,1. De modo análogo pueden desgiosarse la matriz G. El producto de matrices G x it dará lugar, teniendo en cuenta la relación (C.39): cita GR½ ¡ Gr o §TJ Ir ~ :~ ? 1 CC. 45) CC. 46) Por lo que se verificará: ¡~ ~¡ — ‘rl = G ¡ ¡ . ~hia cierta configuración molecular podré ser definida o bien a partir de las coordenadascartesianas de los átomos: 269 Autos o a través del siguiente coordenadas generalizadas: ~ (0.2) It. Y, 2, son las coordenadas de posición molecular. Normalmente se utilizan las coordenadas del centro de masas molecular, aunque en principio cualquier otro punto de referencia definible a partir de las posiciones de los 74 átomos es adecuado. De modo análogo las coordenadas~, e y • son las coordenadasde orientación, asociadas con la rotación molécular. El resto de coordenadas internas se definen ( utilizando notaclón vectorial como: — ¡ R,,, 1j . Y 1 a 71—1 <0.3) r R,,, ~ R,., ¿.2 7 cosO, e [ IR..,,] R,.,,.íi , Y i S N~2 <0.4) cos4, s[ —~::-~---;i-~— ], Vi 174—3 (D. 5> Los vectores auxiliares 5, y 7, se definen como: 270 r ~í. 1 .‘ R,,±‘~a1 — R,.i,~ — ~• IR±.±,,,21 ~ R,,,,, Rl,2. ~.± — [R¡,>,í,:R:, ‘ ~ •.± (D.7) En las anteriores ecuaciones los vectores R1,1 se definen como: ,~m R1 - (D.8) siendo Rj y R, las coordenadasde posición de los centros j e respectivamente. La ecuación (D.5) no permite determinar el ángulo ~,. sino el valor absoluto de éste. La determinación del signo puede realízarse mediante el cálculo del seno del ángulo •,. Por construcción los vectores S~ y T, son ortogonales al vector R1,11.2. Definiendo el vector auxiliar fl1 coso: — x CD. 9) Se obtendrá un vector con la dirección de R±.±±,ssalvo en los casosparticulares en los que coSO, valga ±1 ( ~± y T, con la misma direcciós ). Ea tales casos no existirá aabigúedaden el valor del ángulo •,. Podrá escribirse por tanto: ¡fui Sí ¡ JT~ 1 ¡ sen •~¡ CD. 10) 271 ANEXOS D, E,,, 1.2 sen4, s, ¡ ¡R,., como: — ~ — bw.., cosCo,,.. 2) rs • b,,..~ sen(%2) cos(h,) fi • b,.,sen(e~2) sen(4~~,) 2’ (0.14> Por tanto para la transformaciónde coordenadas: —> (E,.R2 CD. 15) Los elementos de la matriz .7 (anexo C) serán del tipo: (0. 16)~kl — quqí donde los términos xu simbolizan las coordenadasde partida y q± las nuevascoordenadas.De acuerdocon la transformación (0.14) tos elesentos de la matriz 3 en los que intervienen en el denoainadorde la ecuación (0.16) las coordenadasir,verisntes en la transformación serán: ±ulx~ mg1 • O — ~u • Vg1 E Vg, e (D. 17) CD. 19> 273 A±4OOS Por lo que la matriz .3 tendrá la forma: 1 0 0 0 . Bx,, By,, ..fEiiSx, Bx, Bx, OIGO. -2~. By~ _ •oo ío 4 ~ * B,c,, By~ Bz, , .7 0 0 0 1 . B B (D.l9) o o o o . BX~ Sy,, 23.i Bb,<.., ab,,.., Bb,,..í ~ ~ ~ ~ 8y~ BXN Bx~ By~ Bz~ o o O O 84,,.., B4,.., B4,~~ Por la estructura de la matriz .3 se cumple evidentemente: — J,,¡ (0.20) Siendo S~ una submatriz de .1, definida como: Bx,, ay,, ~ 1 =~i 274 Teniendo en cuenta la ecuación (D.14) puede escribirse simbólicamente como: 1< -b cost~ y, b senO b senO cos$ CD.22) sent Donde para simplificar la notación se han suprimido los indices correspondientes a las variables b~..,. e,..2, 4,....s. La matriz de coeficientes no dependede las nuevas coordenadas, por lo que derivando el vector columna, obtenemos la matriz .7 como: producto de dos matrices: fas ‘J~ s~ ‘~ ~f—cosO $~ senO cosa s± i,J ~5enS5en~ b senO b cose coa~ b coso sen~ o —b seno b senO sen~~ cos~J CD. 23) Se cumplirá: — ¡SMI (0.24) Siendo 3,, la matriz transpuesta de .7,,, por tanto el Jacobiano de esta primera transformaciónserá igual a la raíz cuadradadel determinante del la matriz producto de multiplicar .7,, por su transpuesta: a a a ~s. ~ «5 $y «y $, «z st 275 AMEMOS 10.25)— ¡ 3±4 4 1/2 Teniendo en cuenta las relaciones de ortogonalidad entre los ejes auxiliares a, fl y y, el producto de matrices definido anteriormente da lugar a: (—cosO JN~II jbsene 3NJn sene cos4 b cosO cosa —b mene sen~ (-cose senO cos4 (seno sen~ o b2 o senO b cosO b senO sen~’~ costj baena O b coso cos0 -b senO sena[ <0.26) b cosO sena b senO cos#j 0.1 o b2sen2e j <0.27) Por lo que el jacobiano de la transformación indicada en la expresión (0.15) es: 2¡41 — ~—± seno,,. 2 (0.28) LS. modo análogo puede procederse a la realización de transformaciones sucesivas: 276 {E,.ft 2. ..N’...±.bu.,±.On..z.41m..,}> Operando de modo semejante se obtendrá: — b,,..2 senQ~..~ <0.30) Sucesivas transformacionesdarían lugar a un jacobiano global de la forma: Mt ¡ • 4 senehí (0.31) It.’ Para realizar esta transformacionesse han introducido los ángulos auxiliares Oo. ~ y •,~ que por ejemplo pueden de!iuairse utilizando los ejes cartesianos del sistema de 1Sreferencia . Estos ángulos constituyen una de las posibles representaciones de la arieniacion molecular. La posición molecular se define a través de las coordenadas de] primer centro de la cadena. 27? ANEXOS 278 AP~XO E DINAMICA DE SISTEMAS CGI LIGADIRAS En éste anexo se tratará la dinásica de sistemas con ligaduras bolonóaicas’ 4. En estos casos es posible el desarrollo de la mecánica mediante el uso de coordenadas generalizadas. De este modo se puede pasar del tratamiento de un sistema en coordenadas cartesianas con 3M coordenadas y 71. restricciones de ligadura a un problema de dimensionalidad 3>1—71, en coordenadas generalizadas. Serán se vió en el anexo C la transformación de coordenadas cartesianas a coordenadas generalizadas conduce a la expresión de la energía de un cierto sistema conservativo como: fl(LPq> — ti(q) • fl(q,p,,) (¡.1) Donde 3< se expresa como: SN SP SS It.! ‘—± a—a SM SN a. fiq •q k•i I•I e—! 279 Autos Las derivadas con respecto al tiempo < velocidades 1 en ambos sistemas de coordenadasestán relacionadas por: 3M X~ E St (E. 4) g 3M ~ E k, (E. 5) a- i Por otra parte, de acuerdo con la formulación de Lagrange los momentos generalizados p 5 vienen dadis para sistesas conservativos como: aq~ q,q gk LL~ 8~ic.Iqt8~i>q ‘—5 a’t Utilizando las ecuaciones y (E. 10> se deduce: ~ SM SN ax• ~, ag~~ x, E E E II ~ ~ kl i± b~l Ma Xa— ~ ( q,jPg, (E. 12) u—’ Derivando respectoal tiempo la ecuación anterior se obtiene: e, — E ~ E pg,4¡(—~~;] tA~7Jt Ox,> 33 3’. 3% SM aq ‘ 8~d’ • E E E E E~‘ representa las coordenadas rígidas del sistema. Ea sistemasmolecularescon distancias intrasolecularestijas se pueden expresar tales ligaduras mediante coordenadas genetailndas 4,,, Una de las posibles maneras de expresar las coordenadas g,±es: — .4-.. [ 3]si a, Donde la coordenada generalizada rígida q,, mantiene fija la distancia entre los centros a~ y fi,. Derivando respecto a las coordenadas se obtiene: — ~ ~a,$j —8,,,,,, 1 (x~~ X«,> (L29> Donde fi son deltas de Kroneek*r. Sustituyendo este ab resultado en la ecuación (E.26) se obtiene: Nr 35 3M 3M 3M 3 — rrr •NtflqJ,, II~ k.lbtttsSdt CLaC) 285 ‘unos Nr m,i¿,, —L--——-)•Z (4,,fl, «~ 1(x~ —x« 1 (£21) Los valores A vienen dados coleo: AL ax~ AiflLLL c %bH8~r•J~kJj b1 klc±d—t (E. 32) La ecuación (£231) puede escribirse en notación vectorial como: Kl. m 5&~m. ~ <~e.aj~a.$j 1 R«,$, ?tq~ (£4331 Utilizando esta ecuación es posible llevar a cabo la integración de las ecuaciones del movimiento de un sistema con ligaduras empleandocoordenadascartesianas mediante el cálculo de los valores A tales que se verifiquen las condiciones de %. ± ligadura a lo largo de la trayectoria generada. 286 ANEXO r t&ESTREO DE LAS ORIENTACUES MOLECUJLRES £1 muestreo de las orientaciones moleculares Se realizó mediante el siguiente esquema: a> Generación de un eje de giro al azar, uniformemente distribuido en todas las direcciones del espacio. U Generación de un ángulo de giro al azar, <. uniformemente distribuido en el intervalo (—t¿,A4). c) Aplicación del giro así generado a las coordenadas átoticas relativas < respecto del centro de gravedad de la molécula ). El paso • a se consigue mediante el siguiente procedimiento: a.l) Generación de tres números aleatorios (z.y,z) uniformemente distribuidos en el intervalo 1 —1. 1 1. a.2) Cálculo de t, definido coto: 2 2 2 2 CV.. lx Rx•y•z 2 Si R > 1 se retorna al paso a. 1 en cato contrario se prosigue en el paso a.a a.3) El vector director normalizado del eje de giro se define utiiizamdo los cosenos directores 0M’ O’., e 2. 217 ANEXOS g — 1 c,,, c~, c 2 1 c5— ±OR, c,= y/B. c2 zA~ 1 (F. 2) Fi procedimiento utilizado en las simulaciones por Monte Carlo, una vez obtenido el eje de giro g. y el ángulo de rotación ¿ consistió en generar la correspondiente matriz de giro y aplicar ésta a los vectores correspondientes. Sea a — !a,.a,.a1] el vector de un cierto átono con respecto al centro de masas molécular. Las componentes paralela y perpendicular al eje £ serán: a3 — ( ga ) . 5 a1 a- < ga ) g LV.. 3) (F. 4) Definiendo el vector auxiliar b: b — g x a1 — g x a (rS) Se puede expresar el giro como a •a1 •a±cos<•bsenc CV.. 6) 288 Operando se obtiene: a a< + (a—a 11) cose • b sen~ a = a (l—cos¿)+acosg+bsen¿ 2 c1 c,c, a a •!2—cos¿1 c,C, c, 1 c5c1 c,c2 c,c2 )fa5) 100 1Ia~) cyc5 fja, j •cos¿ cío ~ a II ¡ IIIc, >t¾) 1001 o sen¿ ~ O —c2 c, a,1[—c~ c~ O] :: Agrupando términos: a Ga G es la matriz de giro buscada: y 10.2> se obtiene: ________ = (r,,, — r1 3 ( ~ ~í.ol (0.43 (8g20) a (r,.2 — r, 3 1 ¿,.3,a ~íoJ (0.53 291 ANEXOS flonde los a son deltas de Kronecker. Si las masas de todos los grupos de la cadena son iguales los términos podrán escribirse cono: —l 5 ~ ¿j’ 2a ir,,,— r, ¡ a 1r,.~— r 1)(r,.¿— r~,2) o sí ¡1—ii a 2 2 = O 0 Si ¡3—1> 5 3 H¿,.., 11= ~ (~,.± — r, >¿r,,2— r1> —I (0.6) 10.7) (0.8) (0.9) <0.10) (0.11) (0.12) (0. 23) (0.14) 292 —I —7a (r,. 3— r,,2)~1¿a~,2cí.,s a” b 2(cose (1—cose) • sen2O cos (0.24> ligadas al valor b y (0.19> <0.20) (0.21) (0.15) (0.16) (0.17) (0.18) 293 ANEXOS En la ecuación <0.25) se presenta a modo de ejemplo la estructura del determinante ¡11r¡ para el n—hexano con distancias y ángulos de enlace rígidos. H,~ 11,3 11,, 11,, 0 0 0 0 0. 11¿5 H~, 112, 0 112$ 11,6 0 0 0 H,~ 11s, 1135 11,, 11,~ 11,6 0 0 0 H,~ O H,¿ 11<~ 1145 0 H~, 11~ O IHTI= O 11,~ 113$ ~hs 11~g Heé 1157 ~se 0 (0.251 0 1156 ~ 0 (156 1166 ~c, 0 11~ o o 0 11<, 1157 Meí 1177 ~78 3119 O 0 0 ~ ~ O $e 11p 315, o o O O O H~ 11,, ~89 H~ El desarrollo de los determinantes de ¡ 11r1 da lugar, como puede fácilmente intuirse, a expresiones bastante complejas en función de los ángulos de torsión. La probabilidad de una cierta configuración en el subespaciode las posiciones vendria definida para los modelos de ligadura rígida (anexo Cl coco, p(R’.fl’. <*1 )= (0.26) 294 Escribiendo en función de] temor métrico para e] sistema sin ligaduras. ¡G¡. (anexos C y D) y del determinante de la matriz auxiliar 11” se obtiene: ‘a t (0.27) ¡0 1 I~I ¡~<~j ~=1 Donde ¡E’¡,simboliza el determinante de la matriz II asociada con la molécula 1. Una de las consecuenciasmás importantes de la utilización de modelos de ligadura rígida qué incluyen fijación de los ángulos de enlace de la molécula es la complicada dependencia del teMor métrico 1c~¡ con los ángulos de torsión. A menudo para evitar estas dificultades se someten a ligadura rígida exclusivamente las distancias de enlace. incorporándose potenciales de tipo oscilador armónico para la variación de los ángulos de enlace respecto de sus posiciones de equilibrio. Ea tales casos loa deterainantes ¡1(I,tienen una forma más sencilla. por ejemplo para el n-hexano, considerando las masas de todos los grupos iguales: II,, ~ O O O H,¿ H¿~ H2, O O 0 >4,, H,~ 11,. 0 (0.25) o o )f,~ i*~~ )<,s O O O Ng Hgg 295 ANEXOS 2m’ 1,2 (0.26) —I 2 m 1, (1—cosO,) (0.27) Ea este caso no aparece dependencia respecto de los ángulos de torsión en el tensor métrico ¡0 El determinante ¡Mt 1 para el n—butano con ligaduras rígidas. y considerando iguales las masas de todos los grupos tiene la forma: 2 1—c0 c0 c (i—c ) O O O 2 •5 cosa e l—c 4(1—c 3 l—c O c U—c 1 O e e e e 2 •58cO5~ £fl!...I= c0 l—c 2 1—c c (C.2S) 1,10 o o e c0 (l—%) O l—c0 4 i—c0 1: a+a0cos~ 0 c8(l—c9) l—c0 2 +scomt 2 2 Donde c8=coso, 50=sen e y • es el ángulo de rotación interna. Desarrollando el determinante anterior se obtiene: 296 1 lCfr (0.33) (0.34) 297 7. REFERENCIAS (1) P.JFLORY, Statistlcai tfechanics cf Chain Mojecuies. <.JAJiley & Sons, New York 1969) (2] W,L.JOROENSEN, .1. Chem. Phys. . 87, 5304(1923) ¡3) A.R4W4GJ~nTNER, Monte Cario i~rthods Sn Statisticai ?hysics, Toplcs Current Phys. Vol 36. editado por K.BIND~, (Springer—Verlag, Berlin, 1924) (41 M.P.ALLI1~ y D.J.TILDESLEY , Compjter Siwlatlon of Liquids (Clarendon, Oxford, 1987) 5] PO. de GEI4NES, Scaling Concepts in Polymer Physics, (Cornelí Unlverslty Press, lthaca. 1979) 161 Siuaalation of Llquids and Soilds, Molecular Dynamics and Monte Carlo Hethods jo Etatistlcal Physics, editado por G.CICCOTTI, D.fl(ENKEI. y 1.R.NcDONALD (North-Holland, Ajasterdam 1927> 17] .J.P.HANS2~ y (.R,NcDONALD, Theory of Simple Llquids, 2nd edlt ion. (Academic Press, London. 1926) [8] PAMADDEN. Mclecular-Dynaalcs S¡sa¿lation of Statistlcal Mechanlcs systems, Proceedlngs of dic International School of Physics Ronco Fersi, Caurse XCVII. Editado por OCICCOTI Y W.G.HOOVER (North Hollane, Amsterdam, 1926) [9] MALVABEZ, F,J.BERMEJO, W.S,HOWELLS, E.ENCISO, NOALMARZA y M.0ARClA-}iEJR1íAflD~, McI. Phys. 71. 865(1990) 299 (10] F.J.B~HEJO, E.ENCISO, .3.C.DORE, P.C11JEUX, NOARCIA y J.SANT0RO, .1. Che.,. Phys., 87, 7171(198?) H.ALVAREZ, F.J.BERMEJO. WS.HOWELLS, P.CI4IEIJX, BENCíSO. J.ALONSO y N.OARCIA, .1. Ches. ?hys. . 95, 3689(1989) 1113 N.0.AIj*RZA, LENCISO, lAI&MSO, F~J,B~MEID, y 14.AI.VAREZ, Mo). ?hys.. 70, 485(1990) (121 .1.I.SIEPMANN. Mo). Phvs., 70. 1145<1990> 113) D.A.McQUARRIE. Statistical Mechanlos. (Harper & Rov, New York, 1976) ¡14] H.OOLDSTE1N. Mecánica Clásica. Aguilar, MadrId, 1972) (15] L.D.LANDAU y E.M.LIFSHITZ, Mecánica, Vol 1 del Curso de Física Tedrlca. (Reverté, Barcelona, 1970) (161 J.P.RYcJCA~T y 0.CICCOTT. .1. Ches. Phys. , 78, 7362(1983) 1171 J.P.RYcKAERT, Computar Sis*¿la(Icn of Flulds. Polymers erad Solids, RATO AS! Serles C, Vol 293,editado por C.R.CAThOIJ, S.C.?ARKfl y N.P.ALLEN, (Kluver Academio Publishers. Dordrecbt, 1990) (SS] N. 00 y HA. ScllmA0A, Macromolecules, 9. 535(1976) (191 H.C.?F~ENDSEN y W.P.VAYI OIJNSTEREN, Molecular Llquids. Pynamlcs and Interactions, editado por A.J. Barnes, Ni. Orwille—Thomas and J. Yarwood (Reidel. 1984) 300 (20) 1.11>4. GUAÑIR, Computer &iuzaulation of FJulds, Polymers and SofléS. NATO AS! Series C, Vol 292,editado por CRCATLO[4, S.C.P#J1XER y M.P.ALLEN, iNluver Acadeaic Publisbere. Dordrecbt, 19901 [211 0.CHANDL~, !ntroduction lo Modero Statistical Mechanics. [Oxford Universlty Presa, New York, 1987) [221 l.S,ROWLINSON y F.L.S~lNTON, Liquid anó Liquid Mixtures, Thlrd Edil ion, (Butterworth Sclsntific, London 1982> (23) REOBERO. D.J.EVANS and 0.P.HORRISS. .1. ches. Phys.. 84. 6933(1986) [24< l.P.PYCYJflT y A.BELLESIANS, Ches. Phys. Lett., 30. 123(1975) [2S1 IP RYCKAERT y A.EELW4M1S. J. Ches. Soc. Faraday Dlscuss. 66. 95(1978) [261 M.RIGSY, £B.SHITH, W.A.WAXEHAN y G.C.NAITLAND. The Forces between Moiecules. (Clarendon Presa, Oxford, 19861 [271 SL..PR[~E, Compoter Siaajlation of Flulds, Polymers and Sollds, NATO AS! Series C, Vol 293,editado por C.R.CATLOW. S.C.PAÑIE~ y H.P.ALLRN. iKiuwer Academic Publíshere. Dordrecht, 1990> [2(11 MICILLAN, Cosputer Siaujation cf Fluida. Polymers and Solida, NATO AS] Series C, Vol 293.edltado por C.R.CATLOW, S.C.PARKER y M.P.ALLEN, rKluwer Acadesio Rublishers. Dcrdrecht. 1990) [291 ROAR y M.PARRINELLO, Phys. Rey, Lett. , 22, 2471(1985> 301 120) R.CA? y 14.PAlRRlNaLo. SInsp)e Molecular Systems al Very Hig}s Oensity. editado por A. POLIAN [Flenum, 1989) [21) N.PARRINELLO, Conferencia en el curso: NATO AS! on Computar Simulation Sn Materials Soleose, Interatomlc ?otentiais, Techniquesand Applicatlons, Aussois [Francia). 1991 (22] D,IC.Rfl!LER y PA. MADOEN, McI. Ph,s. , 70. 921<1990) 122] R.G.PARR y W.YANG, Density~-Functional Theory of Atoas and Molecules. (Clarendon Presa, Oxford 1989) [34] PHORENBERG y W400N,Phys. Reo. 2, 136. 864(1964> 35] W.K0HN y LISHAN. Phys Reo A. 140. 1232(1965) 361 l.N.LEVINE, Qulsica Cuántica, [Editorial AC, Nadrid 1977) [27) 14, L.IOIIITNSEN. R.C.RIN1I)NO y 3.91001. J. N.a. Ches,. Soc. 103. 4392(1981) [23]J.i.LJLLO y S.YIP, ~1. Che,,. Phys. . 85, 4056(1986) 29] D.FPENXEL, tfolecular—Dvnamlcs Sierúlation of Statlstical Mechanics Systems, editado por C.CICCOTTI y I4.GI400VER (North Holland, Amsterdam 1986 [401 lIFRENREL, Computar Sie,ulation of Fluids, Polymers aoci SoIlds, NATO AS] Series C. Vol 29J.mditmdo por C.R.CATh014, S.C.P.kR] [46>M. PAPRINELLO y A. RAHMA14, .1. AppI. Phys. , 52. 7192>1981) [47) E. LOMBA, Tesis Doctoral. S.NOSE, 2. Ches,. ?hys. 81, 511(19841 [52! S.N0SE, Conferencia en el curso: RATO AS! on Cos,puter Simulation ir Matcrjals Science, Interatomio Potentials. Techniques ar4 Applications, Aussois (Francia), 1991 >53> 1411PRESS, B.P.FLANNERY, S,A.TEUKOLSKY. 14. T.VETTERLING, “Numericai Reclpes, The Art of Scientlfic Cos,puting, [Canbridge University Press, Cambridge 1996) >541 iP VAILEAU y S.0.I4EITTINOTON, Sial istlcal Mechanlcs, pan .4, editado por B.í.BERNE, [Plenus,,New York, 1977) 303 >55) 0.0LOWRY (editor), tlarkcv Chalns anó Monte Carlo Calculatlons lo Fol yaer Selence [ Dekker, New York 1970> (56! NOVAN KAMiFEN. Stcchastlc ?rocesses In Physics .and Chealstry, North—Holland, krnsterdaa 1981> jS?] N.1ffTROPOLIS, A,M.WOSE34BLUTII. *14 ~OSE1~8L\lTM, A.itTElI1~ y E.TELLER, .7. Ches,. Phys. 21, Th87[1953[ (59! tIRAD, C.PAPGALI y B.i.8FRNE. Mc!. Fhys. . 37, 1773(1979) [59] T.N~UT] y NCC, Biopolyaers. 24, 527(1985) [60) 31.0>001 y WL.IORGENSEN, 3. Ches. Phys.. 75. 1944>1981> [61) DRESFRIUS, E,.).BERNE y D.CHANDLFR, 2. Che,,. Phvs. 70, 3395(1979) [62] P.K.MEJ~1OTRA, M.M~E1 y kL BEVERIOGE, 1. Ches,. Ph>’s., 79. 3156(1983) [631 PJ.ROSSKY, í.D.DOLL y HLF’RIF.DMAN, ,1. Che,,. Fhys. . 69, 4628(1978> 64] SGOLDMANN, .1. Ches,. Phys. 79, 3938(1983) [65] tIMEZEI, KA.BENCSATH, SG0LDMAN y 5.9114GB. Molecular Siniulation, 1. 87>19871 (66] l.C.0W[CR] y H.A.SCHERACA, Che,,. Phys. Lett. 47. 600(1977) [67! .JC.OWICXI y 11.A.SCHERAGA, .V les. Ches,. Soc. , 99, 7413(1977> [69] .7.í.S[EPMANN, [Cosnunicaciénprlvada) 304 >69) NC ALNAR2A y E ENCISO (en preparacIón) >70] A.l.C.L#flD, Cespitar Siavlat Ion of Fluid,. Poiymers and So)ids, NATO ASí Series C. Vol 293, editado por C.R.CATLOI4. SC.PARJC~ y M.PALLEN. CK(uwer Acadenio Publishera, Oordrecht, 1 990) [71) ABARANYA[ y 0. .7. EVAMS. Piel. Phys., 70. 53(1990) [72] .7,PRYGKÁERT, C.CICC0TT[, y 11.í.CBERENDSEN. .7. Coaput. Phys. , 23. 327(1977) [73> (f.C ANDERSEN, 1. Conput. Phys., 52. 24(1983) [74) 1. ALONSO. Fi. BERMEJO. >4. OABC[A—11ERNANDEZ, 2. E. MABT[NEZ y VS HOWELLS, .1. Moiec. Struc. (7) 1991 [75) .7 ALONSO(Cornunlcac.tón privada) [76] C.H.BENNET. .7. Comput. ?hys. , 22, 245(1976) [77) 5. 11. NORT11RUP y l.A.HcCAÑION, Eiopolymers, 19, 1001(1980) >78) 0. BROWN y .1 MR. CLAME. .1. Che.,. Phys. , 92. 3062(1990) >79> V.P,SPIRIDINOV y A.A.L0PATKIN, Tratamiento sates,ático de datos flsicc—quialcos, (Editorial MI?. HoscO, 1973) 1801 E. ENCISO, 3. ALONSO, N. 0. ALMABZA, y Y. 3. PEPJ4E.7D, .1. Chem. Phys. , 90. 412>1989), E. ENCISO. 1. ALONSO. tI. 0. AIJ4ARZA. y E. .7. BERMEID, 1. Chem. Phys. , 90. 422(1969), 305 [E)] O. CHMJDLER y LR.PRATT, .7. Ches. Phys. . 65. 2925(1976) (82) A. BÁÑON, Y. SERRANO-ADÁN y 1. SANTAYARIA, .1. chea. Phys. , 83. 297<1985) >831 W.L.JORGENSEN, .1. Osem. ?hys.. 77, 5757<1992> (94] W.L.2OROENSEN, .1. Am. Ches,. Sor.. 1D3. 4721(1981) [SS] W.L,JOROENSEN, J.D,HAOUPA, y CA.SWENSON. .~. Am. Chem. Sor. 106. 6638(1984) (96] T.A.WEBER, 5. Che,,. Phys. , 69, 2347(1978> 87) S.TGXVARRO. .1. Chem. Fhys. . 89, 3008(1999> 1893 R.EBBE3G, G.P.M&RTUSS y D.3.EVAJ4S, .5. Ches,. Phya.. 86, 4555(1987) 89) P. A. W[ELOPOSK[ y EA4.SMITH, .1. Che,,. Phys. 84, 6940(1986) >90) 5 0. LECOEPTER y O.i.7[LDESLEY, Mo]. Phys. , 68, 519>1989> (91) K.0.HONNELL y C.KRALL, en: £quations of State: Theory and Application, ACS Symposiua series, N 300. [American Chemical Society, Washington DC, 1996> >92> HBDW[GRT Tahies of InX estaje and other Piathemat.i ca] Data, [Macmillan Cospany, New York, 1954) [93] M.BISHOF. D.CEPERLEY, H.L.F’R>SCH y M.H.KALOS. ./. Chem. Phys. 72. 3228(2990) 106 [94> M.KKIAiAB. M.VACATE?uLQ y O.Y.YOON. .1. Ches,. Phys. . 89, $20611988 [95>NGALMAPZA, EENC>SO y E.3REBHEJO, 2. Chem. Ph~., [7), 1991 >96> ALOPEZ RODR>GUEZ, 2.VEGA. .LlEREIRE, y S.LAGO,1991, Mol. Phys <71, 1991 [97> FLS..ACTON, Nu,,erlca] Nethods that Wcrk, (Jiarper & Bow, New York, 1970) >931> S.TOXVAERD y P.PÁD>LLA, .1, Ches,. ?hys (1991) [99> CAREXEN, Mathe,,atical flethods for Physicists, (Academic Preas, New York, 1970) >100> RD,HAPRISOM [editor), flook of Dala, Chemlstry. Physical Science, Physics. ttIuffield Advanced Science—Longman, London. 1972> [011 TBOUBLIK, 2. Che.,. ?hys., 63, 4084U975) >102) C. VEGA, Tesis Doctoral, Universidad Complutense, Madrid 19911 >103> >1EROW1~ [Conunicación privada. 1991) >104! .LD.C~ y OPIEGNER, Thers,ocheaistry of Organic and Qrgano,,ete]lic Co,,pounds. (Acadenic Presa, London. tIew York, 1970) [05! 1. t.SIIYMANN y 0EBflJKfl. Mo!. Phys. . (7) 1991 307 >106] M.N.ROSENBLUTI{ y á.W.ROSENBLUTIt.3. Che,,. Phys. , 23, 356>1955) >107> S.TOXVAERD. .7. Ches,, Phys. . 93. 4=90(1990) 1108! W.F.¶JRPHY, l.M.PERNA=JDEZ-SANCHIYy K.RAG>4AVNRACHI. 2. Phys. Ches,., 95, 1124(1991) SIMULACIÓN DE LÍQUIDOS MOLECULARES FLEXIBLES ÍNDICE 1. INTRODUCCIÓN 1.1. MECÁNICA ESTADÍSTICA CLÁSICA 1.2. COORDENADAS lNTERNAS 2. MODELOS DE POTENCIAL 2.1 MODELO SIMPLIFICADO PARA N-ALCAMOS 2.2 MODELO SIMPLIFICADO PARA 1.2-DICLOROETANO 2.3 MODELO REALISTA DE N-BUTANO 3. TÉCNICAS DE SIMULACIÓN 3.1 MÉTODOS DE MONTE CARLO 3.2 DINÁMICA MOLECULAR 3.3 ANÁLISIS DE CONFIGURACIONES. TÉCNICAS UMBRELLA 3.4 ESTIMACIÓN DE ERRORES 4. SIMULACIÓN POR MONTE CARLO DE MOLÉCULAS FLEXIBLES 4.1 SIMULACIÓN EFICIENTE DE N-BUTANO 4.2 MÉTODOS DE REPTACIÓN 5. RESULTADOS 5.1 N-BUTANO LÍQUIDO (MODELO SIMPLIPICADO) 5.2 SIMULACIÓN DE 1.2-DICLOROETANO 5.3 SIMULACIÓN DE N-ALCANOS 5.4 N-BUTANO LÍQUIDO (MODELO DETALLADO) 6. CONCLUSIONES 7. ANEXOS A DERIVACIÓN DE FUERZAS INTRAMOLECULARES EN SISTEMAS CON LIGADURAS RÍGIDAS B. CÁLCULO DE PRESIÓN C. TRANSFORMACIÓN DE LA INTEGRAL DE CONFIGURACIÓN D. TRANSFORMACIONES DE COORDENADAS E. DINÁMICA DE SISTEMAS CON LIGADURAS F. MUESTREO DE ORIENTACIONES MOLECULARES H. CÁLCULO DEL TENSOR MÉTRICO PARA CADENAS CON LIGADURAS RÍGIDAS 8. REFERENCIAS 1: