Introducci on a la Din amica Molecular Curso corto interactivoiqcelaya.itc.mx/~omar/Docencia_files/EjerciciosEOCG60Aniv.pdf · como archivos de estructura o de interacci on par. En

Introduccion a la Dinamica MolecularCurso corto interactivo

Dr. Edgar Omar Castrejon Gonzalez

IT-Celaya

19 de abril de 2018

EO Castrejon-Gonzalez (IT-Celaya) LAMMPS 19 abril 2018 1 / 54

1 IntroduccionDinamica Molecular

2 Ensayo de tension

3 Flujos Couette y Poiseuille 2D

4 NaCl

5 RDF en fluido LJ

6 Micelas en 2D

7 Micelas 3D

8 Equilibrio Vapor–Lıquido

9 Viscosidad de corte en Fluidos Complejos

10 Diferentes modelos de agua


Simulacion Molecular

Simulacion Molecular es una tecnica ampliamente empleada endiversas areas, como: Fısica, Quımica, Biologıa, Alimentos y Cienciade Materiales.

El principal objetivo de la Simulacion Molecular es resolvernumericamente los modelos matematicos que plantea la MecanicaEstadıstica sobre un problema real.

Con Simulacion Molecular se pueden obtener condiciones que sondifıciles de reproducir en un experimento, o entender lo que sucede anivel atomico.


Introduccion a LAMMPS

LAMMPS es un codigo paralelizado de simulacion por DinamicaMolecular, recientemente han incluido calculos por Montecarlo.

Esta codificado en C++ y puede ser compilado manualmente comosoftware independiente o como librerıa dentro de otro programa(Ver Manual).

Tambien se puede compilar en GPUs.

Se requiere software adicional para generar las estructuras iniciales ypara procesar los datos generados.


http://lammps.sandia.gov/doc/Manual.html

Introduccion

Para ejecutar una simulacion en LAMMPS se requiere un archivo deentrada (input) , bastara con ejecutar la siguiente instruccion paracorrer el programa: lmp serial.exe < input , dondelmp serial.exe es el nombre del ejecutable de Lammps.

Dentro del archivo input se pueden invocar otros archivos externos,como archivos de estructura o de interaccion par.

En el archivo input se escriben las instrucciones para calcularpropiedades e imprimir archivos con resultados especıficos.

Una vez terminada la simulacion, se genera por default un archivo desalida llamado log.lammps con informacion de la corrida.


Ensayo de tension

Se trata de un sistema constituido por esferas Lennard-Jones, cuyainteraccion esta dada por:

U(r) = 4ε

[(σijrij

)12

−(σijrij

)6]rc ≤ 2.5σ (1)

donde ε = 1.0 y σij.

Nombre de la carpeta: crack

Archivos requeridos: a) in.crack

Las interacciones estan definidas por los comandos: pair style

lj/cut 2.5 y pair coeff 1 1 1.0 1.0 2.5


¿Como funcionan algunos comandos de LAMMPS?

1 fix avetime Nevery Nrepeat Nfreq . For example, if Nevery=2,Nrepeat=6, and Nfreq=100, then values (6) on timesteps90,92,94,96,98,100 will be used to compute the final average ontimestep 100.


Instrucciones a seguir(a) Ejecutar Lammps usando:

lmp serial.exe < in.crack

(b) Visualizar el archivodump.lammps usando OVITO.

Se puede observar la propa-gacion de la fractura.


Flujos Couette y Poiseuille en un fluido LJ, en 2D

Nuevamente el sistema es constituido por esferas Lennard-Jones, cuyainteraccion esta dada por WCA:

U(r) = 4ε

[(σijrij

)12

−(σijrij

)6]rc ≤ 1.12246σ (2)


Nombre de la carpeta: 2flow

Archivos requeridos: a) in.flow.couette y b) in.flow.pois

El flujo se genera con los comandos: velocity upper set 3.0 0.0

0.0 para flujo Couette y fix 6 flow addforce 0.5 0.0 0.0

para flujo Poiseuille.


Instrucciones a seguirFlujo Couette

