Propiedades-reologicas-e-interfaciales-de-fluidos-complejos-mediante-simulacion-molecular

Vista previa del material en texto

UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO 
PROGRAMA DE MAESTRÍA Y DOCTORADO EN INGENIERÍA 
 QUÍMICA – POLÍMEROS 
 
 
 
PROPIEDADES	
  REOLÓGICAS	
  E	
  INTERFACIALES	
  DE	
  FLUIDOS	
  
COMPLEJOS	
  MEDIANTE	
  SIMULACIÓN	
  MOLECULAR 
 
 
 
Tesis 
que para optar por el grado de: 
DOCTORA EN INGENIERÍA 
 
 
	
  
PRESENTA 
M.	
  en	
  I.	
  PATSY	
  VERÓNICA	
  RAMÍREZ	
  GONZÁLEZ	
  
 
 
Tutor principal: 
Dr.	
  Sergio	
  E.	
  Quiñones	
  Cisneros,	
  IIM 
	
  
Comité tutorial: 
Dr.	
  Enrique	
  Geffroy	
  Aguilar,	
  IIM 
Dr.	
  Ulrich	
  K.	
  Deiters,	
  Universidad	
  de	
  Colonia,	
  Alemania 
	
  
 
	
  
México	
  D.F.,	
  Noviembre	
  de	
  2013	
  
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
 
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
 
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
 
 
2	
   	
  
	
  
	
  
 
 
 
 
 
 
 
JURADO ASIGNADO: 
	
  
	
  
	
  
	
  
Presidente: Dr. José Roberto Zenit Camacho 
	
  
Secretario: Dr. Juan Pablo Aguayo Vallejo 
	
  
Vocal: Dr. Enrique Soto Castruita 
	
  
1 er. Suplente: Dra. Jaqueline Quintana Hinojosa 
	
  
2 d o. Suplente: Dr. Sergio E. Quiñones Cisneros 
	
  
	
  
	
  
	
  
	
  
	
  
Cd. Universitaria, D.F. a de de 2013 
	
  
	
  
	
  
	
  
	
  
	
  
TUTOR DE TESIS: 
	
  
Dr. Sergio E. Quiñones Cisneros 
	
  
	
  
	
  
-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐-­‐	
  
FIRMA 
	
  
	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  3	
  
	
  
Agradecimientos	
  	
  
	
  
A	
   mi	
   asesor	
   el	
   Dr.	
   Sergio	
   Quiñones	
   Cisneros,	
   por	
   enseñarme	
   a	
   ser	
   original	
   y	
   no	
  
permitirme	
  nunca	
  copiar	
  a	
  los	
  demás,	
  siempre	
  hacer	
  las	
  cosas	
  por	
  mi	
  misma	
  y	
  pensar	
  
más	
  allá	
  de	
  la	
  imaginación,	
  por	
  prestarme	
  toda	
  su	
  atención	
  y	
  preocuparse	
  siempre	
  por	
  
mi	
  conocimiento	
  y	
  crecimiento.	
  
To	
  Dr.	
  Ulrich	
  K.	
  Deiters,	
  because	
  you	
  have	
  been	
  a	
  marvelous	
  teacher	
  for	
  me	
  with	
  the	
  
largest	
  humility	
  I	
  have	
  ever	
  met	
  in	
  one	
  person.	
  	
  
To	
  Dr.	
  Thomas	
  Kraska,	
  because	
  every	
   time	
   I	
   required	
  your	
  help	
  you	
  always	
  had	
   the	
  
time	
  and	
  the	
  mood.	
  
Al	
   Dr.	
   Roberto	
   Rozas,	
   porque	
   fuiste	
   el	
   primero	
   que	
   me	
   enseñó	
   a	
   programar	
   y	
   a	
  
aprender	
   dinámica	
   molecular.	
   Me	
   tuviste	
   mucha	
   paciencia	
   y	
   compartiste	
   conmigo	
  
muchos	
  trucos.	
  
Al	
  Dr.	
  Enrique	
  Soto	
  Castruita	
  por	
  compartir	
  conmigo	
  sus	
  experiencias	
  y	
  conocimientos	
  
en	
  reología	
  y	
  ayudarme	
  a	
  tener	
  confianza	
  en	
  los	
  equipos.	
  
Al	
   Dr.	
   Roberto	
   Zenit	
   Camacho,	
   Dr.	
   Juan	
   Pablo	
   Aguayo	
   Vallejo	
   y	
   a	
   la	
   Dra.	
   Jaqueline	
  
Quintana	
  Hinojosa	
  por	
  sus	
  consejos	
  valiosos	
  en	
  la	
  revisión	
  de	
  esta	
  tesis.	
  
Al	
   Dr.	
   Enrique	
   Geffroy	
   porque	
   siempre	
   es	
   muy	
   accesible	
   conmigo	
   y	
   me	
   ha	
   dado	
  
consejos	
  bastante	
  útiles.	
  	
  
Agradezco	
   también	
  el	
   apoyo	
  de	
   la	
   fundación	
  para	
   la	
   investigación	
  alemana	
   (DFG	
  De	
  
391/331),	
   y	
   el	
  Consejo	
  Nacional	
  de	
  Ciencia	
  y	
  Tecnología	
  de	
  México	
   (CONACYT)	
   con	
  
número	
  de	
  beca	
  	
  210692	
  .	
  
	
  
	
  
	
  
	
  
4	
   	
  
	
  
A	
  mi	
  mamá	
  por	
  ser	
  mi	
  apoyo	
  en	
  todo	
  momento	
  y	
  ayudarme	
  cuando	
  más	
  lo	
  necesito.	
  
Por	
  ser	
  mi	
  mejor	
  consejera	
  y	
  mi	
  mejor	
  amiga.	
  	
  
	
  
A	
   mi	
   papá	
   por	
   esa	
   insistencia	
   en	
   mi	
   educación,	
   esa	
   perseverancia	
   por	
   alentarme	
  
siempre	
  a	
  seguir	
  estudiando.	
  Por	
  el	
  gran	
  apoyo	
  sentimental	
  y	
  espiritual.	
  	
  
	
  
A	
  mi	
  hermano,	
  por	
  ser	
  mi	
  guía	
  intelectual	
  y	
  un	
  amigo	
  en	
  quien	
  siempre	
  confiar.	
  
	
  
A	
  mi	
  esposo	
  por	
  alentarme	
  a	
  continuar	
  mis	
  estudios	
  y	
  no	
  permitir	
  que	
  me	
  diera	
  por	
  
vencida	
   en	
   los	
  momentos	
  más	
   difíciles.	
   	
   Y	
   por	
   supuesto	
   a	
  mi	
   hijo	
   que	
   es	
  mi	
  mayor	
  
motivo	
  para	
  seguir	
  adelante.	
  
	
  
A	
  mi	
   familia	
  y	
  amigos,	
  por	
  ese	
  gran	
  apoyo	
   invaluable	
  en	
  algún	
  momento	
  de	
  mi	
  vida	
  
durante	
  el	
  doctorado.	
  	
  	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  5	
  
	
  
CONTENIDO	
  
	
  
RESUMEN	
   7	
  
ANTECEDENTES	
   8	
  
MOTIVACIÓN	
   9	
  
OBJETIVOS	
   10	
  
JUSTIFICACIÓN	
   11	
  
	
  
CAPÍTULO	
  1	
   12	
  
DINÁMICA	
  MOLECULAR	
   12	
  
1.1	
   Qué	
  es	
  la	
  dinámica	
  molecular	
   12	
  
1.2	
   Convención	
  de	
  imagen	
  mínima	
  y	
  condiciones	
  de	
  frontera	
  periódicas	
   14	
  
1.3	
   Potencial	
  de	
  Lennard-­‐Jones	
   15	
  
1.4	
   Posiciones	
  iniciales	
   16	
  
1.5	
   Radio	
  de	
  corte	
  y	
  lista	
  de	
  vecinos	
   18	
  
1.6	
   Velocidades	
  iniciales	
  y	
  temperatura	
   19	
  
1.7	
   Fuerzas	
  intermoleculares	
  y	
  fuerzas	
  de	
  enlace	
   20	
  
1.8	
   Diseño	
  de	
  un	
  programa	
  de	
  simulación	
  molecular	
   21	
  
1.9	
   Ecuaciones	
  de	
  movimiento	
   22	
  
1.10	
   Compresión	
  isotérmica	
   23	
  
1.11	
   Integración	
  completa	
   23	
  
1.12	
   Clasificación	
  de	
  la	
  dinámica	
  molecular	
   26	
  
	
  
CAPÍTULO	
  2	
   28	
  
CÁLCULO	
  DE	
  PROPIEDADES	
  FISICOQUÍMICAS	
  E	
  INTERFACIALES	
  DE	
  CADENAS	
  LINEALES	
  
DE	
  ÁTOMOS	
   28	
  
2.1	
   Equilibrio	
  y	
  estabilidad	
  de	
  los	
  sistemas	
  termodinámicos	
   28	
  
2.2	
   Cálculo	
  de	
  equilibrio	
  de	
  los	
  sistemas	
  termodinámicos	
   34	
  
2.3	
   Equilibrio	
  de	
  fases	
   34	
  
2.3.1	
  Perfil	
  de	
  densidades	
   37	
  
2.3.2	
  Diagrama	
  de	
  fases	
  T-­‐V	
   38	
  
2.4	
   Presión	
   42	
  
2.4.1	
  Diagrama	
  de	
  fases	
  P-­‐V	
   43	
  
2.5	
   Potencial	
  químico	
   49	
  
6	
   	
  
	
  
2.5.1	
  Diagrama	
  de	
  fases	
  µ-­‐V	
   51	
  
2.5.2	
  Diagrama	
  de	
  fases	
  µ-­‐P	
   54	
  
2.6	
   Tensión	
  interfacial	
   63	
  
	
  
CAPÍTULO	
  3	
   66	
  
CÁLCULO	
  DE	
  PROPIEDADES	
  ESTRUCTURALES	
  PARA	
  CADENAS	
   66	
  
3.1	
  	
  Función	
  de	
  distribución	
  radial	
   66	
  
	
  
CAPÍTULO	
  4	
   72	
  
CÁLCULO	
  DE	
  PROPIEDADES	
  REOLÓGICAS	
  PARA	
  CADENAS	
   72	
  
374.1	
  Dinámica	
  molecular	
  de	
  no	
  equilibrio	
   72	
  
4.2	
   Cálculo	
  de	
  viscosidad	
   72	
  
4.3	
   Tensor	
  de	
  esfuerzos	
   82	
  
4.3.1	
   Diferencia	
  de	
  esfuerzos	
  normales	
   87	
  
	
  
CONCLUSIONES	
   92	
  
	
  
APÉNDICE	
  1.	
  Variables	
  adimensionales	
  de	
  interés	
   94	
  
APÉNDICE	
  2.	
  Posiciones	
  iniciales	
   95	
  
APÉNDICE	
  3.	
  Posiciones	
  iniciales	
  de	
  cadenas	
   96	
  
APÉNDICE	
  4.	
  Lista	
  de	
  vecinos.	
   97	
  
APÉNDICE	
  5.	
  Velocidades	
  iniciales	
  ytemperatura.	
   98	
  
APÉNDICE	
  6.	
  Fuerzas	
  intermoleculaes	
  y	
  de	
  enlace.	
   99	
  
APÉNDICE	
  7.	
  Ecuaciones	
  de	
  movimiento	
   100	
  
APÉNDICE	
  8.	
  Compresión	
  isotérmica	
   101	
  
APÉNDICE	
  9.	
  Tablas	
  de	
  presión	
  y	
  potencial	
  químico	
  para	
  cadenas	
   102	
  
	
  
BIBLIOGRAFÍA	
   106	
  
	
  
	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  7	
  
	
  
RESUMEN	
  
	
  
En	
   este	
   trabajo	
   se	
   diseñó	
   y	
   construyó	
   un	
   programa	
   computacional	
   para	
   calcular	
  
propiedades	
   fisicoquímicas,	
   interfaciales	
  y	
   reológicas	
  de	
   cadenas	
   lineales	
  de	
  átomos	
  
de	
  la	
  misma	
  especie.	
  Estas	
  cadenas	
  incluyen	
  tamaños	
  de	
  1,	
  2,	
  4,	
  8	
  y	
  16	
  átomos	
  unidos	
  
por	
  un	
  enlace	
   tipo	
  resorte.	
  Las	
  cadenas	
  son	
   totalmente	
   flexibles	
  y	
  no	
  se	
   tomaron	
  en	
  
cuenta	
   ángulos	
   de	
   torsión.	
   Las	
   propiedades	
   calculadas	
   incluyen	
   equilibrio	
   de	
   fases,	
  
presión,	
   potencial	
   químico,	
   tensión	
   interfacial,	
   función	
   de	
   distribución	
   radial,	
  
viscosidad	
   y	
   tensor	
   de	
   esfuerzos.	
   La	
   técnica	
   utilizada	
   para	
   el	
   cálculo	
   de	
   estas	
  
propiedades	
   es	
   la	
  Dinámica	
  Molecular.	
   Se	
   estudió	
   el	
   efecto	
   que	
   tiene	
   la	
   longitud	
   de	
  
cadena	
   sobre	
   las	
   propiedades	
   fisicoquímicas	
   y	
   se	
   encontró	
   buena	
   concordancia	
   con	
  
los	
   datos	
   reportados	
   en	
   la	
   literatura.	
   En	
   la	
   parte	
   de	
   propiedades	
   reológicas	
   se	
  
encontró	
  que	
   es	
   posible	
   predecir	
   la	
   aparición	
  de	
   la	
   primera	
  diferencia	
   de	
   esfuerzos	
  
normales	
  y	
  el	
  adelgazamiento	
  como	
  una	
  función	
  de	
  la	
  longitud	
  de	
  cadena.	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
8	
   	
  
	
  
ANTECEDENTES	
  
	
  
Durante	
   las	
   pasadas	
   décadas,	
  muchas	
   publicaciones	
   han	
   reportado	
   las	
   propiedades	
  
del	
  fluido	
  de	
  Lennard-­‐Jones,	
  como	
  los	
  libros	
  de	
  Allen	
  Tildesley1	
  y	
  Daan	
  Frenkel2.	
  Este	
  
