|
|||||||||||
COMPUTER SIMULATION OF LIQUIDS M.P. ALLEN D.J. TILDESLEY OXFORD SCIENCE PUBLICATIONS Este es un libro para enseñar a las personas que esten interesadas en la creación y diseño de simulaciones del comportamiento de liquidos atómicos y moleculares. Este libro será útil para lños estudiantes de primer año de graduación. Este tipo de programas de simulaciones se han llevao a cabo alo largo de generaciones de estudiantes y maestros, los cuales se han dado a la tarea de completar y recompilar estos programas. En los primeros seis capítulos del libro se dan a conocer los tacnicas de simulación, llamdas MONTE CARLO Y DE DINÁMICA MOLECULAR. Introducción Una breve historia de la simulacion computacional En los modelos antiguos de líquidos, se sontenía que los líquidos eran una aglomeración de un gran número de moléculas, y hasta la fecha hay algo interesante en el estudio de ensambles de bolas mecánicas, puestan en movimiento por vibración mecánica. De cualquier manera el uso de una gran cantidad de objetos físicos para representar molèculas puede consumir mucho tiempo, hay limitaciones obvias en el tipo de interacción entre ellas, y os efectos de gravedad no pueden ser despreciados. La extensión natural de esta aproximación es el uso de un modelo matemático, en lugar de un modelo físico, y ejecutar el análisis por computadora. Desde hace 50 años se utilizaron computadoras muy poderosas, como la computadora ENIAC de Los Alamos National Laboratories. Loas trabajos más antiguos fundaron lo que ahora se conoce como Método de MonteCarlo , que abreviaremos como MC. Las primeras simulkaciones fueropn bastante idealistas, utilizando para ellos modelos como moléculas esféricas, se podría decir que el primer método llevado a cabo fue el del gas ideal como primera aproximación. Posteriormente se llevaron a cabo simulaciones con una interacción del potencial de Lennard-Jones. Esto hizo posible ue se compararan resultados de simulación por computadora y el gas argon. Una técnica difreente se ha utilizado para obtener las propiedades dinámicas de una gran cantidad de partículas. Dinámica Molécular (MD), es el término usado para describir la solución de la ecuación clásica de movimiento de Newton (newton Equation) de un conjunto de moléculas. Fué Rahman 1964, quien se dedicó a resolver la ecuación de movimiento de Newton bajo la interacción de un potencial de Lennard- Jones, o con artículas de Lennard-Jones. Posteriormente de este desarrollo de la dinámica molecular, se desarrollo ampliamente la simulación computacional. SIMULACION COMPUTACIONAL: MOTIVACIONES Y APLICACIONES Algunos problemas en mecánica estadística son exactamente solubles, uno de los más famosos es el modelo del gas ideal PV = NKT. Algunos problemas en mecánica estadística, que no tienen una solución exacta, pueden resolverse si se utiliz un análisis de aproximaciones. Es importante hacer notar que es de gran ayuda poder encontrar un modelo que se pueda comportar como un líquido (el cual debe ser muy idelizado). Las simulaciones computacionales juegan un rol muy importante para obtener esencialmente, resultados exactos para problemas en mecánica estadística, que de otra manera solo podían resolverse por métodos aproximados, o simplemente no se podían resolver. En el siguiente diagrama se muestra, el desarrollo de la relación entre un modelo computacional y un modelo teórico, así como un diagrama de flujo para comparar los modelos teoricos y experimentales. La simualción computacional provee una ruta directa de los detalles microscópicos de un sistema (las masas de los átomos, la interacción entre ellos, la geometría molecular, etc.) a las propiedades macroscópicas de interés experimental (la ecuación de estado, los coeficientes de transporte, parametros de orden estructural, y algunos otros). Así como son de interés académico, este tipo de información es útil tecnolóogicamente hablando. Sería practicamente imposible llevar a cabo experimentos bajo condiciones extremas de teperatura y presión, cuando una simulación computacional de un plasma a temperatura extrema, de un reactor nuclear, o de un sistema planetario es perfectamente posible. Finalmente, cuando la velocidad de los eventos moleculares es en si mismo una dificultad experimental, este no es problema para el simulador. Un amplio rango de fenómenos físicos desde la escala molecular a la escala galáctica, puede ser estudiada usando la técnica de simulación computacional. En la mayoría del libro, estaremos al tanto de los detalles de llevar a cabo una simulación computacional. En el resto del capítulo, mostraremos como ingresar información a la computadora, yy en el capítulo 2 mostraremos como obtener información de la computadora. MODELACIÓN DE SISTEMAS Y POTENCIALES DE INTERACCIÓN. INTRODUCCIÓN. En la mayoría del libro, el estado microscópico de un sistema puede ser especificado en términos de las posiciones y momentos de un conjunto constitutivo de partículas: los átomos y moléculas. Con la aproximación de Bohrn-Oppenheimer, es posible expresar el hamiltoniano de un sistema en función de las variables nucleares, el movimiento de los electrones han sido promediados. Haciendo la aproximación adicional donde una descripción clásica es adecuada, podemos escribir el hamiltoniano H de un sistema de N moléculas como la suma de ls funciones de energia cinética y potencial de un conjunto de coordenadas generalizadas qi y de momentos pi para cada molécula i. Adoptando una notación condensada: tenemos: Las coordenadas generalizadas q pueden ser simplemente el conjunto de coordenadas cartesianas de cada átomo ( el núcleo) en el sistema, pero, como veremos, es muchas veces más útil tratar moléculas como cuerpos rígidos, en tal caso q consistirá de coordenadas cartesianas de cada centro de masa de cada molécula, en conjunto con variables , que especifican la orientación molecular. En cualquier caso, p se mantiene para un conjunto apropiado de momentos conjugados. Usualmente, la energía cinética H toma la forma: donde mi es la masa molecular, y el índice a varía sobre diferentes (x,y,z) componentes del momento de la molécula i. La energía potencial V contiene información interesante en relación con la interacción intermolecular: asumiendo que V (is fairly sensibly behaved), será posible construir, de H, una ecuación de movimiento (en la forma de Hamiltoniano, lLagrangiano o Newton) que gobernará totalmente la evolución temporal de un sistema y de todas sus propiedades mecánicas. La solución a esta ecuación invoucrará los calculos a partior de V, las fuerzas fi, y torques tao, actuando sobre las moléculas. El hamiltoniano también dictará la funcion de distribución de equilibrio para las posiciones moleculares y los momentos. Entonces, generalmente, es H (o V) es la entrada básica de una simulación computacional. La aproximación usada universalmente en una simulación computacional es el particionar la energía potencial en términos que involucren pares, tripletes, etc., de moléculas. En la siguiente sección se considerará esto en detalle. Antes de abandonar esta sección, mencionaremos brevemente algunas diferentes aproximaciones de los calculos de V. En estos desarrollos, la distribución de electrones en el sistema no se modela por un potencial efectivo V(q), pero es tratado por una forma de la teoría de densidad funcional. En una aproximación, la densidad electrónica se representa por una extensión de la teorúia del gas de electrones. En otra, los grados de libertad electrónica son explicitamente incluídos en la descripción, y los electrones son aceptados a relñajarse durante el curso de la simulación por el proceso conocido como simulación annealing (simulated annealing). Estos dos métodos permiten la división de V into pairwise and terminos de grado mayor. Esto parece prometedor para simulaciones futuras de sólidos y líquidos. SISTEMAS ATÓMICOS. Cnsideremos primero el caso de un sistema que contiene N atomos. La energía potencial puede ser dividida entre los terminos dependiendo de sus coordenadas de átomos individuales, pares, tripletes, etc.: El primer térimino en la ecuación 1.4 representa el efecto de un campo externo, (incluyendo, por ejemplo, las paredes contenedoras) en el sistema. Los terminos restantes representan las interacciones de partículas. El segundo término v2, el potencial par, es el más importante,. El potencial par depende solamente de la magnitud de la separación par, entonces puede ser escrito v2(rij). La figura 1.3 muestra uno de los más recientes estimados para el potencial par entre dos átomos de argon, como una función de la separación. Este potencial BBMS fue derivado considerando una gran cantidad de datos experimentales, incluyendo esparcimiento molecular de rayos, espectroscopía del argón (dimer), y propiedades de estado sólido, junto con los cálculos teóricos de las contribuciones de largo-rango. El potencial es también consistente con los estimados actuales de coeeficientes de transportación en la fase de gas. El potencial BBMS muestra las características típicas de la intercción intermolecular. Se muestra un comportamiento atractivo en grandes separaciones, esecnialmente debido a la correlación entre las nubes de electrones de los atoms (van der walls o london dispersion). En adición a especies de carga, terminos coulombianos estaran presentes. Hay un poso negativo, responsable de la cohesión en fases condensadas. Finalmente, hay un muro repulsivo, ó escalón repulsivo a distancias cortas, debido a la repulsión entre las nubes de neutrones. El término v3 en la ecuación (1.4), involucra tripletes de moléculas, es indudablemente significativo para densidades de líquidos. Estimaciones de la magnitud de la contribución del triple-dipole, y del sistema de tres cuerpos, han sido hechos para gases inertes en su estado sólido. Se ha encontrado que mas del 10 % de la energía (lattice) de argón, (y más en el caso de especies más polarizables), puede ser debido a estos términos no aditivos en el potencial, podemos esperar el mismo orden de magnitud para mantenerse en la fase líquida. Cuatro-cuerpos (y más de cuatro) terminos en la ecuación (1.4) se espera que sean pequeños en comparación con v2 y v3. Debido al tamaño de los terminos debido a los tres cuerpos en el potencial, son raramente incluidos en las simulaciones computacionales, esto es debido a que los calculos sobre la resolución de este problema tomaríoa mucho tiempo para su resolució. Afortunadamente, la aproximación de la interacción par da una buena descripción de las propiedades de los líquidos debido a que el promedio de los efectos de la interacción de tres cuerpos puede ser incluida parcialmente definiendo un par potencial "efectivo". Para hacer esto, reescribimos la ecuación (1.4) en la forma: Los potenciales pares que aparecen en las simulaciones computacionales son generalmente pares potenciales efectivos de este tiopo (regarded), representando todos los efectos de muchos cuerpos, para simplificar, solamente usaremos la notación v(rij) o v(r). A consecuencia de esta aproximación tenemos que el potencial par efectivo necesario para reproducir experimentalmente los datos puede apagarse para depender de la densidad, temperatura, etc, cuando el potencial de dos cuerpos v2(rij) claramente no lo hace. Ahora veamos el más simpple, más idealizado, potencial par comunmente usado en las simualciones computacionales. Esto refleja las características de las interacciones reales de una manera empirica y general. Ilustrado con el BBMS potencial de argon en la figura 3 es un simple potencial Lennard Jones. El cual provee una descripción razonable de las propiedades del argón, via simulación computacional, si los parámetros epsilon y sigma se escogen apropiadamente. El potencial tiene un largo rango de atractive tail de la forma, un pozo negativo de ondo, y un muro negativo (negative well) de ancho, y un muro repulsivo escalonado para distancias menores a . Para nuestro proposto de investigar las propiedades generales de los líquidos, y por comparación con la teoría, alatamente idelizados potenciales parespuede ser de valor. En la figura 1.4, ilustramos tres formas donde, algo irrealistics, son muy simples y convenientes para usar en las simulaciones computacionalesy en la teoríoa de estado líquido. Estas son: potencial de esfera-dura: El potencial de pozo cuadrado y el potencial de esfera suave donde nu es un parámetro, algunas veces entero. El potencial de esfera suave, se vuelve progresivamente suave conforme se va incrementando nu. El potencial de esfera suave contiene una parte no atractiva. Algunas veces e útil dividir los potenciales realmente en sus componentes atractivas y repulsivas, y la separación propuesta por Weeks, Chandler, y Anderson involucra splitting el potencial al mínimo. Para el potencial Lennard-Jones, las partes atractivas y repulsivas son: Esta separación se ilustra en la figura 1.5. En la teoría de perturbación, un fluido hipotético de moleculas interactúa via el potencial repulsivo , es tratado como un sistema de referencia y la parte atractiva es la perturbación. Se deberíua de notar que el potencial es significativamente fuerte en relación al inverso de la 12va potencia del potencial de esdera suave, que es pensado algunas veces como la parte repulsiva de Fig. 1.5 Para iones, como es de esperarse, estos potenciales no son suficiente para representar las interacciones de largo alcance. Una simple aproximación es para suplir uno de los potenciales pares con la interacción carga a carga de Coulomb. (1.11) donde zi, zj, son las cragas en los iones i y j, y epsilon es la permitividad del vacío, no confundirse con la otra epsilon. Para sistemas ionicos, las interacciones de inducción son importantes: la carga del ión induce un dipolo en la vecindad del ion. Ete termino no es aditivo par y es dificil incluirlo en la simulación. El modelo de concha es un intento cruso de tomar esta polarización del ión en cuenta. CALCULANDO EL POTENCIAL Este es el momento apropiado para introducir nuestra primera pieza de codigo computacional, el cual ilustra el cálculo de la energía potencial en un sistema de átomos de Lennard-Jones. Convirtiendo las ecuaciones algebraicas de este capítulo en código de FORTRAN. Suponemos que los vectores coordenados de nuestros átomso estaan situados en tres arreglos de Fortran RX(I), RY(I), RZ(I), con el índice de partícula I vaiando de 1 a N (el núemro de partículas). Para el potencial de Lennard-Jones es útil tener ya precomputarizado el valor de sigma cuadrada, el cual se encuentra en la variable SIGSQ. La energía potencial se pondrá en una variable llamada V, el cual se lleva a cero inicialmente, y posteriormente es acumulada en un doble loop sobre todos los distintos pares de átomos, teniendo cuidado de contar cada par solamente una vez: V = 0.0 DO 100 I=1, N-1 RXI =RX(I) RYI =RY(I) RZI =RZ(I) DO 99 J=I+1, N RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RIJSQ = RXIJ ** 2 + RYIJ ** 2 + RZIJ ** 2 SR2 = SIGSQ/RIJSQ SR6 = SR2*SR2*SR2 SR12 = SR6**2 V = V + SR12-SR6 CONTINUE CONTINUE V = 4.0 * EPSILON *V |
Attachment: maestria2nis.txt (0 KB) [
Download ]
|
Hay varias maneras de
codificar las correciones de mínima imagen. Hemos
mencionado también la posibilidad de conducir simulaciones en condiciones de contorno periódico no-cúbico. El código FORTRAN para implementar la correción mínima de imagen en el caso del octaedro truncado y el dodecaedro rómbico |
Attachment: buenoinpMC.txt (2 KB) [
Download ]
|
enter file
final.dat file final.dat Enter random number integer 2937484 iimc version 0.1 Montecarlo simulation of Lennard-Jones atoms Periodic Boundaries in all directions HECTOR DOMINGUEZ Enter file to keep last running data finalc1.dat file finalc1.dat Enter file to keep g(r) data grfilc.dat file grfilc.dat Enter max. displ. par. 0.3 Enter number of blocks 10 Enter number of cycles 500 Enter interval for update of max. displ. 10 Enter temperature 0.7 Enter potential cutoff distance 3.5 No. points in distr.func. = 100 Thickness for g(r) bins = 0.05 Number of blocks = 10 Number of cycles = 500 Ratio update interval for atoms = 10 Maximum displacement = 0.3 initial temperature = 0.7 Potential cutoff distance = 3.5 -- Virial --= -970.579174 Initial values # Part. = 500 Potential = -3.6020609 Total eng = -2.5520609 Temperature = 0.7 Density = 0.5 Pressure = -0.620579174 start of Markov Chain cycle Potential pressure 1 -3.9087 -0.4861 2 -4.0612 -0.4100 3 -4.1057 -0.4241 4 -4.2001 -0.3787 5 -4.1569 -0.3845 6 -4.1373 -0.3958 7 -4.1988 -0.3953 8 -4.2517 -0.3075 9 -4.2680 -0.4024 10 -4.2723 -0.3557 ** run values ** -4.1561 -0.3940 ** fluctation values ** 0.1215 0.0987 Final values # Particles = 500 Potential = -4.27195038 Total eng = -3.22195038 Temperature = 0.7 pressure = -0.270289173 [roberto@hdom01 MAESTRIA]$ |
Attachment: tesis1.txt (12 KB) [
Download ]
|
LA APROXIMACIÒN SEMICLASICA Un tratamiento completo de mecànica cuàntica sobre el esparcimiento elàstico desde un potencial central W(R) inicia desde la ecuaciòn independiente del tiempo de Schrödinger (en el sistema coordenado del centro de masa) (2.33) el cual tiene que resolverse sujeto a la condiciòn de frontera (2.34) donde E = . El primer término en (2.34) representa una onda plana incidente paralela al eje Z y el segundo representa las ondas esparcidas hacia afuera. La amplitud de la onda esparcida, , determina la secciòn de cruce diferencial dentro de la relaciòn En el tratamiento de onda-parcial (partial-wave) de Faxén y Holtsmark (1927), es desarrolladsa en series de polinomios de Legendre. Se asume que W(R) no es más singular en el origen de lo que es * y no es significativo, en *, más rápido que 1/R. Si Pl(x) es el polinomio de Legendre de oren l, tenemos: (2.36) donde las funciones radiales Ul(R) son las soluciones regulares de: (2.37) Estas soluciones tienen la forma asintótica a lo largo de R: (2.38) donde Cl y * son independientes de R. Las constantes * son conocidas como las fases superiores (phase shifts), y determinan la amplitud de esparcimiento, las cuales tienen una expansión en serie como se muestra a cotinuación: (2.39) La expresión (2.39) para f(*) siempre es convergente para el potencial central satisfaciendo las condiciones específicas, pero la evaluación de la suma es algunas veces impractica, o inconveniente, en el esparcimiento átomo-átomo debido a la gran cantidad de valores de importancia asociados a l, (o de l). Debido a esta dificultad, métodos semi-clásicos han sido desarrollados, los cuales tienen mucha de la simplicidad de la aproximación puramente clásica, pero los cuales permiten para los efectos de interferencia cuántica, los cuales son producidos en el esparcimiento por potenciales reales. En la aproximaciòn semi-clàsica (Ford y Wheeler 1959; Berry y Mount 1972; Child 1974), los términos principales son retenidos en la expansión de el levantamiento de fase * en potencias de h, usando el método WKB (Wentzel-Kramers-Brillouin). Un potencial efectivo de una dimensión Weff(R) puede ser definido como (2.40) donde el segundo término refpreseta la barrera centrífuga. Proveído de este potencial posee un punto clásico de vuelta, la siguiente WKB aproximación se obtiene del levantamiento de fase. (2.41) donde (2.42) Para una atracción fuerte, la variación de Weff con l se muestra en la figura 2.5. Parece mostrar que para valores grandes de E y l existe un punto de vuelta clásico, pero para pequeños valores existen tres puntos de vuelta en R=Ro, R1, y R2. La orbita singular clásuica ocurre cuando la energía E coincide con el máximo en el potencial efectivo donde la puntos de vuelta en Ro y R1 se juntan. (Fig. 5(b)). Fig 5(a) La variación del potencial efectivo Weff(R) definiido por (2.40) para valores diferentes del momento angular cuántico l, con l1>l2>l3>l4>l5. (b) Esta es una ilustración a gran escala. Es una escala mayor a (a) y muestra la variación de Weff(R) con R para un valor dado de l= l', donde l' es suficientemente pequeño para Weff(R) para exibir un máximo y un mínimo. Para energías grandes tal como E = E3 existe un punto de vuelta clásico en R = Ro', pero para energías pequeñas tales como E = E1 existen tres puntos de vuelta en Ro, R1, y R2. Para la energía E = E2 la cual coincide con el máximo en Weff(R) los puntos de vuelta (quiebre) Ro y R1 se juntan. La modificación a la fase semi-clásica de cambio para tres puntos de vuelta se discutirá posteriormente. La condición para la expresión (2.41) parta ser una buena aproximaión para la fase de cambio exácta es que el potencial efectivo debería no variar apreciablemente sobre distancias comparables a la longitud de onda de de Broglie, excepto en la vecidad inmediata del punto de vuelta clásico. esta condición puede ser expresada como sigue: (2.43) En general, bajo condiciones semi-clásicas, los valores importantes de l son grandes, en este caso l puede ser relacionada al parámetro de impacto clásico b por la siguiente relación: (2.44) donde L es el momento angulare clásico. Haciendo esta identificación y tratando a l como una variable continua parece ser que k(R) y f(R), definidas en (2.42) y (2.24), respectivamente, difieren sólo por una constante: y también, de (2.26), (2.25) y (2.41) (2.45) Esto es conocido como la relación clásica de equivalencia. El desarrollo de la aproximación semi-clásica procede de dos pasos nuevos: 1. Asumir en la expresión (2.39) para el esparcimiento de amplitud f(*) se puede reemplazar por una integración, de la siguiente manera: (2.46) 2. Los polinomios de Legendre son reemplazados por el término principal de sus expansiones asintóticas: (2.47a) (2.47b) donde Jo(x) es la función de Bessel de primer tipo de orden cero. Estas aproximaciones para Pl(cos *) pueden ser empleadas para toda l>1. Usando las aproximaciones (2.47a y b) en (2.46), la integral sobre l puede ser evaluada por el método de fase estacionaria o una aproximación similar. 2.5 ESPARCIMIENTO HACIA FUERA DE LA REGIÓN DEL ANGULO PEQUEÑO (SMALL-ANGLE). Debemos considerar primero el esparcimiento en ángulos cuya aproximación (2.47a) se mantenga, esto es, a ángulos no muy cerca de la dirección horizontal (*=0) o (*=*). Usando la relación para *=0 (2.48) se muestra que el término [*] en (2.46) puede reemplazarse por exp(*). Usando (2.47a), econtramos entonces: (2.49a) donde (2.49b) Para la mayoría de valores de l, la variación de * y * con l es más rápida y las partes reales e imaginarias del integrando en (2.49a) oscila rapidamente en relación al cero. Las mayores contribuciones a la integral entonces vienen de regiones en donde tanto * ó * es cero, y estan son reguiones de la fase estacionaria. Usando la relación clásica de equivalencia (2.45) vemos que si * =0, para l =l'entonces: (2.50) Esta relación muestra que el punto de fase estacionaria, en l = l', para un ángulo dado de esparcimientro *, se determina por el parámetro de impacto clásico correspondiente al mismo ángulo de esparcimiento. La relación (2.45) y (2.50) son correctas en el caso que exista solamente un punto de vuelta clásico en el potencial efectivo EWeff(R), en tal caso |*| < . Si el orbitamiento ocurre para el cual |*| > *, una modificación debe llevarse a cabo, tal consideración debe tomar en cuenta lo que se muestra en a sección 2.7. Si *(b,E) se alcanza para una trayectoria repulsiva, *>0 y el punto de fase estacionaria está en *, cuando *(b,E) se levanta de una trayectoria atractiva, *<0, y el punto de fase estacionaria está en *. POTENCIALES REPULSIVOS Vamos a considerar primero un potencial repulsivo W(R)>0 , para el cual *(b,E) decrece monótonamente desde *=* a b=0 a b - *. Hay un punto de fase estacionaria para un ángulo dado de esparcimiento * en l = l', esto es en b=b', y expandiendo * con referencia a este punto encontramos (2.51) Desde (2.52) vemos que * es negativo y (2.53) Usando (2.53) y (49), se encuentra que: (2.54) donde, desde la única contribución importante a la integral viene de la reguión de l = l', el término de variación lenta (l + 1/2)^(1/2) se ha tomado fuera de la iontegral, y el límite inferior se ha extendido desde l=0 a l = -*. Usando el resultado estándar (2.55) obtenemos: (56) donde (57) y donde I(*) es la sección de cruce diferencial clásica (30). Se sigue que para un potencial repulsivo puro la semi-clásica y clásica sección de cruce son identicas. Este resultado no es correcto si sólo unos cuantos cambios de fase de bajo orden contribuyen a la suma (39) y requiere que la región de fase estacionaria , de ancho *l =**, es grande en comparación con la unidad. Es claramente necesario que el potencial W(R) debe variar ligeramente con R en orden de que * debe también variar ligeramente (smoothly) con l. REPULSIÓN INTERNA- ATRACCIÓN EXTERNA Vamos a considerar la forma típica de las funciones de deflección para esparcimiento ión-átomo, en la ausencia de orbitación, mostrado en la Fig. 4, correspondiente al potencial W(R) con una región repulsiva interna y una región xterna atractiva. Parece mostrar que el valor más negativo de * es igual a -*, donde * es el ángulo arcoiris y si Provee la región de fase estacionaria en b = bi están bien separados de la amplitud de esparcimiento f(*) es la suima de las contribuciones de fi(*) para cada región, entonces: (2.58) La contribución de los tres puntos de la fase estacionaria, para *<*. puede ser evluada a lo largo de las mismas líneas para el potencial repulsivo. Se encuentra que (2.59) donde * es la sección de cruce diferencial para b = bi y (60) Para obtener (60) el hecho es que se uso que * es negativa en l = l1 y l =l2 (esto es b = b1 y b =b2) pero positiva en l = l3 (b = b3) Fig 6. Forma típica de una sección de cruce diferencial, ampliada por sin *, para una fucnión de deflección de la forma mostrada en la Fig. 4. El primario y supernumerable arcoiris máximo son A, B, C, ... y el ángulo arcoiris en *=* se marca por una flecha. Las oscilaciones rápidas son debidas a la interferencia entre la amplitud para el esparcimiento arcoiris f(*) con la amplitud fo(*), *>* o f1(*), *<* (ver (21) y (63)). La sección de cruce diferencial (2.61) muestra dos tipos de oscilaciones, como se ilustra en la Fig. 6. Existen oscilaciones lentas en I(*) como una función de *, debidas a la interferencia entre f2(*) y f3(*), controlada por la diferencia de fase *0 constante , entonces la separación angular entre los máximos sucesivos es (2.62) La interferencia de f1 con f2 y f1 con f3, superimpone mucho más rápidas oscilaciones en este patrón de arcoiris supernumerable. La vecindad del ángulo arcoiris *=* requiere un tratamiento especial, como los puntos de fase estacionaria en l2 y l3 se juntan y f2(*) y f3(*) no pueden ser tratados separadamente. Discutiremos esto en el siguiente párrafo. Como vimos la amplitud del esparcimiento arcoiris, fm(*), se extiende en la región *>*, entonces debido a esto también hay solamente un punto de fase estacionaria en esa región, en l = lo correspondiente a la amplitud fo(*), la sección de cruce diferenial (2.63) también muestra oscilaciones de interferencia, tiñendo fuera conforme * se incrementa. 6 ESPARCIMIENTO ARCOIRIS En e punto l = lm(b = bm), la función de deflexión clásica está en un mínimo y puede ser expresda de la siguiente manera: (64) Usando la relación de equivalencia (45) vemos que para l cercano a lm (65) de donde se encuentra que (con * = -*) (66a) donde (66b) Entonces de (49) (67a) con (67b) La integral en (67a) puede expresarse en términos de la función de Airy Ai(x), donde (68) Esta función decae exponencialmente para x>0 y oscila para x<0, con la forma asintótica (69a) (69b) desde (67a) y (68) (70) Para *>*m, la interferencia entre fm(*) y fo(*) produce oscilaciones, la amplitud del cual decrece exponencialmente con el incremento del ángulo (ver Fig. 6). Para *<*m la amplitud fm(*) oscila y reproduce el arcoiris supernumerario (supernumerary rainbow), el cual puede ser descrito alternativamente debido a la interacción entre f2(*) y f3(*). De hecho una aproximación sin provar para fm(*) debida a Berry(1966) reproduce exactamente el resultado fm(*) = f2(*) + f3(*), para valores de * menores que *m 7 ORBITACIÓN. |
Attachment: trabajo1.txt (1 KB) [
Download ]
|
DINÁMICA MOLECULAR ECUACIÓN DE MOVIMIENTO PARA SISTEMAS ATÓMICOS En el capítulo se trata con las técnicas usadas para resolver las ecuaciones clásicas de moviiento para un sistema de N moléculas interactuando vía un potencial V como el de la ecuacion (1.4). Estas ecuaciones pueden escribirse de diferentes maneras. La manera más apropiada es la ecuación de movimiento de Lagrange L = H-V Con la función de Lagrange definida como L(q,dq/dt), se define en terminos de la energia potencial y cinética, dependiente de coordenadas generalizadas q Si consideramos un sistema de átomos, con coordenadas cartesianas ri y las definiciones usuales de H y V, entonces rtenemos: mi dri/dt(dri/dt) = fi. donde mi, es la mas a del átomo i y fi = dL/dt es la fuerza en el átomo. Esta ecuaciones también se aplican al marco de centro de masa, con fi representando la fuerza total sobre la molécula i. |
Attachment: MC segunda prueba.txt (2 KB) [
Download ]
|
[roberto@hdom01
MAESTRIA]$ ./canonico
enter file final.dat file final.dat Enter random number integer 2937484 iimc version 0.1 Montecarlo simulation of Lennard-Jones atoms Periodic Boundaries in all directions HECTOR DOMINGUEZ Enter file to keep last running data final2c.dat file final2c.dat Enter file to keep g(r) data grfilc2.dat file grfilc2.dat Enter max. displ. par. 0.3 Enter number of blocks 10 Enter number of cycles 500 Enter interval for update of max. displ. 10 Enter temperature Interrupt [roberto@hdom01 MAESTRIA]$ ./canonico enter file final.dat file final.dat Enter random number integer 2937484 iimc version 0.1 Montecarlo simulation of Lennard-Jones atoms Periodic Boundaries in all directions HECTOR DOMINGUEZ Enter file to keep last running data finalc2.dat file finalc2.dat Enter file to keep g(r) data grfilc2.dat file grfilc2.dat Enter max. displ. par. 0.3 Enter number of blocks 100 Enter number of cycles 5000 Enter interval for update of max. displ. 10 Enter temperature 0.7 Enter potential cutoff distance 3.5 No. points in distr.func. = 100 Thickness for g(r) bins = 0.05 Number of blocks = 100 Number of cycles = 5000 Ratio update interval for atoms = 10 Maximum displacement = 0.3 initial temperature = 0.7 Potential cutoff distance = 3.5 -- Virial --= -970.579174 Initial values # Part. = 500 Potential = -3.6020609 Total eng = -2.5520609 Temperature = 0.7 Density = 0.5 Pressure = -0.620579174 start of Markov Chain cycle Potential pressure 1 -4.1561 -0.3940 2 -4.4919 -0.2639 |
Attachment: buenoresult.txt (4 KB) [
Download ]
|
[roberto@hdom01
MAESTRIA]$ ./coulomb
enter file titoe.dat file titoe.dat Enter templating option 1) Fluid 1 with non-zero velocities 2) Fluid 1 with zero velocities (frozen molecules) 1 ------------- Couloumb energy= -31.353118kcal/mol Couloumb energy= -0.213479327/KT vad de waal en.= 2.16039057kcal/mol vad de waal en.= 0.0147098202/KT md version 0.1 Molecular Dynamics of soft atoms with electrostatic interactions. If Templating is on, then the first n1 particles are frozen. Periodic Boundaries in all directions Option of having diffusion coefficients HECTOR DOMINGUEZ --- March 2002 --- Enter file to keep last running data titog.dat file titog.dat Enter file to keep g(r), pair correl. func. data titog1.dat file titog1.dat Enter number of blocks 10 Enter number of timesteps per block 20 Enter timestep (ps) .002 Enter constant temperature option .true. Enter required temperature (Kelvin) 303 Enter potential cutoff distance (A) 7.8750E+01 Enter LJ epsilon(J/KT), sigma(A) for 1-1 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 1-2 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 2-2 fluid inter. 0.49207,3.743 Enter mass 1 and mass 2 (uma) 16.04,16.04 Enter lb factor for electrostatic (e*e/4*pi*eps*k*T) 10.350825 Enter kxm,kym,kzm in ewald sum and k-factor 4,4,4,2.5571E-02 Do you want diffusion 1) yes, 2) no 2 point in g(r) distr. func = 100. Number of blocks = 10 Number of cycles = 20 Timestep = 0.002 Constant temperature option = T initial temperature = 303. Potential cutoff distance = 78.75 Epsilon and sigma for the 1-1 = 0.49207 3.743 Epsilon and sigma for the 1-2 = 0.49207 3.743 Epsilon and sigma for the 2-2 = 0.49207 3.743 Mass of 1 and Mass of 2 = 16.04 16.04 kxm,kym,kzm,kf. for Ewald = 4 4 4 0.025571 --- MIXTURE OF + and - CHARGES --- Initial values N1 = 124 N2 = 124 Density of 1= 3.17359417E-05(N/A3) Density of 2= 3.17359417E-05(N/A3) Energy = 1.48225397(J/KT) Kinetic = 1.5(J/KT) Potential = -0.0177460296(J/KT) coulumbic = -0.350127301(J/KT) Temperature = 298.735037(K) Pressure = 0.00227924356(Kbar) Volume = 3907241.87(A3) Block Kinetic Coulomb Potential Momentum Temp Press J/NKT J/NKT J/NKT Kg-m/Ns K Kbar 1 1.5000 -0.3484 -0.0171 0.0040 303.0000 0.0023 2 1.5000 -0.3457 -0.0164 0.0040 303.0000 0.0023 3 1.5000 -0.3441 -0.0162 0.0040 303.0000 0.0023 4 1.5000 -0.3432 -0.0165 0.0040 303.0000 0.0023 5 1.5000 -0.3431 -0.0170 0.0040 303.0000 0.0023 6 1.5000 -0.3436 -0.0175 0.0040 303.0000 0.0023 7 1.5000 -0.3448 -0.0179 0.0040 303.0000 0.0023 8 1.5000 -0.3464 -0.0176 0.0040 303.0000 0.0024 9 1.5000 -0.3482 -0.0166 0.0040 303.0000 0.0024 10 1.5000 -0.3499 -0.0157 0.0040 303.0000 0.0024 ------------------------------------------------------- 1.5000 -0.3457 -0.0169 0.0040 303.0000 0.0023 ------------------------------------------------------- ** fluctation values ** 0.0000 0.0007 0.0000 303.0000 0.0000 Final values # Part. fl-1= 124 # Part. fl-2= 124 Potential = -0.0159187354 coulumbic = -0.350594784 Total eng = 1.48408126 Temperature = 303. Final rcut = 78.75 Density fl-1= 3.17359417E-05 Density fl-2= 3.17359417E-05 pressure = 0.00243153334 Real thickness in distr. func. = 0.7875173 |
Attachment: manual canonico.txt (0 KB) [
Download ]
|
INBSTRUCCIONES PARA INSERTAR EN CANONICO ./canonico final.dat interger numeros de 8 a 9 dígitos -out1.dat -g(r) función de correlación -parámetro de máximo de numero -0.3 buen número - númnero de bloques 50 000 a 100 000 ? apoco, con 500 esta bien -número de ciclos 1000 ciclos intervalo- 10 En temperatura absoluta - mayor a 1.2 distancia al potencial : 3.5 a 4 tener cuidado que el radio de corte no sea más grande que la mitad de la longitud de la caja |
Attachment: manualcoulomb.txt (2 KB) [
Download ]
|
[roberto@hdom01
MAESTRIA]$ ./coulomb
enter file titoe.dat file titoe.dat Enter templating option 1) Fluid 1 with non-zero velocities 2) Fluid 1 with zero velocities (frozen molecules) 1 ------------- Couloumb energy= -31.353118kcal/mol Couloumb energy= -0.213479327/KT vad de waal en.= 2.16039057kcal/mol vad de waal en.= 0.0147098202/KT md version 0.1 Molecular Dynamics of soft atoms with electrostatic interactions. If Templating is on, then the first n1 particles are frozen. Periodic Boundaries in all directions Option of having diffusion coefficients HECTOR DOMINGUEZ --- March 2002 --- Enter file to keep last running data titog.dat file titog.dat Enter file to keep g(r), pair correl. func. data titog1.dat file titog1.dat Enter number of blocks 10 Enter number of timesteps per block 20 Enter timestep (ps) .002 Enter constant temperature option .true. Enter required temperature (Kelvin) 303 Enter potential cutoff distance (A) 7.8|750E+01 Enter LJ epsilon(J/KT), sigma(A) for 1-1 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 1-2 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 2-2 fluid inter. 0.49207,3.743 Enter mass 1 and mass 2 (uma) 16.04/16.04 Enter lb factor for electrostatic (e*e/4*pi*eps*k*T) 10.350825 Enter kxm,kym,kzm in ewald sum and k-factor 4,4,4,2.5571E-02 Do you want diffusion 1) yes, 2) no 2 point in g(r) distr. func = 100. Number of blocks = 10 Number of cycles = 20 Timestep = 0.002 Constant temperature option = T initial temperature = 303. Potential cutoff distance = 7.8 Epsilon and sigma for the 1-1 = 0.49207 3.743 Epsilon and sigma for the 1-2 = 0.49207 3.743 Epsilon and sigma for the 2-2 = 0.49207 3.743 Mass of 1 and Mass of 2 = 16.04 0. kxm,kym,kzm,kf. for Ewald = 4 4 4 0.025571 --- MIXTURE OF + and - CHARGES --- Initial values N1 = 124 N2 = 124 Density of 1= 3.17359417E-05(N/A3) Density of 2= 3.17359417E-05(N/A3) Energy = NAN(J/KT) Kinetic = NAN(J/KT) Potential = -0.015936583(J/KT) coulumbic = -0.24717578(J/KT) Temperature = NAN(K) Pressure = NAN(Kbar) Volume = 3907241.87(A3) Block Kinetic Coulomb Potential Momentum Temp Press J/NKT J/NKT J/NKT Kg-m/Ns K Kbar |
Attachment: prueba1.txt (2 KB) [
Download ]
|
Interrupt
[roberto@hdom01 MAESTRIA]$ ./canonico enter file final.dat file final.dat Enter random number integer 7 iimc version 0.1 Montecarlo simulation of Lennard-Jones atoms Periodic Boundaries in all directions HECTOR DOMINGUEZ Enter file to keep last running data out12.dat file out12.dat Enter file to keep g(r) data out13.dat file out13.dat Enter max. displ. par. 0.3 Enter number of blocks 10 Enter number of cycles 10 Enter interval for update of max. displ. 1 Enter temperature 1.2 Enter potential cutoff distance 3.5 No. points in distr.func. = 100 Thickness for g(r) bins = 0.05 Number of blocks = 10 Number of cycles = 10 Ratio update interval for atoms = 1 Maximum displacement = 0.3 initial temperature = 1.2 Potential cutoff distance = 3.5 -- Virial --= -970.579174 Initial values # Part. = 500 Potential = -3.6020609 Total eng = -1.8020609 Temperature = 1.2 Density = 0.5 Pressure = -0.370579174 start of Markov Chain cycle Potential pressure 1 -3.4528 -0.0164 2 -3.3477 0.1299 3 -3.3400 0.1890 4 -3.3099 0.1418 5 -3.3188 0.1437 6 -3.2694 0.2797 7 -3.3266 0.0817 8 -3.3132 0.1270 9 -3.2858 0.1516 10 -3.3025 0.1775 ** run values ** -3.3267 0.1405 ** fluctation values ** 0.0563 0.1131 Final values # Particles = 500 Potential = -3.30932227 Total eng = -1.50932227 Temperature = 1.2 pressure = 0.276037918 [roberto@hdom01 MAESTRIA]$ |
Attachment: prueba1DM.txt (4 KB) [
Download ]
|
[roberto@hdom01
MAESTRIA]$ ./coulomb
enter file titoe.dat file titoe.dat Enter templating option 1) Fluid 1 with non-zero velocities 2) Fluid 1 with zero velocities (frozen molecules) 1 ------------- Couloumb energy= -31.353118kcal/mol Couloumb energy= -0.213479327/KT vad de waal en.= 2.16039057kcal/mol vad de waal en.= 0.0147098202/KT md version 0.1 Molecular Dynamics of soft atoms with electrostatic interactions. If Templating is on, then the first n1 particles are frozen. Periodic Boundaries in all directions Option of having diffusion coefficients HECTOR DOMINGUEZ --- March 2002 --- Enter file to keep last running data titog.dat file titog.dat Enter file to keep g(r), pair correl. func. data titog1.dat file titog1.dat Enter number of blocks 1000 Enter number of timesteps per block 20 Enter timestep (ps) .002 Enter constant temperature option .true. Enter required temperature (Kelvin) 303 Enter potential cutoff distance (A) 7.8750E+01 Enter LJ epsilon(J/KT), sigma(A) for 1-1 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 1-2 fluid inter. 0.49207,3.743 Enter LJ epsilon(J/KT), sigma(A) for 2-2 fluid inter. 0.49207,3.743 Enter mass 1 and mass 2 (uma) 16.04,16.04 Enter lb factor for electrostatic (e*e/4*pi*eps*k*T) 10.350825 Enter kxm,kym,kzm in ewald sum and k-factor 4,4,4,2.5571E-02 Do you want diffusion 1) yes, 2) no 2 point in g(r) distr. func = 100. Number of blocks = 1000 Number of cycles = 20 Timestep = 0.002 Constant temperature option = T initial temperature = 303. Potential cutoff distance = 78.75 Epsilon and sigma for the 1-1 = 0.49207 3.743 Epsilon and sigma for the 1-2 = 0.49207 3.743 Epsilon and sigma for the 2-2 = 0.49207 3.743 Mass of 1 and Mass of 2 = 16.04 16.04 kxm,kym,kzm,kf. for Ewald = 4 4 4 0.025571 --- MIXTURE OF + and - CHARGES --- Initial values N1 = 124 N2 = 124 Density of 1= 3.17359417E-05(N/A3) Density of 2= 3.17359417E-05(N/A3) Energy = 1.48225397(J/KT) Kinetic = 1.5(J/KT) Potential = -0.0177460296(J/KT) coulumbic = -0.350127301(J/KT) Temperature = 298.735037(K) Pressure = 0.00227924356(Kbar) Volume = 3907241.87(A3) Block Kinetic Coulomb Potential Momentum Temp Press J/NKT J/NKT J/NKT Kg-m/Ns K Kbar 1 1.5000 -0.3484 -0.0171 0.0040 303.0000 0.0023 2 1.5000 -0.3457 -0.0164 0.0040 303.0000 0.0023 RESULTADOS 996 1.5000 -0.3825 -0.0254 0.0040 303.0000 0.0024 997 1.5000 -0.3806 -0.0240 0.0041 303.0000 0.0024 998 1.5000 -0.3789 -0.0217 0.0041 303.0000 0.0024 999 1.5000 -0.3772 -0.0242 0.0041 303.0000 0.0024 1000 1.5000 -0.3757 -0.0275 0.0041 303.0000 0.0023 ------------------------------------------------------- 1.5000 -0.3963 -0.0231 0.0040 303.0000 0.0024 ------------------------------------------------------- ** fluctation values ** 0.0000 0.0058 0.0000 303.0000 0.0000 Final values # Part. fl-1= 124 # Part. fl-2= 124 Potential = -0.0280749931 coulumbic = -0.375318767 Total eng = 1.47192501 Temperature = 303. Final rcut = 78.75 Density fl-1= 3.17359417E-05 Density fl-2= 3.17359417E-05 pressure = 0.00229995958 Real thickness in distr. func. = 0.7875173 enter file |
Attachment: prueba1MC.txt (1 KB) [
Download ]
|
[roberto@hdom01
MAESTRIA]$ ./canonico
enter file final.dat file final.dat Enter random number integer 765342739 iimc version 0.1 Montecarlo simulation of Lennard-Jones atoms Periodic Boundaries in all directions HECTOR DOMINGUEZ Enter file to keep last running data out1.dat file out1.dat Enter file to keep g(r) data out2.dat file out2.dat Enter max. displ. par. 0.3 Enter number of blocks 1000 Enter number of cycles 20 Enter interval for update of max. displ. 10 Enter temperature 303 Enter potential cutoff distance 7.8750E+01 No. points in distr.func. = 100 Thickness for g(r) bins = 0.05 Number of blocks = 1000 Number of cycles = 20 Ratio update interval for atoms = 10 Maximum displacement = 0.3 initial temperature = 303. Potential cutoff distance = 78.75 -- Virial --= -1047.60696 Initial values # Part. = 500 Potential = -3.77735828 Total eng = 450.722642 Temperature = 303. Density = 0.5 Pressure = 150.452393 start of Markov Chain cycle Potential pressure 1 28.2289 221.0384 resultados Final values # Particles = 500 Potential = 28.8529239 Total eng = 483.352924 Temperature = 303. pressure = 222.019915 [roberto@hdom01 MAESTRIA]$ |