(a) Ejecutar Lammps usando:lmp serial.exe <in.flow.couette

(b) Visualizar el archivo dump.flow

usando OVITO.

Se puede observar, en la ani-macion de Ovito, el perfil lin-eal de velocidad.


Flujo Couette: Perfil de velocidadElabore la grafica usando xmgrace

0 5 10 15 20y Axis

0

0.5

1

1.5

2

2.5

3ve

loci

ty

Velocity profile in Couette flow


Instrucciones a seguirFlujo Poiseuille

(a) Ejecutar Lammps usando:lmp serial.exe <in.flow.pois

(b) Visualizar el archivodump.flowP usando OVITO.

Se puede observar, en la an-imacion de Ovito, el perfilparabolico de velocidad.


Flujo Poiseuille: Perfil de velocidadElabore la grafica usando xmgrace

00 0

1.1

1.34 1.32 1.34 1.37 1.35

0.709

00 5 10 15

y Axis

0

0.5

1

1.5

2ve

loci

ty

Velocity profile in Poiseuille Flow


Simulacion de NaCl en 3D

El sistema esta constituido por dos tipos de atomos Na y Cl, elpotencial de interaccion es EIM (embedded-ion method) propuestopor Zhou.

Para mayor informacion referente a la forma matematica delpotencial, visitar el siguiente sitio web.

Nombre de la carpeta: 3eim

Archivos requeridos: a) data.eim , b) in.eim y c) ffield.eim

El acomodo inicial aparece en el archivo de estructura data.eim .


http://lammps.sandia.gov/doc/pair_eim.html

Instrucciones a seguirNaCl

(a) Ejecutar Lammps usando:lmp serial.exe < in.eim

(b) Visualizar el archivo dump.eim

usando OVITO.

Se puede observar, el cambioen la estructura cristalina.


Calculo de la Funcion de Distribucion Radial (FDR) en unsistema LJ a diferentes Temperaturas

Se trata de un sistema en 3D constituido por esferas Lennard-Jones.

Las interacciones estan dadas por:

U(r) = 4ε

[(σijrij

)12

−(σijrij

)6]rc ≤ 2.5σ (3)


Nombre de la carpeta: 4rdf

Archivos requeridos: a) in.lj

Las interacciones estan definidas por los comandos: pair style

lj/cut 2.5 y pair coeff 1 1 1.0 1.0 2.5


Funcion de Distribucion Radial

Es una medida de laprobabilidad deencontrar una partıcula auna distancia r medidadesde la partıcula centralde referencia.

Para diferentes estados de agregacion:


Instrucciones a seguirRDF-LJ

1 Verificar que la temperatura adimensional sea 0.5: vim in.lj

2 Ejecutar Lammps usando: lmp serial.exe < in.lj

3 Depurar la salida y guardarla en un archivo denominado RDF1.rdfusando: tail -100 LJ.rdf > RDF1.rdf.

4 Seleccionar las columnas 2 y 3 y guardar el archivo en RDF2.rdfusando: cut -f 2-3 -d‘‘ ’’ RDF1.rdf > RDF2.rdf.

5 Elaborar la grafica en xmgrace usando: xmgrace RDF2.rdf.

6 Modificar el archivo de entrada in.LJ, cambiando la Temperatura a4.0; el archivo donde se guardaran los datos de RDF seguira siendoLJ.rdf.


Instrucciones a seguir

7 Repetir los pasos del 2 al 5 para generar la grafica de RDF enxmgrace:

1 lmp serial.exe < in.lj.2 tail -100 LJ.rdf > RDF3.rdf.3 cut -f 2-3 -d‘‘ ’’ RDF3.rdf > RDF4.rdf.4 xmgrace RDF4.rdf.

8 Ejecutar Lammps usando: lmp serial.exe < in.lj

9 Modificar el archivo in.LJ , cambiando T = 7.0 y ρ = 0.1. En elarchivo LJ.rdf se guardaran los datos de RDF.

10 Repetir los pasos del 2 al 5 para generar la grafica de RDF enxmgrace:

1 lmp serial.exe < in.lj.2 tail -100 LJ.rdf > RDF5.rdf.3 cut -f 2-3 -d‘‘ ’’ RDF5.rdf > RDF6.rdf.4 xmgrace RDF6.rdf.


Resultados RDF: Solido

0 0.5 1 1.5 2 2.5 3distance

0

1

2

3

4

5

g(r

)


Resultados RDF: Lıquido

0 0.5 1 1.5 2 2.5 3distance

0

0.5

1

1.5

2

g(r

)

Liquid, T = 4.0


Resultados RDF: Gas

0 0.5 1 1.5 2 2.5 3distance

0

0.5

1

1.5

2

g(r)

Gas, T = 7.0, density = 0.1


Micelas en 2D

Se trata de un sistema en 2D constituido por moleculas anfifılicas ensolucion.

Las interacciones entre partıculas no enlazadas afines esta dada por elpotencial L-J 12/6 con rc = 2.5σ.

Las interacciones entre partıculas no enlazadas que no son afines, seusa el potencial 12/6 con rc = 1.1228.

El enlace esta dado por un potencial armonico de la forma:E = k (r− r0)

2

Nombre de la carpeta: 5micelle2D

Archivos requeridos: a) in.micelle , b) def.micelle , c)micelle2d.f



1 Ver el archivo in.micelle usando: less in.micelle.

2 El archivo micelle2d.f es un codigo en Fortran para

generar el archivo de estructura. Se debe compilar

usando: gfortran micelle2d.f -o datagen. con esto se generaun ejecutable llamado datagen.

3 Correr el ejecutable usando como archivo de entrada def.micelle ycomo salida el archivo data.micelle usando: ./datagen <def.micelle > data.micelle

4 Ejecutar Lammps usando: lmp serial.exe < in.micelle

5 Visualizar en Ovito y generar una pelıcula.

6 Generar un archivo de estructura con 200 surfactantes de longitud 6(5 # of tails) en una placa de 40 40 (Modificar el def.micelley guardarlo como def.micelle2) .


Resultados: Parte 1

Inicio Fin


Resultados: Parte 2

Inicio Fin


Micelas 3D

Se trata de un sistema en 3D constituido por moleculas anfifılicas ensolucion.

Las interacciones entre partıculas no enlazadas afines esta dada por elpotencial L-J 12/6 con rc = 2.5σ.

Las interacciones entre partıculas no enlazadas que no son afines, seusa el potencial 12/6 con rc = 1.1228.

El enlace esta dado por un potencial armonico de la forma:E = k (r− r0)

2

Nombre de la carpeta: 6micelle3D

Archivos requeridos: a) in.micelle3D , b) def.micelle3D , c)Micelas3DCurso.f



1 El archivo Micelas3DCurso.f es un codigo en Fortran para

generar el archivo de estructura. Se debe compilar

usando: gfortran Micelas3DCurso.f -o datagen. con esto segenera un ejecutable llamado datagen3d.

2 Generar un archivo con 10 % aceite, para lo cual se debe modificar elarchivo def.micelle3d usando: vim def.micelle3d.

3 Ejecutar el programa: ./datagen3d < def.micelle3d >data.micelle3d

4 Ejecutar Lammps usando: lmp serial.exe < in.micelle3d

5 Visualizar en Ovito y generar una pelıcula.


Resultados

Inicio Fin


Equilibrio Vapor–Lıquido en un fluido LJ

Se trata de un sistema en 3D constituido por esferas LJ.

Las interacciones estan dada por el potencial L-J 12/6 con rc = 2.5σ,σ = ε = 1.0.

Nombre de la carpeta: 7evl

Archivos requeridos: a) inLV.LJ

Abrir el archivo inLV.LJ para analizarlo, la temperatura debe ser 0.7.

Ejecutar el programa: lmp serial.exe < inLV.LJ

Visualizar el archivo de salida dump.nvt con Ovito.

Efectuar la grafica del perfil de densidad usando: tail -504

densityprofx.out > denx.out para depurar el archivo y:

xmgrace -block denx.out -bxy 2:4


Resultados: EVL, T = 0.5