fluido	
   se	
   ha	
   estudiado	
   extensivamente	
   con	
   la	
   ayuda	
   de	
   métodos	
   por	
   simulación	
  
computacional	
   –	
   como	
   técnicas	
   de	
   dinámica	
   molecular	
   o	
   Monte	
   Carlo3	
   –	
   así	
   como	
  
también	
  por	
  métodos	
   estadísticos	
   como	
   la	
   teoría	
   de	
   la	
   perturbación4.	
   La	
   cantidad	
   y	
  
calidad	
  de	
   las	
  propiedades	
   termofísicas	
  obtenidas	
  por	
  estos	
  métodos	
  permitieron	
   la	
  
construcción	
  de	
  ecuaciones	
  de	
  referencia	
  para	
  el	
   fluido	
  de	
  Lennard-­‐Jones,	
   como	
  por	
  
ejemplo	
  en	
  el	
  trabajo	
  de	
  Johnson	
  et	
  al5.	
  
Para	
   el	
   fluido	
   de	
   dímeros	
   de	
   Lennard-­‐Jones,	
   la	
   situación	
   de	
   los	
   datos	
   obtenidos	
   es	
  
igualmente	
   buena,	
   y	
   por	
   lo	
   tanto	
   ecuaciones	
   de	
   estado	
   también	
   se	
   han	
   construido6.	
  
Posteriormente,	
  se	
  han	
  realizado	
  simulaciones	
  computaciones	
  para	
  fluidos	
  de	
  cadenas	
  
flexibles	
  de	
  Lennard-­‐Jones7,8.	
  Sin	
  embargo,	
  el	
  número	
  de	
  publicaciones	
  de	
  este	
  tipo	
  de	
  
fluidos	
  es	
  mucho	
  menor	
  que	
  para	
  monómeros	
  y	
  dímeros.	
  Las	
  propiedades	
  obtenidas	
  
de	
  estas	
  simulaciones	
  son	
  sobretodo	
  datos	
  de	
  pVT	
  y	
  equilibrio	
  de	
  fases	
  y	
  de	
  nuevo,	
  se	
  
han	
  podido	
  construir	
  ecuaciones	
  de	
  estado	
  de	
  los	
  resultados	
  de	
  simulación5.	
  
No	
  obstante,	
  al	
  parecer	
  existen	
  propiedades	
  termofísicas	
  de	
  los	
  fluidos	
  de	
  cadenas	
  de	
  
Lennard-­‐Jones	
  que	
  aún	
  no	
  se	
  han	
  publicado,	
  como	
  el	
  potencial	
  químico,	
  al	
  menos	
  no	
  
para	
  toda	
  la	
  región	
  de	
  interés	
  que	
  va	
  de	
  bajas	
  a	
  altas	
  densidades.	
  	
  
Respecto	
   a	
   las	
  propiedades	
   reológicas,	
   la	
  mayoría	
  de	
   los	
   trabajos	
   encontrados	
  en	
   la	
  
literatura	
  se	
  refieren	
  a	
  las	
  ecuaciones	
  de	
  Green-­‐Kubo1.	
  Más	
  aún,	
  recientemente	
  se	
  ha	
  
publicado	
   una	
   nueva	
   técnica	
   basada	
   en	
   dinámica	
   molecular	
   de	
   no	
   equilibrio	
   que	
  
permite	
  calcular	
  la	
  viscosidad	
  a	
  diferentes	
  rapidez	
  de	
  corte9.	
  Con	
  esta	
  técnica	
  existe	
  un	
  
amplio	
   campo	
   para	
   explorar	
   propiedades	
   reológicas	
   que	
   aún	
   no	
   se	
   han	
   reportado,	
  
como	
  el	
  cálculo	
  del	
  tensor	
  de	
  esfuerzos.	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  9	
  
	
  
MOTIVACIÓN	
  
	
  
El	
  comportamiento	
  reológico	
  de	
  fluidos	
  no	
  newtonianos	
  ha	
  sido	
  siempre	
  de	
  interés	
  a	
  
nivel	
  industrial.	
  Por	
  ejemplo	
  la	
  extrusión	
  de	
  un	
  polímero,	
  la	
  extracción	
  de	
  petróleo	
  o	
  la	
  
preparación	
  de	
  una	
  mayonesa.	
  Estos	
  fluidos	
  se	
  comportan	
  de	
  manera	
  distinta	
  al	
  agua	
  
por	
   ejemplo,	
   que	
   es	
   un	
   fluido	
   newtoniano.	
   Si	
   se	
   agita	
   en	
   un	
   recipiente	
   un	
   fluido	
  
newtoniano	
   y	
   en	
   otro	
   recipiente	
   un	
   fluido	
   no	
   newtoniano	
   se	
   observarán	
  
comportamientos	
  diferentes,	
  como	
  el	
  mostrado	
  en	
  la	
  siguiente	
  figura.	
  
	
  
	
  
	
  
El	
  comportamiento	
  del	
   segundo	
  recipiente	
  se	
  conoce	
  como	
  efecto	
  Weissenberg	
  y	
  ha	
  
sido	
   objeto	
   de	
   muchos	
   estudios	
   reológicos10.	
   El	
   análisis	
   de	
   la	
   deformación	
   de	
   este	
  
fluido	
  se	
  puede	
  llevar	
  a	
  cabo	
  con	
  la	
  ayuda	
  de	
  un	
  reómetro.	
  Sin	
  embargo,	
  qué	
  es	
  lo	
  que	
  
está	
   ocurriendo	
  desde	
  dentro	
  del	
   fluido,	
   a	
   nivel	
  molecular,	
   es	
   la	
  motivación	
  de	
   este	
  
trabajo.	
   Lo	
   que	
   se	
   pretende	
   es	
   entender	
   qué	
   sucede	
   con	
   las	
   cadenas	
   moleculares	
  
cuando	
  son	
  sometidas	
  a	
  esfuerzos	
  y	
  cómo	
  se	
  ve	
  reflejado	
  este	
  comportamiento	
  a	
  nivel	
  
macroscópico.	
  	
  
	
  
	
  
10	
   	
  
	
  
OBJETIVOS	
  
	
  
Con	
  base	
  en	
  los	
  antecedentes	
  y	
  motivaciones,	
  el	
  objetivo	
  principal	
  de	
  este	
  trabajo	
  es	
  
desarrollar	
   una	
   metodología,	
   basada	
   en	
   dinámica	
   molecular,	
   que	
   permita	
   calcular	
  
propiedades	
   fisicoquímicas,	
   interfaciales,	
   estructurales	
   y	
   reológicas	
   de	
   sistemas	
  
constituidos	
  por	
  cadenas	
  lineales	
  de	
  átomos	
  de	
  la	
  misma	
  especie.	
  	
  
	
  
Dentro	
  de	
  estos	
  objetivos,	
   se	
   calculará	
  el	
  potencial	
  químico	
  y	
  el	
   tensor	
  de	
  esfuerzos	
  
para	
   cadenas	
   lineales,	
   resultados	
   que	
   a	
   la	
   fecha	
   no	
   se	
   encuentran	
   reportados	
   en	
   la	
  
literatura.	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  11	
  
	
  
JUSTIFICACIÓN	
  
	
  
Las	
   propiedades	
   fisicoquímicas	
   son	
   fundamentales	
   para	
   la	
   comprensión	
   de	
   los	
  
fenómenos	
  que	
  ocurren	
  en	
   la	
  naturaleza,	
   tanto	
  a	
  nivel	
  de	
   investigación	
   como	
  en	
   las	
  
aplicaciones	
  donde	
  se	
  requieren	
  para	
   lograr	
  un	
  diseño	
  funcional	
  o	
   incluso	
  optimizar	
  
procesos.	
  Pero	
  no	
  siempre	
  es	
  posiblemedir	
  dichas	
  propiedades	
  de	
  manera	
  directa	
  o	
  
es	
   muy	
   costoso	
   realizar	
   las	
   mediciones	
   de	
   forma	
   experimental,	
   por	
   ello	
   se	
   ha	
  
desarrollado	
   la	
   dinámica	
   molecular	
   que	
   por	
   medio	
   de	
   simulaciones	
   numéricas	
  
permite	
   obtener	
   una	
   buena	
   aproximación	
   del	
   comportamiento	
   real	
   de	
   los	
   sistemas	
  
físicos.	
  
	
  
Existen	
  disponibles	
  programas	
  o	
   software	
  que	
  permiten	
  realizar	
  modelado	
  o	
  diseño	
  
molecular	
  a	
  nivel	
  académico	
  y	
  profesional.	
  Sin	
  embargo,	
  este	
  trabajo	
  está	
  basado	
  en	
  la	
  
construcción	
   de	
   un	
   programa	
   propio	
   que	
   tiene	
   como	
   justificación	
   el	
   obtener	
  
información	
   de	
   una	
   naturaleza	
   fundamental	
   y	
   ser	
   una	
   semilla	
   para	
   futuras	
  
investigaciones	
  en	
  el	
  tema,	
  como	
  lo	
  son	
  las	
  cadenas	
  ramificadas	
  de	
  átomos.	
  	
  
	
  
El	
  modelo	
  deberá	
  predecir	
  el	
  comportamiento	
  observado	
  en	
  materiales	
  conocidos	
  con	
  
la	
   finalidad	
   de	
   validar	
   el	
   método.	
   De	
   tal	
   modo	
   que	
   se	
   compararán	
   resultados	
  
obtenidos	
  en	
  la	
  literatura	
  con	
  los	
  obtenidos	
  con	
  nuestro	
  código.	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
12	
   	
  
	
  
CAPÍTULO	
  1	
  
	
  
DINÁMICA	
  MOLECULAR	
  	
  
	
  
En	
   este	
   capítulo	
   se	
   explicará	
   detalladamente	
   el	
   concepto	
   de	
   dinámica	
  molecular.	
   Se	
  
presentarán	
   los	
   principales	
   conceptos	
   que	
   se	
   deben	
   tomar	
   en	
   cuenta	
   para	
   una	
  
simulación,	
   como	
   el	
   potencial	
   de	
   interacción	
   y	
   condiciones	
   de	
   frontera.	
   También	
   se	
  
mostrarán	
   las	
   ecuaciones	
  principales	
  y	
   fuerzas	
  que	
   rigen	
   todo	
  el	
  movimiento	
  de	
   las	
  
partículas	
   durante	
   la	
   simulación.	
   Importantes	
   trucos	
   computacionales	
   se	
   muestran	
  
para	
  ahorrar	
  tiempo	
  de	
  cómputo.	
  Por	
  último	
  se	
  explica	
  cómo	
  se	
  clasifica	
  la	
  dinámica	
  
molecular	
   en	
   términos	
   de	
   equilibrio	
   y	
   se	
   muestra	
   un	
   diagrama	
   del	
   programa	
   de	
  
simulación.	
  
1.1 Qué	
  es	
  la	
  dinámica	
  molecular	
  	
  
	
  
La	
  Dinámica	
  Molecular	
  es	
  una	
   técnica	
  computacional	
  que	
  simula	
  el	
   comportamiento	
  
de	
  un	
  sistema	
  clásico	
  de	
  muchas	
  partículas	
  y	
  a	
  partir	
  del	
  cual,	
  se	
  pueden	
  calcular	
  los	
  
promedios	
  en	
   tiempo	
  de	
   las	
  propiedades	
  del	
   sistema,	
  permitiendo	
  una	
  visualización	
  
del	
   movimiento	
   de	
   partículas.	
   En	
   este	
   contexto,	
   la	
   palabra	
   ‘clásico’	
   significa	
   que	
   el	
  
movimiento	
   fundamental	
   de	
   las	
   partículas	
   constituyentes	
   obedece	
   las	
   leyes	
   de	
   la	
  
mecánica	
  clásica.	
  
La	
   dinámica	
   molecular	
   es	
   un	
   método	
   determinístico11,	
   es	
   decir,	
   que	
   el	
   estado	
   del	
  
sistema	
  en	
  un	
  tiempo	
  dado	
  en	
  el	
   futuro	
  se	
  puede	
  predecir	
  a	
  partir	
  del	
  estado	
  actual.	
  
En	
  cada	
  paso,	
  las	
  fuerzas	
  sobre	
  los	
  átomos	
  se	
  calculan	
  y	
  combinan	
  con	
  las	
  posiciones	
  y	
  
velocidades	
   actuales	
   para	
   generar	
   nuevas	
   posiciones	
   y	
   velocidades	
   a	
   plazos	
  
inmediatos.	
  La	
  fuerza	
  que	
  actúa	
  sobre	
  cada	
  átomo	
  se	
  supone	
  que	
  es	
  constante	
  durante	
  
el	
   intervalo	
   de	
   tiempo.	
   Los	
   átomos	
   luego	
   se	
   mueven	
   a	
   sus	
   nuevas	
   posiciones,	
   se	
  
actualiza	
  el	
  conjunto	
  de	
  fuerzas	
  y	
  un	
  nuevo	
  ciclo	
  de	
  Dinámica	
  Molecular	
  comienza2.	
  
	
  	
  	
  	
  	
  	
  13	
  
	
  
Las	
  trayectorias	
  se	
  obtienen	
  al	
  resolver	
  las	
  ecuaciones	
  diferenciales	
  de	
  la	
  segunda	
  ley	
  
de	
  Newton1:	
  
d 2xi
dt2
=
Fxi
mp, i
	
   	
   	
   	
   	
   (1)	
  
La	
   Dinámica	
   Molecular	
   es	
   una	
   excelente	
   aproximación	
   para	
   un	
   amplio	
   rango	
   de	
  
materiales	
  y	
  en	
  varios	
  aspectos	
  es	
  muy	
  similar	
  a	
   la	
  experimentación	
  real.	
  Cuando	
  se	
  
realiza	
  un	
  experimento	
  real,	
  se	
  prepara	
  una	
  muestra	
  del	
  material	
  a	
  estudiar,	
  se	
  coloca	
  
la	
   muestra	
   en	
   el	
   instrumento	
   de	
   medición	
   (como	
   un	
   termómetro,	
   manómetro	
   o	
  
viscosímetro)	
  y	
  se	
  mide	
  la	
  propiedad	
  de	
  interés	
  durante	
  un	
  cierto	
  intervalo	
  de	
  tiempo.	
  