Configuracion final

Perfil de densidad



Configuracion finalPerfil de densidad



Configuracion final

Perfil de densidad


IntroduccionLos polımeros estan categorizados dentro de los fluidos noNewtonianos.

Los homopolımeros presentan un adelgazamiento de corte(shear-thinning) al someterse a flujo.

La curva de flujo tıpica se muestra a continuacion:


Desarrollo

El calculo de la viscosidad de corte esta dada por:

η = −〈Pxy〉γ

(4)

Donde η es la viscosidad de corte, 〈Pxy〉 es un promedio de ensamblede la componente xy del tensor de presion, calculado por laEcuacion 5, y γ es la tasa de corte.

Pkl = −1

V

∑i

pkiplimi

+1

2

∑i

∑j>i

rkijFlij

(5)


Ejercicio

Calcular la viscosidad de corte en LAMMPS [1] para un sistema constituidopor 30 cadenas de 100 sitios de acuerdo al siguiente modelo [2]:



Archivos requeridos: a) def.chain b) chain.f c) in.fene

Generar un archivo de estructura con 100 cadenas, para lo cual esnecesario modificar el def.chain

(a) Compilar el programa chain.f usando: gfortran chain.f

-o chain

(b) Ejecutar el programa: ./chain < def.chain >data.chain3

(c) Ejecutar Lammps usando: lmp serial.exe < in.fene

(d) Calcular la viscosidad de corte con η = − 〈Pxy〉γ , para

este caso = -(-0.7152/0.03) = 23.84


Continuacion

(e) Construir una curva de viscosidad de corte vs tasa de corte (γ)(ejercicio grupal).

(f) Construir una curva de Rg vs γ (ejercicio grupal). Valores de tasa decorte: 0.0001, 0.005, 0.001, 0.01, 0.05, 0.07, 0.1, 0.15, 0.27.

(g) Visualizar las funciones de distribucion radial usando:I tail -100 rdfEq.out > rdf.outI xmgrace -block rdf.out -bxy 2:3


Resultados:

Configuracion final EMDConfiguracion final NEMD


Resultados:


Diferentes modelos de aguaSPCE, TIP3P, TIP4P, TIP5P


Ejercicio: Comparar RDF, MSD y densidad de tres modelosde agua

En la carpeta 9water encontrara 3 carpetas: a) SPCE, b) TIP3P y c)TIP4P

Instrucciones:

1 Ingresar a la carpeta SPCE

2 Ejecutar LAMMPS usando: lmp serial.exe < input.dat .3 Elaborar grafica de RDF:

a) tail -1000 wat.rdf > rdf.out

b) xmgrace -block rdf.out -bxy 2:3

4 Elaborar grafica de MSD: xmgrace wat.msd y calcular el coeficientede difusion. Ans. 0.000782A/fs = 7.82× 10−9cm2/s

5 Calcular la densidad promediando valores del archivo wat.dens Ans.1.03g/cm3


Resultados SPCE: Configuraciones

Configuracion Inicial Configuracion Final


Resultados SPCE: RDF


Resultados SPCE: MSD


Resultados modelo TIP3P

Instrucciones:

1 Ingresar a la carpeta TIP3P







Resultados TIP3P: Configuraciones



Resultados TIP3P: RDF


Resultados TIP3P: MSD


Resultados modelo TIP4P

Instrucciones:

1 Ingresar a la carpeta TIP4P







Resultados TIP4P: Configuraciones



Resultados TIP3P: RDF


Resultados TIP4P: MSD


Referencias

S. J. Plimpton, (1995). J. Comp. Phys. 117, 1.

Castrejon-Gonzalez O., J. Castillo - Tejas, O. Manero & J. F. J. Alvarado, (2013).Structure factor and rheology of chain molecules from molecular dynamics, J.Chem. Phys., 138, 184901-1-11.


Documents

Introducci on a la Din amica Molecular Curso corto interactivoiqcelaya.itc.mx/~omar/Docencia_files/EjerciciosEOCG60Aniv.pdf · como archivos de estructura o de interacci on par. En