Si	
  las	
  mediciones	
  están	
  sujetas	
  a	
  ruido	
  estadístico	
  (como	
  lo	
  son	
  la	
  mayoría),	
  entonces	
  
mientras	
   más	
   resultados	
   obtenidos	
   se	
   promedien,	
   más	
   precisa	
   se	
   convierte	
   la	
  
medición.	
  En	
  una	
  simulación	
  molecular	
  se	
  sigue	
  exactamente	
  el	
  mismo	
  procedimiento.	
  
Primero	
  se	
  prepara	
  una	
  muestra:	
  se	
  selecciona	
  un	
  sistema	
  modelo	
  que	
  consiste	
  de	
  N	
  
partículas	
   y	
   se	
   resuelven	
   las	
   ecuaciones	
   de	
   Newton	
   de	
   movimiento	
   hasta	
   que	
   las	
  
propiedades	
   del	
   sistema	
   no	
   cambien	
   más	
   con	
   el	
   tiempo	
   (es	
   decir,	
   se	
   calibra	
   el	
  
sistema).	
  Después	
  de	
  la	
  calibración,	
  se	
  realizan	
  las	
  mediciones.	
  	
  De	
  hecho,	
  algunos	
  de	
  
los	
  errores	
  más	
  comunes	
  que	
  se	
  pueden	
  cometer	
  en	
  un	
  experimento	
  por	
  computadora	
  
son	
  muy	
   parecidos	
   a	
   los	
   errores	
   que	
   se	
   pueden	
   hacer	
   en	
   experimentos	
   reales	
   (por	
  
ejemplo,	
  la	
  muestra	
  no	
  está	
  preparada	
  correctamente,	
  las	
  mediciones	
  son	
  muy	
  pocas	
  o	
  
no	
   se	
   está	
   midiendo	
   lo	
   que	
   se	
   piensa).	
   Para	
   medir	
   una	
   cantidad	
   observable	
   en	
  
Dinámica	
   Molecular,	
   se	
   debe	
   primero	
   ser	
   capaz	
   de	
   expresar	
   esta	
   cantidad	
   como	
  
función	
  de	
  las	
  posiciones	
  y	
  momentos	
  de	
  las	
  partículas	
  del	
  sistema.	
  
Un	
   concepto	
   importante	
   en	
   dinámica	
  molecular	
   es	
   el	
   ensamble.	
   Un	
   ensamble	
   es	
   un	
  
conjunto	
   de	
   sistemas	
   termodinámicos	
   que	
   permiten	
   hacer	
   un	
   análisis	
   estadístico.	
  
Existen	
   diversos	
   tipos	
   de	
   ensambles	
   estadísticos,	
   pero	
   los	
   más	
   utilizados	
   son	
   el	
  
microcanónico,	
   el	
   canónico	
   y	
   el	
   gran	
   canónico2,12.	
   El	
   ensamble	
  microcanónico	
   es	
   un	
  
sistema	
  termodinámico	
  donde	
  no	
  se	
  intercambia	
  energía	
  ni	
  materia	
  con	
  el	
  ambiente,	
  
el	
   ensamble	
   canónico	
   es	
   aquel	
   donde	
   se	
   intercambia	
   energía	
   pero	
   no	
   materia	
   y	
   el	
  
macrocanónico	
   donde	
   se	
   intercambian	
   ambas	
   cosas.	
   Para	
   el	
   programa	
   de	
  Dinámica	
  
14	
   	
  
	
  
Molecular	
   que	
   se	
   realizará	
   en	
   este	
   trabajo,	
   será	
   de	
   especial	
   interés	
   el	
   ensamble	
  
canónico.	
   Para	
   este	
   ensamble	
   las	
   variables	
   fijas	
   en	
   el	
   sistema	
   son:	
   el	
   número	
   de	
  
partículas	
  (Np),	
  el	
  volumen	
  (V	
  =	
  L3)	
  y	
  la	
  temperatura	
  (T).	
  Con	
  estas	
  tres	
  variables,	
  se	
  
podrán	
  calcular	
  propiedades	
  fisicoquímicas,	
  interfaciales	
  y	
  reológicas.	
  	
  
Con	
  el	
  objetivode	
  facilitar	
  el	
  manejo	
  de	
  las	
  variables	
  que	
  se	
  ocupan	
  en	
  un	
  programa	
  
de	
   Dinámica	
  Molecular,	
   es	
   conveniente	
   convertirlas	
   en	
  magnitudes	
   adimensionales.	
  
En	
  el	
  Apéndice	
  1	
  se	
  presenta	
  el	
  procedimiento	
  para	
   la	
  conversión	
  de	
  magnitudes	
  de	
  
las	
  principales	
  variables	
  adimensionales13.	
  
1.2 Convención	
  de	
  imagen	
  mínima	
  y	
  condiciones	
  de	
  frontera	
  periódicas	
  	
  
	
  
Se	
   considera	
   como	
   la	
   celda	
   principal	
   o	
   caja	
   central	
   al	
   volumen	
   que	
   contiene	
   las	
   N	
  
partículas	
   y	
   se	
   puede	
   imaginar	
   que	
   ésta	
   se	
   reproduce	
   infinitas	
   veces	
  mediante	
   una	
  
traslación	
   rígida	
   en	
   las	
   tres	
   direcciones	
   cartesianas	
   (imagen	
   mínima),	
   como	
   se	
  
muestra	
  en	
  la	
  Figura	
  12.	
  
	
  
	
  
Figura	
  	
  1.	
  Convención	
  de	
  imagen	
  mínima	
  en	
  un	
  sistema	
  de	
  tres	
  dimensiones	
  
	
  
En	
   una	
   simulación	
   de	
   Dinámica	
   Molecular	
   es	
   frecuente	
   utilizar	
   Condiciones	
   de	
  
Frontera	
  Periódicas	
  (CFP)	
  para	
  eliminar	
  efectos	
  de	
  superficie1.	
  
Celda&
Principal&
Celda&
Principal&
	
  	
  	
  	
  	
  	
  15	
  
	
  
	
  
Figura	
  	
  2.	
  Representación	
  esquemática	
  de	
  las	
  Condiciones	
  de	
  Frontera	
  Periódicas.	
  
	
  
Durante	
  la	
  simulación,	
  cuando	
  una	
  molécula	
  se	
  mueve	
  en	
  la	
  celda	
  principal,	
  su	
  imagen	
  
periódica	
   en	
   cada	
   una	
   de	
   las	
   otras	
   celdas	
   se	
   mueve	
   exactamente	
   con	
   la	
   misma	
  
orientación	
  y	
  de	
  la	
  misma	
  forma.	
  De	
  este	
  modo,	
  conforme	
  una	
  molécula	
  abandona	
  la	
  
caja	
   central,	
   una	
   de	
   sus	
   imágenes	
   entra	
   en	
   la	
   cara	
   opuesta.	
   Es	
   decir,	
   no	
   existen	
  
fronteras	
  o	
  paredes	
  en	
  la	
  celda	
  principal	
  y	
  por	
  lo	
  tanto	
  el	
  sistema	
  no	
  tiene	
  superficie.	
  
En	
  la	
  Figura	
  2	
  se	
  muestra	
  una	
  representación	
  esquemática	
  de	
  las	
  CFP.	
  
Cuando	
   se	
   utilizan	
   las	
   condiciones	
   de	
   frontera	
   periódicas	
   la	
   distancia	
   en	
   cualquier	
  
dirección	
   entre	
   la	
   partícula	
   i	
   y	
   la	
   imagen	
  más	
   cercana	
   de	
   j	
   siempre	
   será	
  menor	
   (en	
  
valor	
  absoluto)	
  a	
   la	
  mitad	
  de	
  la	
  caja	
  (L/2).	
  La	
  manera	
  de	
  computar	
  esta	
  distancia	
  en	
  
fortran	
  es	
  usando	
  la	
  función	
  nint(x)	
  (en	
  inglés,	
  nearest	
  integer	
  function).	
  Esta	
  función	
  
redondea	
  un	
  número	
  real	
  al	
  entero	
  más	
  próximo14.	
  
1.3 Potencial	
  de	
  Lennard-­‐Jones	
  
	
  
Un	
  par	
  de	
  átomos	
  o	
  moléculas	
  neutros	
  están	
  sujetos	
  a	
  dos	
  fuerzas	
  distintas:	
  una	
  fuerza	
  
atractiva	
  que	
  actúa	
  cuando	
  están	
  separados	
  a	
  grandes	
  distancias	
  (fuerza	
  de	
  Van	
  Der	
  
Waals)	
   y	
   una	
   fuerza	
   repulsiva	
   que	
   actúa	
   cuando	
   están	
   separados	
   a	
   pequeñas	
  
distancias	
   (principio	
   de	
   exclusión	
   de	
   Pauli).	
   El	
   potencial	
   de	
   Lennard-­‐Jones	
   es	
   un	
  
modelo	
  matemático	
  sencillo	
  que	
  representa	
  este	
  comportamiento.	
  Fue	
  propuesto	
  en	
  
1924	
  por	
  John	
  Lennard-­‐Jones	
  y	
  es	
  de	
  la	
  forma:	
  
	
  
	
  
	
  	
  
	
  
	
  
	
  
16	
   	
  
	
  
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎠
⎞⎜
⎝
⎛−⎟
⎠
⎞⎜
⎝
⎛=
612
4)(
rr
rU σσε 	
   	
   	
   	
   (2)	
  
donde	
   ε	
   es	
   la	
   profundidad	
   del	
   potencial	
   (con	
   unidades	
   de	
   temperatura	
   según	
   la	
  
relación	
  𝜀/𝑘!),	
  σ	
  es	
  la	
  distancia	
  finita	
  (en	
  unidades	
  de	
  longitud)	
  en	
  la	
  que	
  el	
  potencial	
  
entre	
  partículas	
  es	
  cero	
  y	
  r	
  es	
  la	
  distancia	
  entre	
  partículas.	
  El	
  término	
   !
!
!"
	
  describe	
  la	
  
repulsión	
   entre	
   partículas	
   y	
   el	
   término	
   !
!
!
la	
   atracción	
   entre	
   ellas.	
   Los	
   valores	
   de	
  
𝜀/𝑘! ≈ 120𝐾	
   y	
   𝜎 ≈ 0.34  𝑛𝑚	
   concuerdan	
   razonablemente	
   con	
   las	
   propiedades	
  
experimentales	
  del	
  argón	
  líquido.	
  	
  
2 	
  
Figura	
  	
  3.	
  Potencial	
  de	
  Lennard-­‐Jones	
  entre	
  dos	
  partículas	
  
	
  
En	
  la	
  Figura	
  3	
  se	
  muestra	
  la	
  gráfica	
  del	
  potencial	
  de	
  Lennard-­‐Jones	
  para	
  dos	
  partículas	
  
de	
  Argón.	
  Los	
  valores	
  de	
  los	
  ejes	
  tienen	
  unidades	
  adimensionales.	
  
1.4 Posiciones	
  iniciales	
  	
  
	
  
Para	
  iniciar	
  una	
  simulación,	
  se	
  deben	
  asignar	
  las	
  posiciones	
  y	
  velocidades	
  iniciales	
  de	
  
todas	
   las	
   partículas	
   en	
   el	
   sistema.	
   Las	
   posiciones	
   iniciales	
   de	
   las	
   partículas	
   deben	
  
seleccionarse	
  de	
  tal	
  manera	
  que	
  no	
  se	
  traslapen	
  unas	
  con	
  otras.	
  Con	
  frecuencia,	
  esto	
  
0.5 1.0 1.5 2.0 2.5
r
�1
0
1
2
3
U�r⇥
⇤
⇥
	
  	
  	
  	
  	
  	
  17	
  
	
  
se	
   logra	
   colocando	
   las	
  moléculas	
   inicialmente	
   en	
   las	
   caras	
   de	
   un	
   cubo,	
   y	
   corriendo	
  
muchos	
  pasos	
  de	
  simulación	
  hasta	
  que	
  se	
  logre	
  desaparecer	
  la	
  configuración	
  inicial15.	
  
Sin	
  embargo,	
  en	
  este	
  trabajo	
  se	
  comprobó	
  que	
  se	
  ahorra	
  más	
  tiempo	
  y	
  es	
  aún	
  más	
  real	
  
la	
  simulación	
  si	
  se	
  comienza	
  con	
  una	
  configuración	
  inicial	
  de	
  partículas	
  aleatoria.	
  	
  Esto	
  
se	
  logra	
  colocando	
  las	
  partículas	
  en	
  una	
  caja	
  o	
  celda	
  más	
  grande	
  que	
  la	
  caja	
  objetivo	
  
(alrededor	
   de	
   4	
   veces	
  mayor)	
   y	
   con	
   una	
   condición	
   de	
   distancia	
  mínima	
   entre	
   cada	
  
partícula;	
  como	
  la	
  caja	
  es	
  grande,	
  la	
  probabilidad	
  de	
  traslape	
  es	
  baja.	
  En	
  el	
  Apéndice	
  2	
  
se	
   muestra	
   el	
   algoritmo	
   en	
   fortran	
   para	
   el	
   acomodo	
   inicial	
   de	
   partículas.	
  
Consecuentemente,	
   se	
   comprime	
   la	
   caja	
   isotérmicamente	
   hasta	
   lograr	
   la	
   densidad	
  
deseada.	
  En	
  la	
  sección	
  1.10	
  se	
  hablará	
  del	
  procedimiento	
  de	
  compresión	
  para	
  alcanzar	
  
las	
  dimensiones	
  deseadas.	
  
Para	
   formar	
   cadenas	
   es	
   importante	
   saber	
   qué	
   modelo	
   utilizar.	
   La	
   mayoría	
   de	
   las	
  
simulaciones	
  de	
  polímeros	
  usan	
  modelos	
  de	
  átomos	
  unidos	
  con	
  esferas	
  de	
  Lennard-­‐
Jones16,	
  enlaces	
  estrechos	
  o	
  constreñidos	
  y	
  ángulos	
  de	
  torsión.	
  Pero	
  en	
  este	
  trabajo	
  se	
  
utilizará	
   el	
  modelo	
   de	
   cadenas	
   flexibles	
   con	
   enlace	
   tipo	
   resorte	
   (Figura	
   4).	
   En	
   este	
  
trabajo	
  se	
  define	
  como	
  cadena	
  a	
  la	
  unión	
  de	
  átomos	
  de	
  la	
  misma	
  especie,	
  donde	
  m	
  es	
  
el	
  número	
  de	
  átomos	
  en	
  la	
  cadena.	
  	
  
	
  	
  
	
  
Figura	
  	
  4.	
  Modelo	
  de	
  un	
  enlace	
  tipo	
  resorte	
  entre	
  dos	
  esferas	
  de	
  Lennard-­‐	
  Jones.	
  	
  
	
  
Los	
   átomos	
   en	
   cada	
   cadena	
   se	
   modelan	
   con	
   el	
   potencial	
   de	
   Lennard-­‐Jones	
   y	
   el	
  
potencial	
  de	
  enlace	
  entre	
  ellos	
  se	
  modela	
  con	
  resortes	
  armónicos	
  aplicando	
  la	
  Ley	
  de	
  
Hooke.	
  Como	
  las	
  cadenas	
  son	
  flexibles	
  no	
  hay	
  necesidad	
  de	
  calcular	
  ángulos	
  de	
  enlace	
  
ni	
  de	
  torsión5.	
  
Elprimer	
  paso	
  es	
  saber	
  generar	
  un	
  enlace	
  con	
  dirección	
  aleatoria17:	
  
r"
18	
   	
  
	
  
1. Primero	
   se	
   deben	
   generar	
   los	
   valores	
   aleatorios	
   V1	
   y	
   V2	
   tomados	
   de	
   una	
  
distribución	
  uniforme	
  de	
  (-­‐1,	
  1)	
  de	
  tal	
  manera	
  que:	
  S=(V12	
  +	
  V22)	
  <1.	
  
2. El	
  vector	
  aleatorio	
  es	
  entonces:	
   2V1 1− S, 2V2 1− S, 1− 2S( ) 	
  
Respecto	
  al	
   tamaño	
  del	
  vector,	
  éste	
  se	
  considera	
  arbitrario	
  en	
  un	
  valor	
   igual	
  a	
  σ.	
  De	
  
este	
  modo,	
   se	
   pueden	
  unir	
   cuantos	
   átomos	
   se	
   deseé	
   para	
   formar	
   una	
   cadena.	
   En	
   la	
  
Figura	
  5	
  se	
  muestra	
  un	
  ejemplo	
  de	
  posiciones	
  aleatorias	
  para	
  100	
  cadenas	
  formadas	
  
por	
  8	
  átomos	
  cada	
  una.	
  En	
  el	
  Apéndice	
  3	
  se	
  muestra	
  el	
  algoritmo	
  en	
   fortran	
  para	
   la	
  
formación	
  de	
  cadenas.	
  
	
  
Figura	
  	
  5.	
  Posiciones	
  aleatorias	
  para	
  una	
  cadena	
  de	
  8	
  átomos.	
  
1.5 Radio	
  de	
  corte	
  y	
  lista	
  de	
  vecinos	
  
	
  
El	
   cálculo	
   de	
   fuerzas	
   que	
   actúa	
   sobre	
   cada	
   partícula	
   es	
   la	
   parte	
   que	
   más	
   tiempo	
  
consume	
   en	
   la	
  mayoría	
   de	
   las	
   simulaciones	
   de	
  Dinámica	
  Molecular,	
   ya	
   que	
   se	
   tiene	
  
que	
  considerar	
   la	
  contribución	
  a	
   la	
   fuerza	
  de	
   todas	
   las	
  partículas	
  que	
  se	
  encuentran	
  
dentro	
  de	
   la	
  caja	
  sobre	
   la	
  partícula	
   i,	
  y	
  esto	
  mismo	
  para	
   todas	
   las	
  partículas.	
  Si	
  esto	
  
fuera	
   de	
   este	
   modo,	
   implicaría	
   que	
   el	
   tiempo	
   necesario	
   para	
   la	
   evaluación	
   de	
   las	
  
fuerzas	
  sería	
  del	
  orden	
  de	
  N	
  x	
  N	
  (N2)2.	
  Para	
  el	
  ahorro	
  de	
  tiempo	
  de	
  cómputo,	
  se	
  han	
  
desarrollado	
  varios	
  “trucos”	
  que	
  permiten	
  disminuir	
  el	
  orden	
  de	
  N2	
  a	
  simplemente	
  N.	
  	
  
Uno	
   de	
   los	
  métodos	
   para	
   ahorrar	
   tiempo	
   computacional,	
   es	
   truncar	
   el	
   potencial	
   de	
  
Lennard-­‐Jones	
  a	
  una	
  distancia	
  límite	
  que	
  se	
  denomina	
  radio	
  de	
  corte.	
  Por	
  lo	
  general	
  el	
  
-200
20
-20 0 20
-20
-10
0
10
20
	
  	
  	
  	
  	
  	
  19	
  
	
  
radio	
  de	
  corte	
  rc	
  =	
  2.5	
  σ.	
  Esto	
  significa	
  que	
  todas	
  las	
  partículas	
  que	
  se	
  encuentren	
  fuera	
  
de	
  este	
  radio	
  de	
  corte	
  no	
  contribuyen	
  a	
  la	
  fuerza	
  sobre	
  la	
  partícula	
  i:	
  	
  
Si r ≤ rc U(r) = 4ε
σ
r
"
#
$
%
&
'
12
−
σ
r
"
#
$
%
&
'
6)
*
+
+
,
-
.
.
Si r > rc 0
	
   	
   	
   	
   (9) 	
  
La	
   lista	
   de	
   vecinos	
   es	
   otro	
   ingenioso	
   truco	
   para	
   ahorrar	
   tiempo	
   de	
   cómputo.	
   Cada	
  
átomo	
  en	
  el	
  sistema	
  tiene	
  alrededor	
  una	
  cantidad	
  de	
  vecinos	
  que	
  son	
  potencialmente	
  
probables	
  de	
  estar	
  dentro	
  del	
  radio	
  de	
  corte16.	
  De	
  este	
  modo,	
  se	
  investiga	
  cuáles	
  son	
  
estos	
   vecinos	
   y	
   se	
   calculan	
   únicamente	
   las	
   fuerzas	
   entre	
   ellos,	
   ya	
   no	
   con	
   todos	
   los	
  
átomos,	
   ahorrándose	
   así	
   hasta	
   un	
   80%	
   del	
   cálculo.	
   La	
   lista	
   de	
   vecinos	
   se	
   debe	
  
actualizar	
  cada	
  cierto	
  número	
  de	
  pasos,	
  dependiendo	
  de	
  la	
  precisión	
  deseada.	
  	
  El	
  radio	
  
de	
   la	
   lista	
  de	
  vecinos	
  en	
  este	
   trabajo	
  es	
   rv	
  =	
  3.6	
  σ.	
  El	
   radio	
  de	
   la	
   lista	
  de	
  vecinos	
  así	
  
como	
  el	
  radio	
  de	
  corte	
  se	
  muestran	
  ejemplificados	
  en	
  la	
  Figura	
  6.	
  En	
  el	
  Apéndice	
  4	
  se	
  
muestra	
  el	
  algoritmo	
  en	
  fortran	
  para	
  la	
  formación	
  de	
  la	
  lista	
  de	
  vecinos.	
  
	
  
Figura	
  	
  6.	
  Ejemplo	
  de	
  radio	
  de	
  corte	
  y	
  radio	
  de	
  la	
  lista	
  de	
  vecinos.	
  
1.6 Velocidades	
  iniciales	
  y	
  temperatura	
  
	
  
Respecto	
  a	
  las	
  velocidades	
  iniciales,	
  éstas	
  también	
  son	
  aleatorias	
  y	
  valen	
  de	
  -­‐0.5	
  a	
  0.5.	
  
Sin	
   embargo	
   se	
   debe	
   realizar	
   un	
   escalamiento	
   de	
   tal	
   forma	
   que	
   el	
   valor	
   de	
   las	
  
velocidades	
  resultantes	
  se	
  ajuste	
  al	
  promedio	
  de	
  energía	
  cinética	
  del	
  valor	
  deseado	
  y,	
  
por	
  equilibrio	
  térmico	
  se	
  cumpla	
  la	
  siguiente	
  relación:	
  
rc#
rv#
20	
   	
  
	
  
	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
   vx
2 + vy
2 + vz
2 =
kBT
mp
	
  	
   	
  	
  	
  	
  	
  	
  	
   	
  	
   	
  	
  	
  	
   (3)	
  
donde	
  v	
  es	
  la	
  velocidad	
  en	
  la	
  componente	
  x,	
  y	
  o	
  z	
  de	
  una	
  partícula,	
  kB	
  es	
  la	
  constante	
  
de	
  Boltzmann,	
   T	
   es	
   la	
   temperatura	
   y	
  mp	
   es	
   la	
  masa	
   de	
   la	
   partícula.	
   Si	
   se	
   despeja	
   la	
  
temperatura	
  y	
  se	
  cuenta	
  para	
  todas	
  las	
  partículas,	
  se	
  tiene	
  que:	
  
	
   	
   T0 =
mp
3NpkB
(vx,i
2 + vy,i
2 + vz,i
2 )
i=1
Np
∑ 	
  	
   	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
   	
   (4)	
  
donde	
  el	
  tres	
  se	
  debe	
  a	
  las	
  tres	
  dimensiones	
  utilizadas.	
  De	
  esta	
  forma	
  se	
  puede	
  ajustar	
  
la	
   temperatura	
   instantánea	
  T0	
   para	
   igualar	
   la	
   temperatura	
  deseada	
   escalando	
   todas	
  
las	
  velocidades	
  con	
  la	
  ecuación:
	
  	
  
0
s T
Tf =
	
   	
   	
   	
   	
   	
  
(5) 	
  
	
  La	
   velocidad	
   del	
   centro	
   de	
   masa	
   deberá	
   ser	
   igual	
   a	
   cero.	
   Cabe	
   mencionar	
   que	
   la	
  
mayoría	
  de	
  las	
  personas	
  que	
  trabajan	
  con	
  dinámica	
  molecular	
  acostumbran	
  a	
  utilizar	
  
termostatos	
  como	
  Berendsen	
  o	
  Nosé-­‐Hoover18,	
  pero	
  en	
  este	
  trabajo	
  el	
  termostato	
  que	
  
regula	
   la	
   velocidad	
   es	
   el	
   escalamiento	
   simple	
   de	
   velocidades.	
   En	
   el	
   Apéndice	
   5	
   se	
  
muestra	
  el	
  algoritmo	
  en	
  fortran	
  para	
  el	
  escalamiento	
  de	
  velocidades.	
  
1.7 Fuerzas	
  intermoleculares	
  y	
  fuerzas	
  de	
  enlace	
  
	
  
Para	
  el	
   cálculo	
  de	
   fuerzas	
   se	
   tienen	
  que	
   tomar	
  en	
  cuenta	
   las	
   fuerzas	
  de	
  enlace	
  y	
   las	
  
fuerzas	
   intermoleculares.	
   Para	
   las	
   fuerzas	
   intermoleculares,	
   el	
   mismo	
   modelo	
   de	
  
potencial	
  de	
  Lennard-­‐	
  Jones	
  es	
  utilizado.	
  La	
  fuerza	
  a	
  la	
  que	
  están	
  sujetas	
  las	
  partículas	
  
es	
  el	
  negativo	
  del	
  gradiente	
  del	
  potencial	
  (ver	
  Ecuación	
  6):	
  
F(r) = −∇U(r) = − dU(r)
dr
r = 4ε 12σ
12
r13
− 6σ
6
r7
#
$
%
&
'
(r =
48ε
r
σ
r
)
*
+
,
-
.
12
−
1
2
σ
r
)
*
+
,
-
.
6#
$
%
%
&
'
(
(
r 	
  	
   (6)	
  	
  	
  	
  	
  	
  	
  
	
  	
  	
  	
  	
  	
  21	
  
	
  
Para	
  las	
  fuerzas	
  de	
  enlace,	
  se	
  sigue	
  el	
  modelo	
  de	
  enlace	
  tipo	
  resorte	
  que	
  obedece	
  a	
  la	
  
Ley	
  de	
  Hooke:	
  
U enlace = 1
2
k r − r0( )
2 	
   	
   	
   	
   	
   (7)	
  	
  	
  	
  	
  	
  	
  	
  	
  
F enlace = k r − r0( ) 	
  	
   	
   	
   	
   	
   (8)	
  	
  	
  	
  	
  	
  	
  
donde	
   k	
   es	
   igual	
   a	
   3000	
   ε/σ2	
   	
   y	
   r0	
   es	
   igual	
   a	
   1σ 5.	
   La	
   constante	
   del	
   resorte	
   en	
   el	
  
potencial	
  armónico	
  es	
  igual	
  a	
  este	
  valor	
  porque	
  esel	
  límite	
  donde	
  se	
  asegura	
  que	
  las	
  
cadenas	
  están	
  efectivamente	
  unidas	
  con	
  un	
  enlace	
  de	
  longitud	
  σ;	
  arriba	
  de	
  este	
  valor,	
  
los	
  resultados	
  son	
  idénticos5.	
  
	
  
	
  
Figura	
  	
  7.	
  Fuerzas	
  de	
  enlace	
  y	
  fuerzas	
  intermoleculares.	
  
	
  
En	
  la	
  Figura	
  7	
  se	
  distinguen	
  las	
  fuerzas	
  de	
  enlace,	
  que	
  se	
  encuentran	
  entre	
  dos	
  átomos	
  
unidos	
  y	
   las	
   fuerzas	
   intermoleculares	
   localizadas	
   al	
   interactuar	
   átomos	
  de	
  diferente	
  
cadena	
  lineal.	
  En	
  el	
  Apéndice	
  6	
  se	
  muestra	
  el	
  algoritmo	
  en	
  fortran	
  para	
  el	
  cálculo	
  de	
  
fuerzas	
  de	
  enlace	
  y	
  el	
  cálculo	
  de	
  fuerzas	
  intermoleculares.	
  
1.8 Diseño	
  de	
  un	
  programa	
  de	
  simulación	
  molecular	
  
	
  
Para	
   diseñar	
   una	
   simulación	
   de	
   dinámica	
   molecular	
   se	
   debe	
   tomar	
   en	
   cuenta	
   la	
  
capacidad	
  computacional	
  disponible,	
  el	
  tamaño	
  de	
  partículas,	
  el	
  tamaño	
  de	
  paso	
  y	
  el	
  
tiempo	
   total	
   de	
   duración,	
   de	
   tal	
   forma	
   que	
   el	
   cálculo	
   pueda	
   terminar	
   dentro	
   de	
   un	
  
período	
  razonable.	
  	
  Sin	
  embargo,	
  las	
  simulaciones	
  deben	
  durar	
  lo	
  suficiente	
  para	
  que	
  
sean	
   de	
   interés	
   para	
   las	
   escalas	
   de	
   tiempo	
   de	
   los	
   procesos	
   naturales	
   que	
   se	
   están	
  
estudiando.	
  Para	
  tener	
  conclusiones	
  estadísticas	
  válidas	
  de	
  las	
  simulaciones,	
  el	
   lapso	
  
Fuerza'de'enlace'
Fuerza'intermolecular'
(Lennard2Jones)'
Fuerza'de'enlace'
22	
   	
  
	
  
de	
   tiempo	
   simulado	
   debe	
   coincidir	
   con	
   la	
   cinética	
   del	
   proceso	
   natural.	
   Además	
   el	
  
tamaño	
   de	
   paso	
   se	
   elige	
   lo	
   suficientemente	
   pequeño	
   como	
   para	
   evitar	
   errores	
   de	
  
discretización,	
   es	
   decir,	
   menor	
   que	
   la	
   frecuencia	
   de	
   vibración	
   más	
   rápida	
   en	
   el	
  
sistema.	
  Los	
  tamaños	
  de	
  paso	
  típicos	
  para	
  dinámica	
  molecular	
  clásica	
  son	
  del	
  orden	
  de	
  
1	
  femtosegundo	
  (10-­‐15	
  s),	
  valor	
  utilizado	
  en	
  este	
  trabajo.	
  
El	
   número	
   de	
   moléculas	
   se	
   eligió	
   de	
   la	
   siguiente	
   manera:	
   1000	
   para	
   cadenas	
  
constituidas	
  por	
  un	
   solo	
  átomo,	
  800	
  para	
  m=2,	
  400	
  para	
  m=4,	
  300	
  para	
  m=8	
  y	
  200	
  
para	
  m=16.	
  El	
  número	
  total	
  de	
  pasos	
  varía	
  según	
  la	
  propiedad	
  que	
  se	
  requiere	
  calcular	
  
y	
   el	
   tiempo	
   necesario	
   para	
   alcanzar	
   la	
   condición	
   de	
   equilibrio	
   en	
   general	
   va	
   de	
   los	
  
500,000	
  pasos	
  a	
  los	
  2	
  millones	
  de	
  pasos,	
  según	
  sea	
  el	
  caso.	
  	
  
1.9 Ecuaciones	
  de	
  movimiento	
  
	
  
Una	
  vez	
  que	
  se	
  han	
  calculado	
  las	
  fuerzas	
  intermoleculares	
  y	
  de	
  enlace	
  se	
  suman	
  y	
  se	
  
calcula	
   la	
   aceleración	
   con	
   la	
   segunda	
   ley	
   de	
   Newton	
   (Ecuación	
   1).	
   Con	
   esta	
  
aceleración,	
  con	
  las	
  posiciones	
  y	
  las	
  velocidades	
  se	
  pueden	
  integrar	
  las	
  ecuaciones	
  de	
  
movimiento	
   de	
   Newton.	
   El	
   algoritmo	
   propuesto	
   en	
   este	
   trabajo	
   es	
   el	
   llamado	
  
algoritmo	
   de	
   Verlet–velocidad19.	
   Este	
   algoritmo	
   almacena	
   posiciones,	
   velocidades	
   y	
  
aceleraciones	
  todas	
  al	
  mismo	
  tiempo	
  t	
  y	
  también	
  minimiza	
  errores	
  de	
  redondeo	
  que	
  
son	
   mayores	
   en	
   otros	
   algoritmos	
   como	
   Verlet	
   simple	
   y	
   Leapfrog.	
   	
   La	
   forma	
   del	
  
algoritmo	
  de	
  Verlet–velocidad	
  es:	
  
x (t +Δt) = x(t) + v(t) Δt + 1
2
a(t) Δt2
v (t +Δt) = v(t) + 1
2
a(t) + a(t +Δt)( ) Δt
	
   	
   	
   (10)	
  
Donde	
  a	
  es	
  la	
  aceleración	
  que	
  se	
  calcula	
  con	
  la	
  Ecuación	
  1,	
  y	
  la	
  Fuerza	
  es	
  la	
  obtenida	
  
de	
   sumar	
   las	
   fuerzas	
   intermoleculares	
   y	
   las	
   fuerzas	
   de	
   enlace.	
   En	
   el	
   Apéndice	
   7	
   se	
  
muestra	
  el	
  algoritmo	
  en	
  fortran	
  para	
  la	
  solución	
  de	
  ecuaciones	
  de	
  movimiento.	
  
	
  	
  	
  	
  	
  	
  23	
  
	
  
1.10 Compresión	
  isotérmica	
  
	
  
En	
  la	
  sección	
  1.4	
  se	
  habló	
  del	
  posicionamiento	
  de	
  las	
  moléculas	
  en	
  una	
  caja	
  mayor	
  a	
  
las	
  dimensiones	
  deseadas	
  con	
   la	
   finalidad	
  de	
  un	
  rápido	
  y	
   fácil	
  acomodo.	
  Para	
   lograr	
  
una	
   densidad	
   específica	
   se	
   debe	
   realizar	
   una	
   compresión	
   isotérmica	
   por	
   pasos.	
   Se	
  
calcula	
  un	
  factor	
  de	
  compresión	
  que	
  relaciona	
  la	
  longitud	
  deseada	
  l	
  entre	
  la	
  longitud	
  
inicial	
  l0	
  y	
  el	
  número	
  de	
  pasos	
  fijos	
  n	
  para	
  llegar	
  a	
  la	
  longitud	
  deseada:	
  
n
l
lf
1
0
⎟⎟⎠
⎞
⎜⎜⎝
⎛
=
	
   	
   	
   	
   	
  
(11)	
  
	
  Este	
   factor	
  se	
  multiplica	
   tanto	
  a	
   las	
   longitudes	
  de	
   la	
  caja	
  como	
   las	
  posiciones	
  de	
   las	
  
partículas,	
   y	
   en	
   cada	
   paso	
   es	
   importante	
   actualizar	
   las	
   fuerzas.	
   En	
   el	
   Apéndice	
   8	
   se	
  
muestra	
  el	
  algoritmo	
  en	
  fortran	
  para	
  la	
  compresión	
  isotérmica.	
  
1.11 Integración	
  completa	
  	
  	
  
	
  
Una	
   vez	
   teniendo	
   las	
   partículas	
   dentro	
   de	
   una	
   caja	
   de	
   simulación	
   con	
   dimensiones	
  
deseadas	
   se	
   procede	
   a	
   estabilizar	
   el	
   sistema	
   con	
   un	
   cierto	
   número	
   de	
   pasos.	
   En	
   la	
  
Figura	
  8	
  se	
  muestran	
  las	
  energías	
  cinética,	
  potencial	
  y	
  total	
  de	
  una	
  simulación	
  sencilla	
  
de	
  dinámica	
  molecular.	
  Se	
  puede	
  observar	
  que	
   la	
  energía	
  total	
  es	
  constante	
  después	
  
de	
   muy	
   pocos	
   pasos	
   de	
   simulación.	
   La	
   temperatura	
   es	
   constante	
   a	
   lo	
   largo	
   de	
   la	
  
simulación	
  de	
  dinámica	
  molecular,	
   lo	
  cual	
   indica	
  que	
  el	
  escalamiento	
  de	
  velocidades	
  
funciona	
  como	
  termostato	
  para	
  mantener	
  la	
  temperatura	
  constante..	
  
24	
   	
  
	
  
	
  
Figura	
  	
  8.	
  Energías	
  resultantes	
  de	
  la	
  integración	
  de	
  una	
  simulación	
  de	
  dinámica	
  molecular	
  
	
  
Una	
   vez	
   calculadas	
   las	
   propiedades	
   de	
   interés	
   es	
   importante	
   calcular	
   también	
   la	
  
incertidumbre	
   o	
   error	
   de	
   cálculo.	
   Para	
   todas	
   las	
   propiedades	
   calculadas	
   en	
   este	
  
trabajo	
   se	
   hicieron	
   los	
   cálculos	
   de	
   la	
   siguiente	
   forma:	
   primero	
   se	
   realiza	
   toda	
   la	
  
prueba	
   y	
   se	
   obtienen	
   tantos	
   valores	
   como	
   sea	
   posible	
   de	
   la	
   misma	
   propiedad.	
   Se	
  
realiza	
  un	
  promedio	
  de	
  todos	
  los	
  valores	
  obtenidos	
  y	
  ese	
  es	
  el	
  valor	
  de	
  la	
  propiedad.	
  
Al	
  mismo	
  tiempo	
  se	
  realiza	
  una	
  desviación	
  estándar	
  de	
  todos	
  los	
  valores	
  obtenidos	
  y	
  
ese	
  es	
  el	
  valor	
  de	
  su	
  incertidumbre.	
  
En	
   la	
   Figura	
   9	
   se	
   muestra	
   un	
   diagrama	
   de	
   una	
   integración	
   completa	
   de	
   dinámica	
  
molecular.	
  
-­‐10000	
  
-­‐5000	
  
0	
  
5000	
  
10000	
  
15000	
  
20000	
  
25000	
  
30000	
  
35000	
  
0	
   5000	
   10000	
   15000	
   20000	
   25000	
  
E	
  
Npasos	
  
Ek	
  
Ep	
  
Et25	
  
	
  
	
  
Figura	
  	
  9.	
  Esquema	
  de	
  la	
  integración	
  completa	
  de	
  un	
  programa	
  de	
  dinámica	
  molecular.	
  
DINÁMICA(MOLECULAR(
Posiciones(iniciales(
aleatorias(en(una(caja(((
4(veces(más(grande(
Formación(de(cadenas.(
Tamaño(de(enlace(fijo,(
vector(aleatorio(
Lista(de(vecinos(
Velocidades(iniciales(
aleatorias(en(una(caja(((
4(veces(más(grande(
Escalamiento(de(
velocidades(para(
alcanzar(T(deseada(
Cálculo(de(fuerzas(
intermoleculares(y(de(
enlace(
Ecuaciones(de(
movimiento(con(x0,(v0(
y(a(calculada(
Compresión(isotérmica(
por(pasos(para(llegar(a(
la(densidad(deseada(
Estabilización(del(
sistema(
Cálculo(de(propiedades(
Variables(fijas:(
(Np,(V,(T(
26	
   	
  
	
  
1.12 Clasificación	
  de	
  la	
  dinámica	
  molecular	
  
	
  
En	
  los	
  últimos	
  años	
   la	
  aplicación	
  de	
   la	
  dinámica	
  molecular	
  a	
  problemas	
   importantes	
  
en	
   ciencia	
  de	
   los	
  materiales	
   en	
  general	
   y	
   reología	
   en	
  particular	
   se	
  ha	
   convertido	
  en	
  
una	
  herramienta	
  de	
  análisis	
  indispensable.	
  La	
  dinámica	
  molecular	
  se	
  puede	
  clasificar	
  
en	
   dinámica	
   molecular	
   de	
   equilibrio	
   y	
   dinámica	
   molecular	
   de	
   no	
   equilibrio.	
   La	
  
primera	
   como	
   su	
   nombre	
   lo	
   indica	
   está	
   basada	
   en	
   sistemas	
   termodinámicos	
   en	
  
equilibrio	
   y	
   por	
   no	
   equilibrio	
   se	
   entiende	
   un	
   fluido	
   sobre	
   el	
   cual	
   actúa	
   un	
   campo	
  
externo	
  que	
  lo	
  lleva	
  fuera	
  del	
  equilibrio.	
  	
  
Desde	
   la	
   pasada	
   década,	
   la	
   técnica	
   de	
   dinámica	
   molecular	
   de	
   no	
   equilibrio	
   se	
   ha	
  
convertido	
   en	
   una	
   poderosa	
   herramienta	
   para	
   el	
   estudio	
   de	
   propiedades	
   de	
  
transporte	
   de	
   fluidos	
   simples	
   y	
   moleculares20,21,22.	
   Las	
   investigaciones	
   más	
  
prometedoras	
   incluyen	
   el	
   cálculo	
   de	
   propiedades	
   de	
   transporte	
   como	
   la	
   viscosidad,	
  
conductividad	
  térmica	
  y	
  difusión.	
  
Para	
  el	
  cálculo	
  de	
  la	
  viscosidad,	
  que	
  es	
  el	
  interés	
  de	
  este	
  trabajo,	
  se	
  había	
  utilizado	
  la	
  
dinámica	
  molecular	
   en	
   equilibrio	
   usando	
   las	
   relaciones	
   de	
   Einstein	
   o	
   Green-­‐Kubo1.	
  
Pero	
   estas	
   ecuaciones	
   podían	
   calcular	
   únicamente	
   viscosidades	
   en	
   el	
   límite	
   de	
   cero	
  
rapidez	
  de	
  corte	
  (primera	
  región	
  newtoniana).	
  Si	
  se	
  desea	
  calcular	
  las	
  viscosidades	
  en	
  
el	
   régimen	
   no	
   newtoniano	
   entonces	
   el	
   método	
   de	
   dinámica	
   de	
   no	
   equilibrio	
   debe	
  
utilizarse.	
  	
  
 
A	
  diferencia	
  de	
   los	
  fluidos	
  newtonianos,	
   los	
  polímeros	
  y	
  otros	
  tipos	
  de	
  materiales	
  se	
  
comportan	
   de	
   una	
   manera	
   no-­‐newtoniana	
   y	
   carecen	
   de	
   una	
   ecuación	
   constitutiva	
  
universalmente	
  aceptada	
  para	
  modelar	
  su	
  comportamiento.	
  La	
  dinámica	
  molecular	
  de	
  
no	
   equilibrio	
   se	
   utilizó	
   inicialmente	
   para	
   estudiar	
   la	
   reología	
   de	
   esferas	
   suaves	
   de	
  
Lennard-­‐Jones,	
   sin	
   embargo	
   moléculas	
   más	
   complejas,	
   tales	
   como	
   moléculas	
   de	
  
cadenas	
  cortas	
  y	
  ramificadas	
  son	
  ahora	
  ampliamente	
  simuladas	
  con	
  este	
  método.	
  La	
  
dinámica	
   molecular	
   de	
   no	
   equilibrio	
   tiene	
   por	
   objeto	
   la	
   investigación	
   de	
   las	
  
propiedades	
  de	
   los	
  sistemas	
  de	
   fluidos	
   fuera	
  del	
  equilibrio	
  bajo	
   la	
  acción	
  de	
  campos	
  
externos.	
  
	
  	
  	
  	
  	
  	
  27	
  
	
  
Por	
   lo	
   tanto	
   se	
   puede	
   explicar	
   en	
   el	
   siguiente	
   diagrama	
   las	
   propiedades	
   que	
   se	
  
calcularán	
  y	
  con	
  qué	
  método.	
  	
  
	
  
Figura	
  	
  10.	
  Clasificación	
  de	
  la	
  dinámica	
  molecular	
  y	
  propiedades	
  que	
  se	
  pueden	
  calcular.	
  	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
	
  
28	
   	
  
	
  
CAPÍTULO	
  2	
  
	
  
CÁLCULO	
  DE	
  PROPIEDADES	
  FISICOQUÍMICAS	
  E	
  INTERFACIALES	
  DE	
  CADENAS	
  
LINEALES	
  DE	
  ÁTOMOS	
  
	
  
En	
   este	
   capítulo	
   se	
   dará	
   una	
   explicación	
   termodinámica	
   de	
   equilibrio,	
   estabilidad	
   y	
  
transiciones	
   de	
   fase.	
   Se	
   mostrarán	
   los	
   cálculos	
   y	
   resultados	
   de	
   propiedades	
  
fisicoquímicas	
   para	
   átomos	
   de	
   argón	
   y	
   comparaciones	
   con	
   ecuaciones	
   de	
   estado	
   y	
  
referencias	
   de	
   la	
   literatura.	
   También	
   se	
   mostrarán	
   los	
   resultados	
   de	
   propiedades	
  
fisicoquímicas	
  e	
  interfaciales	
  para	
  cadenas	
  lineales	
  de	
  átomos	
  de	
  la	
  misma	
  especie.	
  	
  
2.1 Equilibrio	
  y	
  estabilidad	
  de	
  los	
  sistemas	
  termodinámicos	
  
	
  
Un	
   sistema	
   se	
   encuentra	
   en	
   equilibrio	
   termodinámico	
   si	
   todas	
   las	
   propiedades	
  
macroscópicas	
  se	
  mantienen	
  sin	
  cambio	
  al	
  pasar	
  el	
   tiempo.	
   	
  Para	
  que	
  un	
  sistema	
  se	
  
encuentre	
  en	
  equilibrio	
  termodinámico	
  se	
  emplea	
  la	
  segunda	
  ley	
  de	
  la	
  termodinámica.	
  
Esta	
   ley	
   establece	
   que	
   cualquier	
   proceso	
   espontáneo	
   debe	
   ir	
   acompañado	
   de	
   un	
  
aumento	
  de	
  la	
  entropía	
  total,	
  esto	
  implica	
  una	
  condición	
  para	
  el	
  equilibrio	
  ya	
  que	
  un	
  
sistema	
  en	
  equilibrio	
  no	
  puede	
  experimentar	
  un	
  cambio	
  espontáneo,	
  dS=0.	
  Se	
  dice	
  que	
  
un	
  sistema	
  heterogéneo	
  compuesto	
  por	
  p	
  fases	
  y	
  m	
  componentes	
  está	
  en	
  equilibrio	
  si	
  
se	
  cumple:	
  
a) Equilibrio	
  térmico	
   	
   T(1)	
  =T(2)	
  =…….=	
  T(p)	
  
b) Equilibrio	
  mecánico	
  	
  	
   P(1)	
  =P(2)	
  =…….=	
  P(p)	
  
c) Equilibrio	
  químico	
  	
   	
   µ1(1)	
  =µ1(2)	
  =…….=	
  µ1(p)	
  
µ2(1)	
  =µ2(2)	
  =…….=	
  µ2(p)	
  
	
  	
  	
  	
  .	
   .	
   	
   .	
  
µm(1)	
  =µm(2)	
  =…….=	
  µm(p)	
  
donde	
  T	
  es	
  la	
  temperatura,	
  P	
  la	
  presión	
  y	
  µ	
  es	
  el	
  potencial	
  químico.	
  	
  
	
  	
  	
  	
  	
  	
  29	
  
	
  
En	
   este	
   trabajo	
   se	
   tiene	
   sólo	
   un	
   componente	
   y	
   un	
   sistema	
   de	
   dos	
   fases	
   (líquido	
   –	
  
vapor).	
  Para	
  que	
  exista	
  un	
  equilibrio	
  termodinámico	
  en	
  este	
  sistema	
  se	
  deben	
  cumplir	
  
entonces	
  las	
  siguientes	
  tres	
  condiciones:	
  
	
  
	
   	
   (12)	
  
	
  
	
  
Con	
   estos	
   criterios	
   se	
   representa	
   el	
   equilibrio	
   térmico,	
   mecánico	
   y	
   químico	
   de	
   la	
  
especie	
   en	
   dos	
   fases,	
   pero	
   esto	
   no	
   garantiza	
   que	
   la	
   fase	
   en	
   equilibrio	
   sea	
   estable.	
  
Adicionalmente	
  a	
   estos	
   criterios	
  para	
  que	
  un	
   sistema	
   sea	
  estable	
   se	
   requiere	
  que	
   la	
  
entropía	
   esté	
   a	
   su	
   máximo,	
   d2S<0.	
   Para	
   que	
   un	
   sistema	
   permanezca	
   en	
   una	
   fase	
  
(estado	
   homogéneo	
   de	
   equilibrio)	
   se	
   deben	
   cumplir	
   criterios	
   de	
   estabilidad.	
   El	
  
primero	
  es	
  la	
  estabilidad	
  térmica:	
  
∂2U
dS2
=
∂T
∂S
≥ 0 ó Cp = T ∂S
∂T p,ni
> 0 	
   	
   	
   (13)	
  
	
  
la	
  estabilidad	
  	
  mecánica:	
  
	
  
∂2U
∂S2
= −
∂P
∂V
≥ 0 ó κ = − 1
V
∂V
∂p T,ni
> 0 	
   	
   (14)	
  
	
  
y	
  la	
  estabilidad	
  química:	
  	
  
0
jnp,T,
>
∂
∂
i
i
n
µ
	
   	
   	
   	
   	
  	
  	
  	
  	
  	
  (15)	
  
Sise	
   satisfacen	
   las	
   Ecuaciones	
   13,	
   14	
   y	
   15	
   la	
   homogeneidad	
   del	
   sistema	
   es	
   estable	
  
frente	
   a	
   perturbaciones	
   internas	
   o	
   externas.	
   	
   Si	
   no	
   se	
   satisfacen	
   los	
   criterios	
   de	
  
T(liq)	
  =	
  T(vap)	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
   
P(liq)	
  =	
  P(vap) 
µ
(liq)	
  =	
  µ(vap) 
30	
   	
  
	
  
estabilidad,	
   el	
   sistema	
   se	
   separa	
   en	
   dos	
   o	
   más	
   partes.	
   Esta	
   separación	
   recibe	
   el	
  
nombre	
  de	
  transiciones	
  de	
  fase23.	
  Existen	
  múltiples	
  diagramas	
  que	
  representan	
  bien	
  
el	
   equilibrio	
   y	
   la	
   transición	
  de	
   fases,	
   y	
   con	
   cada	
  uno	
   se	
  puede	
  visualizar	
  de	
  manera	
  
diferente	
   el	
   mismo	
   problema.	
   Estos	
   diagramas	
   de	
   fases	
   incluyen	
   por	
   ejemplo	
  
diagramas	
   P-­‐T,	
   T-­‐V,	
   P-­‐V,	
   µ-­‐P,	
   etc.	
   A	
   continuación	
   se	
   explicará	
   cada	
   uno	
   con	
   más	
  
detalle.	
  	
  
El	
   primero	
   es	
   el	
   diagrama	
   P-­‐T	
   (Figura	
   11),	
   en	
   el	
   cual	
   se	
   observa	
   una	
   curva	
   de	
  
separación	
  entre	
  el	
  estado	
  líquido	
  y	
  el	
  estado	
  vapor	
  o	
  gas.	
  Esta	
  curva	
  se	
  conoce	
  como	
  
curva	
   de	
   saturación	
   y	
   es	
   el	
   punto	
   exacto	
   en	
   donde	
   ocurre	
   la	
   transición	
   de	
   fase.	
   La	
  
terminación	
  de	
   la	
   curva	
  de	
   fases	
   recibe	
   el	
   nombre	
  de	
  punto	
   crítico.	
   La	
   temperatura	
  
crítica	
  y	
   la	
  presión	
  crítica	
  son	
  los	
  valores	
  de	
   la	
  temperatura	
  y	
   la	
  presión	
  en	
  el	
  punto	
  
crítico.	
  
	
  
	
  
Figura	
  	
  11.	
  Diagrama	
  de	
  fases	
  P-­‐T.	
  
	
  
En	
  la	
  Figura	
  12	
  se	
  puede	
  apreciar	
  un	
  diagrama	
  T	
  –	
  V	
  con	
  tres	
  regiones:	
  la	
  región	
  del	
  
lado	
   izquierdo	
   de	
   la	
   campana	
   que	
   corresponde	
   a	
   la	
   fase	
   líquida,	
   la	
   región	
   del	
   lado	
  
derecho	
  de	
  la	
  campana	
  que	
  corresponde	
  a	
  la	
  fase	
  vapor	
  y	
  la	
  región	
  que	
  se	
  halla	
  dentro	
  
de	
   la	
   campana	
   que	
   corresponde	
   a	
   una	
  mezcla	
   de	
   líquido	
   +	
   vapor.	
   	
   A	
   la	
   línea	
   de	
   la	
  
campana	
   que	
   baja	
   hacia	
   la	
   izquierda	
   del	
   punto	
   crítico	
   se	
   le	
   llama	
   línea	
   de	
   líquido	
  
saturado	
  y	
  a	
   la	
   línea	
  que	
  baja	
  hacia	
   la	
  derecha	
  del	
  punto	
  crítico	
  se	
   le	
   llama	
   línea	
  de	
  
vapor	
  saturado.	
  	
  Con	
  ayuda	
  de	
  este	
  diagrama	
  se	
  puede	
  saber	
  con	
  precisión	
  cuál	
  es	
  el	
  
	
  	
  	
  	
  	
  	
  31	
  
	
  
volumen	
  del	
   líquido	
  y	
  cuál	
  es	
  el	
  volumen	
  del	
  vapor	
  en	
  equilibrio	
  a	
  una	
   temperatura	
  
dada.	
  	
  
	
  
	
  
Figura	
  	
  12.	
  Diagrama	
  de	
  fases	
  T	
  -­‐	
  V.	
  
	
  
Otro	
   ejemplo	
   de	
   separación	
   de	
   fases	
   se	
   muestra	
   en	
   el	
   diagrama	
   presión	
   contra	
  
volumen	
  P-­‐V	
  (Figura	
  13).	
  	
  
	
  
Figura	
  	
  13.	
  Diagrama	
  de	
  fases	
  P-­‐V.	
  
	
  
Este	
  diagrama	
  puede	
  dividirse	
  trazando	
  una	
  curva	
  a	
  través	
  de	
  los	
  puntos	
  de	
  transición	
  
de	
   fase	
   de	
   cada	
   isoterma.	
   Según	
   esto,	
   un	
   sistema	
   con	
   presión	
   P	
   y	
   volumen	
   V	
  
correspondientes	
  a	
  la	
  porción	
  inferior	
  derecha	
  se	
  encuentra	
  en	
  la	
  fase	
  gaseosa,	
  en	
  la	
  
porción	
   superior	
   de	
   la	
   izquierda,	
   estará	
   en	
   la	
   fase	
   líquida	
   y	
   dentro	
   del	
   lugar	
  
32	
   	
  
	
  
geométrico	
  en	
  forma	
  de	
  parábola	
  invertida	
  y	
  abajo	
  del	
  punto	
  crítico	
  estará	
  constituido	
  
por	
   una	
   mezcla	
   de	
   fases	
   líquida	
   y	
   gaseosa	
   cuya	
   composición	
   sigue	
   la	
   regla	
   de	
   la	
  
palanca.	
  Una	
  isoterma	
  obtenida	
  de	
  este	
  diagrama	
  es	
  la	
  que	
  se	
  muestra	
  en	
  la	
  Figura	
  14.	
  	
  
	
  
Figura	
  	
  14.	
  Isoterma	
  de	
  un	
  diagrama	
  de	
  fases	
  	
  de	
  P-­‐V.	
  	
  
	
  
El	
  primer	
  aspecto	
  a	
  notar	
  de	
   la	
  Figura	
  14	
  es	
  que	
   la	
  Ecuación	
  14	
  no	
  se	
   cumple	
  en	
  el	
  
tramo	
  C-­‐D-­‐E	
  y	
  por	
  lo	
  tanto	
  no	
  representa	
  una	
  zona	
  estable,	
  lo	
  cual	
  quiere	
  decir	
  que	
  en	
  
esta	
   sustancia	
   debe	
   producirse	
   una	
   transición	
   de	
   fase.	
   Los	
   estados	
   previos	
   a	
   B	
   y	
  
posteriores	
  a	
  F	
  (por	
  ejemplo	
  A	
  y	
  G)	
  corresponden	
  a	
  situaciones	
  de	
  equilibrio	
  estable.	
  
En	
   los	
  estados	
  C	
  y	
  E	
  se	
  encuentra	
  un	
  mínimo	
  y	
  un	
  máximo	
  respectivamente.	
  Y	
  a	
   los	
  
puntos	
   D,	
   B	
   y	
   F	
   corresponde	
   un	
   mismo	
   valor	
   de	
   presión	
   Pt.	
   Esto	
   implica	
   que	
   el	
  
segmento	
  B-­‐C-­‐D	
  con	
  área	
  sombreada	
  es	
  igual	
  al	
  área	
  sombreada	
  del	
  segmento	
  D-­‐E-­‐F	
  
(regla	
  de	
  la	
  palanca)	
  y	
  justo	
  esta	
  presión	
  es	
  donde	
  ocurre	
  la	
  transición	
  de	
  fases.	
  	
  	
  	
  
Otra	
  representación	
  de	
  los	
  estados	
  correspondientes	
  a	
  cada	
  una	
  de	
  las	
  fases	
  posibles	
  
de	
  un	
  sistema	
  es	
  el	
  diagrama	
  de	
  potencial	
  químico	
  contra	
  presión	
  µ	
  –	
  P	
  (Figura	
  15).	
  En	
  
esta	
  gráfica	
   se	
  observan	
  varias	
   isotermas	
  donde	
  a	
  medida	
  que	
   crece	
   la	
   temperatura	
  
desaparecen	
  las	
  funciones	
  con	
  valor	
  triple	
  de	
  P.	
  	
  
	
  
	
  	
  	
  	
  	
  	
  33	
  
	
  
	
  
Figura	
  	
  15.	
  Diagrama	
  de	
  fases	
  	
  µ-­‐P.	
  	
  
	
  
En	
  la	
  Figura	
  16	
  se	
  muestra	
  la	
  primer	
  isoterma	
  del	
  diagrama	
  anterior.	
  	
  	
  
	
  	
  
Figura	
  	
  16.	
  Isoterma	
  del	
  diagrama	
  de	
  fases	
  	
  µ-­‐P.	
  	
  
	
  
A	
   medida	
   que	
   la	
   presión	
   aumenta	
   se	
   hacen	
   asequibles	
   al	
   sistema	
   tres	
   estados	
   con	
  
valores	
  iguales	
  de	
  P	
  y	
  T,	
  como	
  por	
  ejemplo,	
  para	
  el	
  punto	
  G	
  de	
  la	
  Figura	
  16.	
  De	
  estos	
  
tres	
  estados,	
  el	
  de	
  hasta	
  arriba	
  es	
  inestable,	
  recordando	
  que	
  la	
  zona	
  C-­‐D-­‐E	
  en	
  la	
  Figura	
  
14	
  no	
  cumple	
  el	
  criterio	
  de	
  estabilidad.	
  Los	
  dos	
  valores	
  de	
  potencial	
  químico	
  para	
  el	
  
punto	
  G,	
   el	
  de	
  hasta	
  abajo	
  y	
  el	
  de	
  en	
  medio	
   son	
  estados	
  estables,	
  pero	
  de	
  estos	
  dos	
  
valores	
   el	
   sistema	
   en	
   realidad	
   selecciona	
   el	
   estado	
   más	
   bajo	
   que	
   corresponde	
   a	
   la	
  
34	
   	
  
	
  
curva	
  del	
   gas.	
   Cuando	
   sigue	
   elevándose	
   la	
   presión	
  del	
   sistema,	
   llega	
   a	
   alcanzarse	
   el	
  
punto	
  único	
  D.	
  En	
  este	
  punto,	
  la	
  superficie	
  µ	
  se	
  corta	
  a	
  sí	
  misma	
  y	
  al	
  igual	
  que	
  la	
  Figura	
  
14	
  corresponde	
  al	
  punto	
  en	
  que	
  se	
  presenta	
  la	
  transición	
  de	
  fases.	
  	
  
2.2 Cálculo	
  de	
  equilibrio	
  de	
  los	
  sistemas	
  termodinámicos	
  
	
  
Para	
  muchas	
  sustancias	
  simples	
  que	
  sufren	
  transiciones	
   líquido	
  –	
  vapor,	
   la	
  ecuación	
  
de	
  estado	
  de	
  van	
  der	
  Waals	
  representa	
  bastante	
  bien	
  sus	
  propiedades	
  por	
  encima	
  y	
  
por	
  debajo	
  del	
  punto	
  deebullición.	
  La	
  ecuación	
  empírica	
  interpolada	
  es:	
  
P = NRT
V − bN
−
N 2a
V 2
	
   	
   	
   	
   	
   (16)	
  
	
  
donde	
  a	
  y	
  b	
  son	
  constantes	
  a	
  determinar	
  empíricamente	
  para	
  el	
  sistema	
  particular	
  de	
  
que	
  se	
  trate23.	
  
Otra	
   opción	
  de	
   calcular	
   estas	
   propiedades	
   es	
   la	
   ya	
  mencionada	
   técnica	
  de	
  dinámica	
  
molecular.	
   En	
   este	
   capítulo	
   se	
   mostrará	
   la	
   construcción	
   de	
   todos	
   los	
   diagramas	
  
estudiados	
   en	
   la	
   sección	
   2.1	
   para	
   átomos	
   de	
   argón	
   con	
   dinámica	
   molecular	
   y	
   en	
  
algunos	
   casos	
   se	
   realizará	
   una	
   comparación	
   con	
   datos	
   obtenidos	
   de	
   la	
   página	
   del	
  
National	
  Institute	
  of	
  Standards	
  and	
  Technology	
  (NIST)24.	
  	
  Los	
  mismos	
  diagramas	
  pero	
  
para	
  cadenas	
  lineales	
  de	
  átomos	
  de	
  argón	
  también	
  serán	
  mostrados.	
  
2.3 Equilibrio	
  de	
  fases	
  	
  
	
  
Las	
  simulaciones	
  de	
  dinámica	
  molecular	
  proveen	
  suficiente	
  información	
  para	
  calcular	
  
propiedades	
   en	
   equilibrio	
   de	
   fases7.	
  Más	
   aún,	
   de	
   calcular	
   propiedades	
   en	
   equilibrio	
  
para	
   cadenas	
   de	
   Lennard-­‐Jones8,25.	
   Esto	
   se	
   logra	
   creando	
   una	
   película	
   líquida	
   que	
  
coexiste	
   con	
   su	
   vapor	
   a	
   las	
   mismas	
   condiciones	
   de	
   temperatura	
   y	
   presión26.	
   Las	
  
densidades	
   de	
   bulto	
   para	
   ambas	
   fases	
   se	
   pueden	
   calcular	
   mediante	
   un	
   perfil	
   de	
  
densidades.	
   Con	
   la	
   finalidad	
   de	
   tener	
   una	
   mejor	
   visualización,	
   y	
   un	
   mayor	
   espacio	
  
	
  	
  	
  	
  	
  	
  35	
  
	
  
para	
  que	
  el	
  vapor	
  se	
  desarrolle,	
   se	
  utiliza	
  una	
  caja	
  con	
   longitud	
  en	
  dirección	
  x	
   cinco	
  
veces	
   mayor	
   que	
   las	
   demás.	
   Para	
   obtener	
   una	
   película	
   líquida	
   inicial	
   se	
   utiliza	
   un	
  
campo	
  gravitacional.	
  Para	
   lograr	
  esto,	
  a	
   la	
  aceleración	
  en	
  dirección	
  x	
  se	
   le	
  añade	
  un	
  
término	
  formado	
  por	
  una	
  constante	
  de	
  “gravedad”	
  multiplicada	
  por	
  el	
  signo	
  contrario	
  
de	
  la	
  posición	
  en	
  x	
  de	
  la	
  partícula	
  (Figura	
  17).	
  
	
  
Figura	
  	
  17.	
  Película	
  líquida	
  obtenida	
  con	
  un	
  campo	
  gravitacional.	
  
	
  
	
  
Figura	
  	
  18.	
  Película	
  líquida	
  al	
  centro	
  coexistiendo	
  con	
  su	
  vapor	
  en	
  equilibrio.	
  
Se	
  corren	
  varios	
  pasos	
  de	
  simulación	
  hasta	
  lograr	
  un	
  film	
  definido	
  y	
  una	
  vez	
  formada	
  
la	
  película	
  se	
  apaga	
  el	
  campo	
  gravitacional	
  y	
  se	
  corrige	
  la	
  velocidad	
  en	
  dirección	
  x,	
  con	
  
la	
  finalidad	
  de	
  que	
  el	
  sistema	
  por	
  sí	
  mismo	
  mantenga	
  la	
  película	
  líquida	
  coexistiendo	
  
con	
  su	
  vapor	
  hasta	
  llegar	
  al	
  equilibrio.	
  También	
  se	
  estabiliza	
  la	
  temperatura	
  y	
  se	
  hace	
  
!gx$!gx$
x$
y$
z$
-5 0 5-20 0 20
-5
0
5
x"
y"
z"
-505
-20 0 20
-5
0
5
36	
   	
  
	
  
una	
   corrección	
   debido	
   al	
   desplazamiento	
   del	
   centro	
   de	
   masa.	
   En	
   la	
   Figura	
   18	
   se	
  
muestra	
  la	
  película	
  líquida	
  al	
  centro	
  y	
  su	
  vapor	
  coexistiendo	
  en	
  equilibrio	
  en	
  la	
  caja	
  de	
  
simulación.	
  
Para	
   el	
   caso	
   de	
   cadenas	
   lineales,	
   se	
   sigue	
   el	
   mismo	
   procedimiento	
   y	
   se	
   obtienen	
  
películas	
  líquidas	
  en	
  equilibrio	
  con	
  su	
  vapor.	
  En	
  la	
  Figura	
  19	
  se	
  muestran	
  los	
  sistemas	
  
líquido	
  –	
  vapor	
  para	
  cadenas	
  de	
  2,	
  4,	
  8	
  y	
  16	
  átomos.	
  	
  
	
  
Figura	
  	
  19.	
  Equilibrio	
  líquido	
  –	
  vapor	
  para	
  2,	
  4,	
  8	
  y	
  16	
  átomos	
  en	
  la	
  cadena	
  lineal.	
  
	
  
	
  
	
  
	
  
	
  
	
  	
  	
  	
  	
  	
  37	
  
	
  
2.3.1	
  Perfil	
  de	
  densidades	
  
	
  
Después	
   de	
   obtener	
   la	
   película	
   líquida	
   y	
   el	
   vapor	
   en	
   equilibrio,	
   se	
   deja	
   correr	
   el	
  
programa	
   y	
   cada	
   mil	
   pasos	
   se	
   calculan	
   perfiles	
   de	
   densidades,	
   que	
   al	
   final	
   se	
  
promedian	
   para	
   poder	
   calcular	
   las	
   densidades	
   de	
   bulto26.	
   Para	
   realizar	
   el	
   perfil	
   de	
  
densidades	
  se	
  divide	
  la	
  caja	
  de	
  simulación	
  en	
  100	
  partes	
  iguales	
  y	
  se	
  cuenta	
  cuántas	
  
partículas	
   hay	
   en	
   cada	
   sección.	
   El	
   programa	
   se	
   termina	
   hasta	
   que	
   el	
   sistema	
   sea	
  
estable,	
   es	
   decir,	
   que	
   los	
   valores	
   de	
   densidad	
   del	
   líquido	
   y	
   densidad	
   del	
   vapor	
   no	
  
cambien.	
   En	
   la	
   Figura	
   20	
   se	
   muestra	
   el	
   perfil	
   de	
   densidades	
   adimensional	
   que	
  
corresponde	
  a	
  la	
  Fig.	
  19	
  (para	
  recordar	
  cómo	
  se	
  calcula	
  la	
  densidad	
  adimensional	
  ver	
  
el	
   Apéndice	
   1).	
   Se	
   puede	
   observar	
   que	
   la	
   parte	
   superior	
   corresponde	
   a	
   la	
   parte	
   de	
  
mayor	
  densidad,	
  que	
  en	
  este	
  caso	
  es	
  el	
  film	
  líquido	
  	
  y	
  las	
  dos	
  partes	
  de	
  abajo	
  que	
  son	
  
prácticamente	
  iguales,	
  corresponden	
  a	
  la	
  parte	
  del	
  vapor	
  a	
  ambos	
  lados	
  del	
  film.	
  	
  	
  
	
   	
  	
  
Figura	
  	
  20.	
  Perfil	
  de	
  densidades	
  correspondiente	
  a	
  la	
  Figura	
  18.	
  
	
  
Para	
  calcular	
  el	
  valor	
  de	
   la	
  densidad	
  del	
   líquido	
  y	
  el	
  valor	
  de	
   la	
  densidad	
  del	
  gas,	
  se	
  
realiza	
  un	
  ajuste	
  tangencial	
  hiperbólico27,28:	
  
( ) ( ) ⎥⎦
⎤
⎢⎣
⎡ −−−+
d
)(2
2
1
2
1
vapliqvapliq
lxtanhρρρρ 	
   	
   	
   (17)	
  
-15 -10 -5 0 5 10 15
x
0.1
0.2
0.3
0.4
0.5
0.6
0.7
r
38	
   	
  
	
  
Donde	
   ρliq	
   y	
   ρvap	
   son	
   la	
   densidad	
   del	
   líquido	
   y	
   del	
   vapor	
   respectivamente,	
   l	
   es	
   el	
  
espesor	
   de	
   la	
   película	
   y	
  d	
   el	
   espesor	
   de	
   la	
   interface.	
   En	
   la	
   Figura	
   21	
   se	
  muestra	
   el	
  
ajuste	
  sobre	
  el	
  perfil	
  de	
  densidades	
  de	
  la	
  Figura	
  20.	
  	
  
	
  
Figura	
  	
  21.	
  Ajuste	
  al	
  perfil	
  de	
  densidades	
  de	
  la	
  Figura	
  20.	
  
	
  
2.3.2	
  Diagrama	
  de	
  fases	
  T-­‐V	
  	
  
	
  
En	
  la	
  sección	
  anterior	
  se	
  explicó	
  cómo	
  se	
  construyen	
  los	
  perfiles	
  de	
  densidades	
  con	
  los	
  
datos	
   obtenidos	
   de	
   simulación	
   molecular.	
   Con	
   base	
   en	
   estos	
   perfiles	
   se	
   explicó	
  
también	
   cómo	
   calcular	
   las	
   densidades	
   de	
   bulto	
   de	
   las	
   fases	
   líquido	
   y	
   vapor	
   a	
   una	
  
determinada	
   temperatura.	
   Ahora	
   se	
   necesitan	
   hacer	
   las	
   mismas	
   simulaciones	
  
anteriores	
  pero	
  variando	
  únicamente	
  el	
  valor	
  de	
  la	
  temperatura.	
  Se	
  obtendrán	
  datos	
  
de	
  densidad	
  de	
  líquido	
  y	
  vapor	
  en	
  equilibrio	
  para	
  diferentes	
  temperaturas	
  y	
  con	
  estos	
  
datos	
   se	
  puede	
  graficar	
  un	
  diagrama	
  de	
   fases	
   como	
  el	
  mostrado	
  en	
   la	
  Figura	
  12.	
  Es	
  
importante	
  recordar	
  que	
  el	
  volumen,	
  para	
  fines	
  prácticos,	
  es	
  únicamente	
  el	
  inverso	
  de	
  
la	
  densidad	
  y	
  viceversa.	
  En	
  la	
  Figura	
  22	
  se	
  muestra	
  el	
  diagrama	
  T-­‐	
  ρ	
  y	
  en	
  laFigura	
  23	
  
se	
  muestra	
  el	
  diagrama	
  T-­‐	
  V	
  para	
  un	
  solo	
  átomo	
  de	
  argón	
  en	
  la	
  cadena.	
  Asimismo	
  se	
  
muestra	
   la	
   comparación	
   a	
   las	
   mismas	
   condiciones	
   con	
   resultados	
   reportados	
   en	
   el	
  
NIST24.	
  	
  
	
  
-20 -10 0 10 20
x
0.1
0.2
0.3
0.4
0.5
0.6
0.7
r
	
  	
  	
  	
  	
  	
  39	
  
	
  
	
  
Figura	
  	
  22.	
  Diagrama	
  de	
  fases	
  T	
  –	
  ρ	
  para	
  un	
  átomo	
  de	
  argón	
  (m=1).	
  Los	
  símbolos	
  son	
  los	
  
resultados	
  calculados	
  con	
  DM	
  y	
  las	
  líneas	
  son	
  los	
  valores	
  reportados	
  por	
  el	
  NIST.	
  
	
  
Figura	
  	
  23.	
  Diagrama	
  de	
  fases	
  T	
  –V	
  para	
  un	
  átomo	
  de	
  argón	
  (m=1).	
  Los	
  símbolos	
  son	
  los	
  
resultados	
  calculados	
  con	
  DM	
  y	
  las	
  líneas	
  son	
  los	
  valores	
  reportados	
  por	
  el	
  NIST.	
  
0.60$
0.70$
0.80$
0.90$
1.00$
1.10$
1.20$
1.30$
1.40$
1.50$
0.0$ 0.1$ 0.2$ 0.3$ 0.4$ 0.5$ 0.6$ 0.7$ 0.8$
T*#
rho*#
m=1#
Curva$del$liquido$
Curva$del$vapor$
Curva$liquido$NIST$
Curva$vapor$NIST$
0.80$
0.90$
1.00$
1.10$
1.20$
1.30$
0.0$ 20.0$ 40.0$ 60.0$ 80.0$ 100.0$ 120.0$ 140.0$ 160.0$ 180.0$ 200.0$
T*#
V*#
m=1#
Curva$del$liquido$
Curva$del$vapor$
Curva$liquido$NIST$
Curva$vapor$NIST$
40	
   	
  
	
  
En	
  la	
  Figura	
  24	
  se	
  muestra	
  el	
  diagrama	
  de	
  fases	
  T	
  –	
  ρ	
  para	
  todas	
  las	
  cadenas	
  lineales	
  y	
  
en	
  la	
  Figura	
  25	
  se	
  muestra	
  el	
  diagrama	
  de	
  fases	
  T	
  –	
  V	
  también	
  para	
  todas	
  las	
  cadenas.	
  
	
  
Figura	
  	
  24.	
  Diagrama	
  de	
  fases	
  T	
  –	
  ρ	
  para	
  las	
  diferentes	
  cadenas	
  lineales.	
  
	
  
Figura	
  	
  25.	
  Diagrama	
  de	
  fases	
  T	
  –	
  V	
  para	
  las	
  diferentes	
  cadenas	
  lineales.	
  	
  
0.0#
0.5#
1.0#
1.5#
2.0#
2.5#
3.0#
3.5#
0.0# 0.1# 0.2# 0.3# 0.4# 0.5# 0.6# 0.7# 0.8# 0.9#
T*#
rho*#
m=1#
m=2#
m=4#
m=8#
m=16#
0.0#
0.5#
1.0#
1.5#
2.0#
2.5#
3.0#
3.5#
0.0# 50.0# 100.0# 150.0# 200.0# 250.0# 300.0#
T*#
V*#
m=1#
m=2#
m=4#
m=8#
m=16#
	
  	
  	
  	
  	
  	
  41	
  
	
  
De	
  las	
  Figuras	
  23	
  y	
  24	
  se	
  puede	
  afirmar	
  que	
  el	
  programa	
  que	
  se	
  diseñó	
  para	
  el	
  cálculo	
  
de	
   equilibrio	
   de	
   fases	
   es	
   correcto	
   porque	
   tiene	
   buena	
   concordancia	
   con	
   los	
   datos	
  
publicados	
  por	
  una	
  fuente	
  confiable	
  como	
  es	
  el	
  NIST24.	
  	
  
De	
  las	
  Figuras25	
  y	
  26	
  se	
  puede	
  concluir	
  que	
  a	
  mayor	
  número	
  de	
  átomos	
  en	
  la	
  cadena	
  
lineal,	
  mayor	
  es	
  la	
  temperatura	
  de	
  la	
  curva	
  de	
  saturación	
  y	
  por	
  ende	
  de	
  la	
  temperatura	
  
crítica.	
  	
  
En	
  la	
  Tabla	
  1	
  se	
  muestran	
  los	
  valores	
  de	
  las	
  gráficas	
  anteriores	
  obtenidos	
  para	
  todas	
  
las	
  cadenas	
  y	
  para	
  varias	
  temperaturas.	
  La	
  temperatura	
  crítica	
  reportada	
  es	
  un	
  simple	
  
estimado	
   y	
   se	
   observa	
   experimentalmente	
   en	
   la	
   práctica,	
   porque	
   al	
   correr	
   la	
  
simulación	
   con	
   esta	
   temperatura,	
   se	
   obtienen	
   valores	
   absurdos	
   que	
   no	
   concuerdan	
  
con	
   la	
   tendencia	
   de	
   datos	
   reportados.	
   En	
   este	
   momento	
   se	
   cree	
   que	
   se	
   alcanzó	
   la	
  
temperatura	
  crítica.	
  	
  	
  
	
  
Tabla	
  1.	
  Valores	
  de	
  densidad	
  de	
  líquido	
  y	
  densidad	
  del	
  vapor	
  en	
  equilibrio	
  para	
  diferentes	
  
cadenas	
  y	
  diferentes	
  temperaturas.	
  	
  
T* T* T*
0.818 ±&0.008 0.0037 ±&0.0002 0.75 0.740 ±&0.002 0.0072 ±&0.0003 1.2 0.682 ±&0.003 0.0042 ±&0.0001 1.6
0.795 ±&0.002 0.0062 ±&0.0002 0.8 0.721 ±&0.002 0.0102 ±&0.0005 1.25 0.666 ±&0.005 0.0058 ±&0.0003 1.65
0.771 ±&0.005 0.0098 ±&0.0005 0.85 0.702 ±&0.006 0.0141 ±&0.0002 1.3 0.650 ±&0.006 0.0078 ±&0.0002 1.7
0.748 ±&0.001 0.0148 ±&0.0008 0.9 0.682 ±&0.008 0.0190 ±&0.0001 1.35 0.633 ±&0.007 0.0103 ±&0.0006 1.75
0.723 ±&0.002 0.0214 ±&0.0001 0.95 0.661 ±&0.003 0.0251 ±&0.0008 1.4 0.616 ±&0.003 0.0135 ±&0.0009 1.8
0.697 ±&0.008 0.0300 ±&0.0003 1 0.639 ±&0.004 0.0328 ±&0.0006 1.45 0.598 ±&0.003 0.0173 ±&0.0008 1.85
0.669 ±&0.004 0.0412 ±&0.0002 1.05 0.616 ±&0.001 0.0422 ±&0.0002 1.5 0.579 ±&0.004 0.0219 ±&0.0001 1.9
0.639 ±&0.006 0.0558 ±&0.0006 1.1 0.590 ±&0.007 0.0538 ±&0.0004 1.55 0.559 ±&0.006 0.0274 ±&0.0002 1.95
0.604 ±&0.006 0.0750 ±&0.0008 1.15 0.561 ±&0.009 0.0683 ±&0.0001 1.6 0.538 ±&0.007 0.0340 ±&0.0003 2
0.563 ±&0.005 0.1013 ±&0.0003 1.2 0.528 ±&0.003 0.0867 ±&0.0002 1.65 0.516 ±&0.009 0.0419 ±&0.0006 2.05
0.507 ±&0.008 0.1409 ±&0.0004 1.25 0.489 ±&0.002 0.1108 ±&0.0009 1.7 0.492 ±&0.002 0.0514 ±&0.0004 2.1
1.3 1.8 2.4
T* T*
0.559 ±&0.002 0.0061 ±&0.0001 2.2 0.520 ±&0.005 0.0014 ±&0.0002 2.5
0.543 ±&0.001 0.0080 ±&0.0005 2.25 0.506 ±&0.003 0.0019 ±&0.0005 2.55
0.527 ±&0.004 0.0103 ±&0.0005 2.3 0.492 ±&0.007 0.0027 ±&0.0007 2.6
0.510 ±&0.008 0.0131 ±&0.0005 2.35 0.478 ±&0.001 0.0036 ±&0.0002 2.65
0.493 ±&0.005 0.0165 ±&0.0002 2.4 0.464 ±&0.003 0.0049 ±&0.0003 2.7
0.475 ±&0.003 0.0205 ±&0.0008 2.45 0.449 ±&0.004 0.0064 ±&0.0001 2.75
0.456 ±&0.002 0.0253 ±&0.0001 2.5 0.433 ±&0.005 0.0084 ±&0.0008 2.8
0.437 ±&0.008 0.0311 ±&0.0002 2.55 0.418 ±&0.007 0.0108 ±&0.0009 2.85
0.416 ±&0.009 0.0380 ±&0.0003 2.6 0.401 ±&0.008 0.0137 ±&0.0002 2.9
0.393 ±&0.004 0.0462 ±&0.0004 2.65 0.384 ±&0.003 0.0173 ±&0.0003 2.95
0.369 ±&0.002 0.0563 ±&0.0005 2.7 0.368 ±&0.002 0.0216 ±&0.0005 3
2.9 3.3Temperatura&crítica&aproximada: Temperatura&crítica&aproximada:
m=1 m=2 m=4
ρliq ρvap
Temperatura&crítica&aproximada:
ρliq ρvap ρliq ρvap
Temperatura&crítica&aproximada: Temperatura&crítica&aproximada:
m=8 m=16
ρliq ρvap ρliq ρvap
42	
   	
  
	
  
2.4 Presión	
  
	
  
Como	
  se	
  explicó	
  en	
  la	
  sección	
  2.1,	
  la	
  construcción	
  de	
  diagramas	
  de	
  fases	
  nos	
  facilita	
  la	
  
visualización	
  de	
  los	
  estados	
  de	
  equilibrio	
  estables	
  y	
  no	
  estables;	
  y	
  la	
  mayoría	
  de	
  estos	
  
diagramas	
  de	
  fases	
  son	
  funciones	
  de	
  la	
  presión.	
  Es	
  por	
  eso	
  que	
  es	
  necesario	
  explicar	
  
cómo	
  se	
  calcula	
  la	
  presión	
  con	
  la	
  técnica	
  de	
  dinámica	
  molecular.	
  
El	
  tensor	
  de	
  presión	
  de	
  un	
  fluido	
  en	
  reposo	
  es:	
  
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
−
−
−
=
zz
yy
xx
00
00
00
P
p
p
p
	
  	
   	
   	
   	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
   	
  	
  	
  	
  	
  (18)	
  
donde	
  cada	
  componente	
  del	
  tensor	
  de	
  presión	
  se	
  puede	
  calcular	
  en	
  DM	
  como	
  la	
  suma	
  
de	
   dos	
   contribuciones,	
   una	
   parte	
   cinética	
   y	
   la	
   otra	
   virial	
   (o	
   de	
   interacciones	
   entre	
  
partículas):	
  	
  
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+= ∑ ∑∑
>i i ij j
2
ij
ij
2
αi,αα )(
1
ir
rFmv
V
p
α
	
  	
  	
  	
   	
   	
  	
  	
  	
  	
  	
  	
  	
  
(19)
 
donde	
  α puede	
  ser	
   x,	
  y	
  o	
  z,	
  V	
  es	
  el	
  volumen	
  de	
   la	
  caja	
  y	
  vi,α es	
  el	
   componente	
  α de	
   la	
  
velocidad	
   de	
   la	
   partícula	
   i.	
   Cuando	
   el	
   sistema	
   está	
   en	
   equilibrio	
   y	
   no	
   existen	
   otros	
  
esfuerzos,	
  los	
  elementos	
  de	
  la	
  diagonal	
  del	
  tensor	
  son	
  iguales29.	
  
La	
  llamada	
  presión	
  hidrostática	
  instantánea30	
  se	
  puede	
  calcular	
  mediante	
  la	
  Ecuación:	
  
Peq = 1
3
Tr(P) =
Pxx
eq +Pyy
eq +Pzz
eq
3
=
1
V
NKT + 1
3
F(rij ) ⋅ rij
j>i
∑
i
∑
#
$
%%
&
'
((
	
  	
  	
  	
  	
  	
  	
  	
  	
  (20)	
  
La	
  presión	
  hidrostática	
  se	
  debe	
  corregir	
  debido	
  al	
  potencial	
  truncado	
  por	
  el	
  radio	
  de	
  
corte31.	
  	
  La	
  corrección	
  para	
  la	
  presión	
  es: