Elementos De Simulacion Computacional.pdf

ELEMENTOS DE SIMULACION COMPUTACIONAL Din´ amica Molecular y M´ etodo de Monte Carlo

Gonzalo Guti´errez Departamento de F´ısica, Universidad de Santiago de Chile. e–mail: [email protected]

Abril 2001 Santiago, Chile

Indice

i

´Indice 1. Introducci´ on 2. Elecci´ on del Modelo 2.1. Ensembles de la Simulaci´ on . 2.2. El potencial interat´ omico . . 2.2.1. La forma del potencial 2.3. Implementaci´ on . . . . . . . .

1

. . . .

1 2 2 3 3

3. Inicializaci´ on 3.1. Condiciones iniciales y condiciones de borde . . . . . . . . . . . . . . . . . . . . .

5 5

4. Generaci´ on de las configuraciones 4.1. Evaluaci´ on de la energ´ıa potencial y las fuerzas . . . . . . . . . . . . . . . . . . . 4.2. Din´ amica Molecular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 6 7 9

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5. An´ alisis de los resultados 11 5.1. Propiedades termodin´ amicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.2. Propiedades est´ aticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 5.3. Propiedades din´ amicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A. Ap´ endice 14 A.1. Comentario Bibliogr´ afico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.2. Direcciones en Internet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 A.3. Simulaci´ on computacional en Chile . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Lista de Abreviaciones MC : Monte Carlo DM : Din´ amica Molecular CBP: Condiciones de Borde Peri´odicas FDP: Funci´ on de Distribuci´ on de pares

Acerca de este documento Notas de clase con motivo del curso “Simulaci´on en Materia Condensada”, Depto. F´ısica, Usach, 2001.

ii c

Gonzalo Guti´errez, 2da. versi´ on, Abril 2001. Enviar cr´ıticas, comentarios, sugerencias a [email protected]

Introducci´on

1.

1

Introducci´ on

El desarrollo de los computadores digitales a partir de la d´ecada de los ’50, y su aplicaci´on a la resoluci´ on de problemas cient´ıficos, ha introducido lo que algunos han llamado “una tercera metodolog´ıa” a la investigaci´ on cient´ıfica: la simulaci´on computacional [1]. Este m´etodo, de car´acter complementario y muchas veces alternativo a los modos convencionales de hacer ciencia, el experimental y el te´ orico, ha ejercido un fuerte impacto en pr´acticamente todos los campos de la ciencia (ver por ejemplo [1, 2]). El objetivo de la simulaci´on computacional es resolver los modelos te´ oricos en su total complejidad, mediante la resoluci´on num´erica de las ecuaciones involucradas, haciendo uso intensivo (y extensivo) de computadores. En el ´ area de la f´ısica, la simulaci´on computacional fue introducida como una herramienta para tratar sistemas de muchos cuerpos a comienzo de la d´ecada de los ’50, con el trabajo pionero de Metropolis et al. [3]. M´ as tarde, auspiciosos resultados iniciales obtenidos en mec´anica estad´ıstica cl´ asica, en particular en el estudio de l´ıquidos, dieron credibilidad a la simulaci´on computacional, extendi´endose r´ apidamente su uso a temas tan diversos como cromodin´amica cu´antica, f´ısica de flu´ıdos, relatividad general, f´ısica del plasma, materia condensada, f´ısica nuclear y ciencia de materiales. Actualmente, gracias al vertiginoso desarrollo de la tecnolog´ıa de los computadores, cuya velocidad crece aproximadamente un factor 2 cada dieciocho meses, la simulaci´on computacional se ha constitu´ıdo en una herramienta de c´alculo esencial, tanto para experimentalistas como para te´ oricos. Mediante un buen modelo computacional no s´olo se pueden reproducir experimentos de laboratorio, sino que adem´ as, gracias a que se pueden variar libremente los par´ametros usados, permite probar (o falsar) modelos te´oricos existentes en rangos de par´ametros imposibles de alcanzar experimentalmente por ahora, resolviendo as´ı conflictos entre explicaci´on te´orica y observaci´ on. Un papel fundamental tambi´en lo juega hoy d´ıa la visualizaci´on de los resultados obtenidos. No s´ olo obtenemos datos num´ericos que pueden ser contrastados con los experimentos, sino tambi´en obtenemos una imagen gr´afica del proceso en cuesti´on. Los dos m´etodos de simulaci´ on computacional m´as usados en f´ısica actualmente son el de la Din´ amica Molecular [8, 9, 10, 11, 12], que es de car´acter determinista, y el de Montecarlo, que es de car´ acter probab´ılistico [6]. Ambos pueden considerarse como m´etodos para generar configuraciones diferentes de un sistema de part´ıculas1 , es decir puntos en el espacio de fases compatibles con las condiciones externas. El m´etodo de la Din´ amica Molecular y el de Montecarlo ha sido empleado con ´exito para simular gases, l´ıquidos y s´ olidos [14], ampli´andose tanto su uso como el desarrollo de t´ecnicas espec´ıficas en forma paralela al avance tecnol´ogico de los computadores. Los sistemas estudiados van desde cientos a miles y u ´ltimamente incluso a decenas de millones de ´atomos. Los aspectos estudiados incluyen propiedades estructurales, termodin´amicas, mec´anicas y cin´eticas. Cabe se˜ nalar que estas t´ecnicas de simulaci´on tienen aplicaciones mucho m´as amplias que las aqu´ı esbozadas. En particular el m´etodo de MC se emplea con ´exito en F´ısica de Part´ıculas (por ejemplo lattice gauge theory), as´ı como para calcular sistemas cu´anticos, donde se ocupan t´ecnicas tales como Path integral MC y otros [7].

2.

Elecci´ on del Modelo

El punto de partida para simular un sistema f´ısico es definir con claridad el problema en cuesti´ on: que tipo de propiedades nos interesa estudiar, dentro de que rango de par´ametros, 1 Por

part´ıculas entendemos ´ atomos o mol´ eculas

2

Elecci´ on del modelo

con que precisi´ on. En funci´ on de ello debemos decidir el n´ umero de part´ıculas a usar, cuales ser´an las variables de control, que potencial interat´omico usar, que tipo de promedios debemos calcular, en que ensemble conducir la simulaci´on.

2.1.

Ensembles de la Simulaci´ on

La informaci´ on que genera una corrida de Din´amica Molecular es la posici´on y la velocidad de cada part´ıcula del sistema en cada instante de tiempo. Por su parte, en el caso de MC lo que obtenemos es la posici´ on de las part´ıculas en cada paso de simulaci´on. Empleando las t´ecnicas tradicionales de la Mec´ anica Estad´ıstica podemos pasar de esta informaci´on microsc´opica a la obtenci´on de magnitudes macrosc´ opicas que nos permitan conectar con el experimento, v´ıa la termodin´ amica. Supongamos que estamos tratando un sistema puro compuesto por N part´ıculas, encerrado en un volumen V y con una energ´ıa fija E. Las posiciones y velocidades definen un espacio de fases de 6N dimensiones. Obtener la posici´on y la velocidad de cada una de las part´ıculas, en cada instante, significa obtener la trayectoria de un punto Γ del espacio de fase en funci´on del tiempo, esto es, Γ(t). Denotemos por A el valor inst´antaneo de un cierto observable. El promedio de esta cantidad A est´ a dado por hAiobs = hAitiempo =

τobs 1 X

τobs

A(Γ(τ )) ,

(1)

τ =1

donde τ representa un tiempo discreto (‘pasos’ de din´amica molecular) y τobs son los pasos totales de la corrida. Suponiendo que el sistema es erg´odico, podemos asociar directamente este promedio con el promedio usual sobre ensemble de la Mec´anica Estad´ıstica, hAiobs = hAitiempo = hAiens .

(2)

En otras palabras, por medio del formalismo de la simulaci´on lo que se hace es generar una sucesi´on de diferentes estados (puntos) del espacio de fases (que se suponen no–correlacionados), compatibles con las condiciones externas ((N, V, E) es este caso), sobre los cuales se toman los promedios. La elecci´ on del ensemble bajo el cual llevar a cabo la simulaci´on est´a dictada fundamentalmente por el tipo de problema a tratar. Los promedios estad´ısticos pueden llevar a peque˜ nas diferencias en los diferentes ensembles, pero ´estas desaparecen en el l´ımite termodin´amico, que se alcanza incluso con unos pocos cientos de part´ıculas [9, 13]. Sin embargo la elecci´on del ensemble si influye al momento de calcular las fluctuaciones cuadr´aticas medias de las magnitudes termodin´ amicas. Estas permiten calcular, por ejemplo, la capacidad cal´orica o el m´odulo de elasticidad. El ensemble convencional en Din´ amica Molecular es el ensemble microcan´onico, (N, V, E), mientras que en Montecarlo es el can´ onico, (N, V, T ), donde T es la temperatura. Tambi´en se han desarrollado t´ecnicas para llevar a cabo simulaciones de MC en otros ensembles, as´ı como simulaciones de DM en el ensemble can´ onico [16], en el ensemble isot´ermico–isob´arico (N, P, T ) [17], a presi´on y temperatura constante y en el ensemble isoent´alpico–isotensi´on (H, t, N ) [18], a tensi´on externa t constante, entre otros. Para una revisi´on, con la referencia a los trabajos originales, de los diferentes m´etodos existentes remitimos al lector a la Ref. [13].

2.2.

El potencial interat´ omico

Un punto de importancia central tanto en Din´amica Molecular como en Montecarlo es la elecci´on del potencial interat´ omico del sistema a simular. De la fidelidad con que ´este represente

El potencial interat´omico

3

las interacciones reales entre las part´ıculas depender´a la calidad de los resultados: la conclusi´on inmediata es que mientras m´ as detalles de la interacci´on posea el potencial, mejor ser´a la simulaci´ on. La contrapartida de esto es que mientras mayor sea la complejidad funcional del potencial, mayor tambi´en ser´ a el tiempo de computaci´on requerido. Evidentemente, si lo que se busca es s´ olo probar ciertos aspectos de un modelo te´orico, lo mejor ser´a emplear un potencial lo m´as simple posible que reproduzca la esencia de ese modelo. Diferente es la situaci´on si lo que se desea es simular materiales reales: entonces el potencial deber´a contener el m´aximo de informaci´ on posible de modo de reproducir los resultados no s´olo cualitativamente, sino tambi´en cuantitativamente. Sin duda en un s´ olido el mejor m´etodo para obtener las fuerzas que act´ uan sobre los ´atomos es por medio de la Mec´ anica Cu´ antica, resolviendo la ecuaci´on de Schr¨odinger para un sistema de N part´ıculas interactuantes (ve´ ase [22] para una interesante introducci´on). De hecho, ya se han desarrollado m´etodos para realizar esta tarea desde primeros principios, como el M´etodo de Car– Parrinello, que combina DM con teor´ıa del funcional de la densidad [23, 24, 25]. Sin embargo, el costo computacional de esto es alto, pudiendo realizarse en la actualidad simulaciones con a lo m´as cientos de part´ıculas. Si se quiere ir m´as all´a, se debe establecer un compromiso entre la calidad del potencial y las posibilidades de c´alculo. Esto es lo que mantiene viva la vigencia de los llamados potenciales emp´ıricos y semi–emp´ıricos y la b´ usqueda de nuevos m´etodos para mejorarlos (ve´ ase, por ejemplo, el MRS (Material Research Society) Bulletin del mes de febrero de 1996 dedicado al tema [26]). 2.2.1.

La forma del potencial

En general, la energ´ıa potencial V de un sistema de N ´atomos puede expresarse en funci´on de las coordenadas de los ´ atomos individuales, de la distancia entre dos de ellos, de la posici´on relativa entre tres ´ atomos, etc: V=

N X

v1 (ri ) +

i=1

N X N X i=1 j>i

v2 (ri , rj ) +

N X N N X X

v3 (ri , rj , rk ) + · · ·

(3)

i=1 j>i k>j>i

donde el primer t´ermino v1 representa las interacciones de un cuerpo (fuerza externa), v2 las interacciones de dos cuerpos, o de pares, v3 interacciones de tres cuerpos y as´ı sucesivamente. El t´ermino de dos cuerpos, v2 , s´olo depende del m´odulo de la distancia inter´atomica rij = |ri − rj |. Este t´ermino es muy importante pues se ha demostrado que ´el por s´ı solo describe muy bien ciertos sistemas f´ısicos, como es el caso del potencial de Lennard–Jones para los gases nobles. El resto de los t´erminos v3 (ri , rj , rk ), v4 (ri , rj , rk , rl ) . . . son la llamada interacci´on de muchos cuerpos. Estos t´erminos toman en cuenta los efectos de cluster sobre un ´atomo causado por tener m´ as de un ´ atomo alrededor de ´el. Por ejemplo, el t´ermino de tres cuerpos v3 (ri , rj , rk ) es de mucha importancia en el caso de s´olidos covalentes, por los enlaces direccionales que ellos poseen. En el caso de metales, el potencial se puede separar en un t´ermino de dos cuerpos y uno de muchos cuerpos, en la forma de un funcional que depende de la densidad electr´onica local alrededor del ´ atomo en cuesti´ on.

2.3.

Implementaci´ on

Una simulaci´ on t´ıpica, tanto de Din´amica Molecular como por medio de Monte Carlo, implica la elaboraci´ on de un programa de computaci´on cuyos elementos centrales se indican en la Figura 1. Ellos son:

4

Elecci´ on del modelo

Elección del ensemble Número de partículas Condiciones iniciales: (x, v), T, P

Inicialización

DM

Generación de Config.

Elección del potencial Cond.de borde Lista vecinos, etc

MC

DM

Cálculo de fuerzas Integrac. Ec. Newton

Calculo energía V(r) Mét. Estocático para mover partículas (Metropolis)

xi(t), vi(t)

Resultados MC

Análisis de resultados

xi(t)

Cálculo de prop. físicas: Prop. temodinámicas Prop. estáticas Prop. dinámicas

Figura 1: Elementos centrales de un programa t´ıpico de simulaci´on computacional.

Condiciones iniciales y condiciones de borde

5

1.

Inicializaci´ on: una vez realizada la elecci´on del modelo, es decir, el ensemble del sistema y el potencial de interacci´ on, entre otros, se debe hacer una descripci´on molecular del sistema en cuesti´ on, especificando las condiciones iniciales de la variables involucradas, tales como posiciones de las part´ıculas, temperatura, volumen, densidad, etc.

2.

Generaci´ on de las configuraciones: en DM se obtienen integrando num´ericamente las ecuaciones cl´ asicas de movimiento que gobiernan el sistema. En MC las mismas se obtienen mediante un m´etodo de car´ acter estoc´astico.

3.

An´ alisis de los resultados. Se trata de evaluar las propiedades f´ısicas con la informaci´on adquirida. Esto se realiza tomando promedios temporales sobre las diferentes configuraciones. Los valores medios as´ı obtenidos, considerados durante un tiempo suficientemente largo, corresponden a los promedios termodin´amicos suponiendo un comportamiento erg´ odico del sistema.

Como se observa, la diferencia fundamental entre DM y MC est´a en el m´odulo 2, es decir, en la forma como generar las trayectorias del sistema: mientras en DM se hace por medio de la 2da Ley de Newton, lo cual permite calcular propiedades din´amicas, en MC se sigue un m´etodo probabilistico y por tanto no es posible un seguimiento a trav´es del tiempo del sistema, excluyendo as´ı la posibilidad de calcular propiedades din´amicas2 . En lo que sigue, analizaremos cada uno de estos m´odulos de la manera m´as general posible, especificando en cada caso lo que es v´alido para DM y que es v´alido para MC.

3. 3.1.

Inicializaci´ on Condiciones iniciales y condiciones de borde

La especificaci´ on de las condiciones iniciales para la posici´on y la velocidad de cada part´ıcula puede realizarse en una variedad de formas, dependiendo de las caracter´ısticas del sistema a simular. En el caso de s´ olidos perfectos es costumbre poner las part´ıculas inicialmente en sus posiciones de equilibrio en la red, tomando como caja de simulaci´on un m´ ultiplo de la celda unitaria en cada una de las direcciones x, y, z. Si se trata de un s´olido con defectos en la red cristalina, por ejemplo un borde de grano, entonces las part´ıculas se ponen inicialmente en posiciones que se suponen cercanas al equilibrio y se ocupa alguna t´ecnica de minimizaci´on [21], por ejemplo el m´etodo del descenso m´ as r´apido, simulated annealing, quenching molecular dynamics u otro, para minimizar la energ´ıa del sistema, que al inicio de la simulaci´on (a temperatura T = 0 K) corresponde a la energ´ıa de cohesi´on. Los velocidades iniciales se especifican generalmente asign´andoles a cada part´ıcula una velocidad escogida al azar de una distribuci´on de Maxwell–Boltzmann. Estas velocidades iniciales pueden ser escaladas para obtener la temperatura deseada. El momentun lineal y angular del sistema se iguala a cero. Para el caso de sistemas l´ıquidos y amorfos (vidrios) se puede proceder de forma similar. Calentando gradualmente, i.e., escalando las velocidades de los ´atomos, y variando las dimensiones de la caja de simulaci´ on obtenemos la temperatura y densidad deseada. La correcta elecci´ on de la condiciones de borde es otro aspecto que debe considerarse en la simulaci´ on. Para simular el bulk es costumbre emplear condiciones de borde peri´odicas (CBP). En el caso de superficies libres u otros defectos es necesario considerar otro tipo de condiciones de borde. 2 Notemos que en ciertos casos particulares es posible encontrar una prescripci´ on que permite asociar al paso de MC un tiempo real.

6

Generaci´ on de las configuraciones

Figura 2: Supercelda de simulaci´on.

4. 4.1.

Generaci´ on de las configuraciones Evaluaci´ on de la energ´ıa potencial y las fuerzas

La evaluaci´on de la energ´ıa potencial (en MC) y las fuerzas (en DM) es la parte m´as costosa, en t´erminos de tiempo de computaci´ on, en una simulaci´on. De hecho, el tiempo que se ocupa en integrar las ecuaciones es casi despreciable al lado de ´este. Para un sistema de N part´ıculas, evaluar en forma directa la interacci´ on de dos cuerpos requiere O(N 2 ) de operaciones, mientras que evaluar la parte de tres cuerpos requiere en principio O(N 3 ) de operaciones. De all´ı la necesidad de elaborar t´ecnicas que permitan ahorrar tiempo en esta tarea. Truncamiento del Potencial Supongamos que tenemos un sistema de N part´ıculas interactuando a trav´es de potencial de pares con CBP y necesitamos evaluar la energ´ıa potencial y la fuerza sobre una cierta part´ıcula i. La situaci´on la podemos representar por una “supercelda”que consiste en la caja de simulaci´on rodeada por sus im´agenes, como se muestra en la Figura 2. La fuerza sobre la part´ıcula i corresponde a la suma sobre sus N − 1 vecinos de la caja. Pero, debido a las CBP, deber´ıamos tambi´en sumar sobre sus im´agenes. Una manera de hacer esto es mediante la convenci´ on de m´ınima imagen: a partir de la part´ıcula i se construye una caja imaginaria de igual dimensi´on y forma que la caja de simulaci´on, y se suma s´olo sobre las part´ıculas dentro de ella. Claro est´a que para usar la convenci´on de m´ınima imagen se requiere que el alcance del potencial sea menor que la mitad de la longitud de la caja. Esto

7

u ´ltimo permite dar una idea de las dimensiones m´ınimas de la caja de simulaci´on (y por tanto del n´ umero de part´ıculas) a emplear en relaci´on al potencial usado. La contribuci´ on principal a la energ´ıa y a las fuerzas sobre una part´ıcula provienen de sus vecinos m´ as cercanos. Por ello es costumbre ocupar potenciales de corto alcance, que generalmente no van m´ as all´ a de quintos vecinos, e introducir un corte (cutoff), rc , m´as all´a del cual el potencial es nulo. Para asegurar que las fuerzas y la energ´ıa potencial tiendan suavemente a cero en r = rc se puede usar un ajuste del tipo dV V(r) → V(r) − V(r)|r=rc − (r − rc ) . (4) dr r=rc De este modo, la suma ya no se realiza sobre los N − 1 vecinos, sino que queda restringida a los que est´ an dentro de la ‘esfera de influencia del potencial’, como se muestra en la Figura 2. Lista de vecinos Hasta ahora el tiempo de computaci´on que hemos ahorrado consiste en que la evaluaci´ on de la energ´ıa y la fuerza no se hace para N − 1 part´ıculas, sino para un n´ umero mucho menor. Sin embargo, para saber cuales part´ıculas son las que est´an a distancia mayor del corte, y por tanto no contribuyen ni a la energ´ıa ni a la fuerza, debemos examinar, en cada paso de computaci´ on, la distancia entre todos los pares de part´ıculas. El tiempo de esta operaci´on es proporcional a N 2 . Para reducir este tiempo de computaci´on, Verlet [27] ide´o un ingenioso sistema de lista de vecinos de cada part´ıcula , que se renueva cada cierto n´ umero de pasos. El m´etodo supone que los vecinos con los cuales interact´ ua la part´ıcula i, o sea aquellos que est´an dentro de la esfera de radio rc (ver Figura 2) no var´ıan mucho entre paso y paso de integraci´on. A lo m´as cada cierto n´ umero de pasos, digamos 30 por ejemplo, algunas part´ıculas entran y otras salen, quedando a distancias menores que rs . Lo que propuso Verlet fue hacer una lista, para cada part´ıcula, de todos los vecinos que est´ an dentro de su esfera de radio rs , y as´ı en vez de examinar la distancia de la part´ıcula i con todas las N − 1 restantes, se examinan esas distancias s´olo con las part´ıculas de su lista. Esta lista se contruye cada cierto n´ umero de pasos. El ahorro con este m´etodo es significativo para sistemas de entre 500 a 5000 part´ıculas, para los cuales el tiempo por paso de simulaci´on baja pr´acticamente a la mitad (ver [9], p´ag.149). Para sistema de entre 100 a 200 part´ıculas los cambios no son sustanciales, mientras que para sistemas de m´ as de 5000 part´ıculas se han ideado m´etodos m´as eficientes, como el link–cell– list [9].

4.2.

Din´ amica Molecular

Una parte central de todo programa de DM lo constituye el algoritmo de integraci´on. Las ecuaciones de movimiento de Newton dadas por la Ec.( 6) son ecuaciones diferenciales ordinarias acopladas, no–lineales, de segundo orden. Ellas deben ser resueltas num´ericamente. Dadas las posiciones y velocidades iniciales a un tiempo inicial t0 , la tarea del algoritmo es entregar las posiciones y velocidades al tiempo t0 + ∆t. En lo que sigue revisaremos el formalismo b´asico de la DM en el ensemble microcan´onico. En este caso las variables termodin´ amicas que se mantienen constantes son el n´ umero de part´ıculas N , elPvolumen V y la energ´ıa interna E. La Lagrangeana viene dada por L = K − V, donde N K = i=1 (mi /2)q˙ 2i es la energ´ıa cin´etica y V es la energ´ıa potencial. La din´ amica de este sistema est´ a gobernada por las ecuaciones de Euler–Lagrange   ∂L d ∂L − =0, i = 1, . . . , N (5) ∂qi dt ∂ q˙i

8

Generaci´ on de las configuraciones

donde qi y q˙i son las posiciones y las velocidades respectivamente. Esto da origen a las ecuaciones de movimiento mi¨ ri = −∇ri V , i = 1, . . . , N . (6) Cuando las fuerzas entre las part´ıculas son conservativas la Hamiltoniana H es una constante de movimiento y la energ´ıa total se conserva: H=K+V =E .

(7)

Cabe hacer notar que generalmente en Din´amica Molecular el momentum tambi´en es una cantidad conservada cuando las paredes del recipiente son reemplazadas por condiciones de borde peri´odicas. Esto significa que existe una ligadura adicional, que se reflejar´a en las propiedades termodin´ amicas calculadas en la simulaci´on. Sin embargo, se ha demostrado que tales efectos son despreciables para sistemas de m´ as de cien part´ıculas [19, 20]. Exiten numerosos algoritmos para integrar las ecuaciones de Newton. Todos ellos convierten las ecuaciones diferenciales en ecuaciones de diferencias finitas. En DM la elecci´on del algoritmo es (nuevamente) un compromiso entre el grado de precisi´on requerido y el costo computacional. Los algoritmos m´ as usados son el de Verlet [27], el velocity Verlet y el algoritmo de Beeman. Otros algoritmos muy usados en DM son los del tipo corrector–predictor como el de Gear [9]. Revisemos por ejemplo el algoritmo de Verlet. Para deducirlo, partimos del desarrollo en serie de Taylor de r(t), r(t + ∆t) = r(t) + v(t)∆t +

a(t) (∆t)2 + · · · 2

(8)

donde la aceleraci´ on es a(t) = F(t)/m. Del mismo modo, r(t − ∆t) = r(t) − v(t)∆t +

a(t) (∆t)2 + · · · 2

(9)

Sumando ambos desarrollos, obtenemos r(t + ∆t) + r(t − ∆t) = 2r(t) + a(t)(∆t)2 + O((∆t)4 )

(10)

As´ı, la nueva posici´ on r buscada, en el tiempo t + ∆t, viene dada por r(t + ∆t) ≈ 2r(t) − r(t − ∆t)a(t)(∆t)2

(11)

El error estimado que contiene la nueva posici´on r es del oreden de ∆t4 , donde ∆t es el paso de tiempo (‘time step’) en la simulaci´ on de DM. N´otese que para evaluar la nueva posici´on r(t + ∆t) s´ olo necesitamos conocer la posici´on anterior (en t − ∆t) y la aceleraci´on en el tiempo t; no se necesita la velocidad. Sin embargo, ´esta la podemos calcular a partir de r(t + ∆t) − r(t − ∆t) = 2v(t)∆t + O((∆t)3 ) de donde v(t) =

r(t + ∆t) − r(t − ∆t) + O((∆t)2 ) 2∆t

(12)

(13)

Como se ve, el error en la velocidad es del orden ∆t2 , adem´as que no se trata en el mismo nivel que la posici´ on. Un algoritmo que supera este hecho es el velocity Verlet, donde la posici´on y la velocidad se obtiene al mismo tiempo (t + ∆t).

9

Notemos que la forma mostrada de conducir la din´amica, esto es, integrando las ecuaciones de movimiento y avanzado con un tiempo subdividido en porciones δt es eficiente para el caso de potenciales continuos, pero si tenemos un potencial tipo esferas duras de radio ro , o sea  ∞ si r < ro V (r) = 0 si r > ro es mejor usar lo que se denomina din´ amica dirigida por eventos. Supongamos que tenemos un conjunto de n esferas duras en un caja c´ ubica y damos a cada esfera una cierta velocidad inicial (‘temperatura’). Entonces el pr´ oximo evento que ocurrir´a ser´a i) el choque de una esfera con la pared, o ii) el choque entre dos esferas. Dado que sabemos el radio de las esferas, as´ı como la velocidad y posici´ on inicial, se puede calcular exactamente cual ser´a el pr´oximo evento y el tiempo te en el cual ocurrir´ a. Entonces en vez de integrar las ecuaciones de movimiento, es m´as eficiente mover directamente cada una de las esferas desde su posici´on ri a la nueva posici´on r’i = vi te . Tanto para el choque de una esfera contra la pared o el choque entre dos esferas entre ellas, se puede suponer una colisi´on el´astica y por tanto solo cambiar el sentido y direcci´on de la velocidad seg´ un la ley de choques el´asticos. En resumen, un programa de DM trabaja de la siguiente forma: 1.

Se leen los par´ ametros que especifican las condiciones de la corrida tales como la temperatura inicial, el n´ umero de part´ıculas, la posici´on de las part´ıculas (estructura fcc, bcc, etc) la densidad, el paso de tiempo ∆t, tiempo total de simulaci´on, etc.

2.

Se inicializa el sistema, esto es, se asignan las posiciones y las velocidades iniciales.

3.

Se calculan las fuerzas sobre todas las part´ıculas.

4.

Se integran las ecuaciones de movimiento de Newton. Este paso as´ı como el anterior conforman el loop central de la simulaci´on. Ellos son repetidos hasta haber calculado la evoluci´ on temporal del sistema durante el tiempo total de simulaci´on deseado. Se van guardando las posiciones, velocidades, fuerzas, etc, durante cada paso en un archivo para luego ser procesadas.

5.

Despu´es de haber completado lo anterior, se calculan y se imprimen los diferentes promedios relevantes. Entonces STOP.

4.3.

Monte Carlo

El otro m´etodo usado para generar configuraciones independientes en el espacio de fases es el llamado “M´etodo de Montecarlo”. Este es un m´etodo de car´acter probabil´ıstico o estoc´astico, que hace uso intensivo de un generador a n´ umeros aleatorios en su funcionamiento. Precisamente su nombre es en honor a la ciudad de Montecarlo, famosa por sus casinos, ruletas y juegos de azar, entre otras diversiones. Supongamos que deseamos calcular una cierta propiedad A de un sistema en el ensemble can´onico (N V T ). Seg´ un la mec´ anica estad´ıstica, el valor promedio de Atotal esta dado por R hAitotal =

dpN drN exp [−βH(pN , rN )]A(pN , rN ) R dpN drN exp [−βH(pN , rN )]

(14)

10

Generaci´ on de las configuraciones

En general, como la parte de energ´ıa cin´etica K del Hamiltoniano depende cuadr´aticamente de los momenta p, estos grados de libertad se pueden integrar anal´ıticamente, quedando s´olo una integral que depende de las coordenadas r Z hAi = drN A(rN )Z(rN ) (15) con Z(rN ) = R

exp[−βU(r)] drN exp −βU(r)

(16)

Nuestro problema entonces se reduce a que dada la energ´ıa potencial U(rN ), debemos calcular la funci´on de probabilidad Z(rN ) y luego hacer la integral. Para ello podr´ıamos utilizar el m´etodo tradicional de cuadraturas. Sin embargo, a poco andar uno se da cuenta que se trata de una integral multidimensional (en 6N dimensiones, donde N es el n´ umero de part´ıculas) donde los m´etodos de cuadraturas son prohibitivos por la lentitud de ellos. Una alternativa a esto es hacer la integral por el M´etodo de Montecarlo. Supongamos que por alg´ un procedimiento podemos generar aleatoriamente puntos del espacio de configuraciones de acuerdo a la distribuci´on de probabilidades Z(rN ). Esto significa que en promedio, el n´ umero de puntos ni generados por unidad de volumen alrededor del punto riN N es igual a LZ(ri ), donde L es el n´ umero total de puntos generados. O sea L

hAi =

1X ni A(riN ) L i=1

(17)

Claramente generar aleatoriamente los puntos en el espacio de configuraciones no es el mejor m´etodo, pues como Z(rN ) es proporcional al factor de Boltzmann exp(−βU ), los puntos que solo tienen baja energ´ıa contribu´ıran significativamente, mientras que los de alta energ´ıa tendr´an un bajo peso relativo. La clave entonces esta en idear un m´etodo que s´olo genere puntos que tengan un peso relativo alto3 . Esto fue precisamente lo que resolvi´o Metropolis en los a˜ nos cincuenta, con el algoritmo que lleva su nombre. A continuaci´ on damos exposici´ on heur´ıstica de este m´etodo. Consideremos la transici´on entre estados a y b, donde la probabilidad que el sistema est´e en el estado a es pa y en el estado b es pb . En equilibrio se cumple la condici´on del balance detallado, y por tanto Wab pa = Wba pb

(18)

donde Wab es la probabilidad de que ocurra la transici´on a → b y Wba es la probabilidad de que ocurra la transici´ on de b → a. Por otro lado, sabemos que en equilibrio la probabilidad de encontrar al sistema en el estado a(b) es proporcional a exp(−βEa(b) ), luego pa = exp(−β∆Eab ). pb

(19)

Wba = exp(−β∆Eab ) Wab

(20)

Reemplazando en ( 18) queda

As´ı, vemos que tambi´en la probabilidad de transici´on entre dos estados posibles est´a controlada por el factor de Boltzmann. Esto nos da la clave de como escoger la configuraci´on (o sea, como 3 Esto

es lo que se llama importance sampling.

11

hacer el sampling) para tener siempre un peso relativo alto: si al ir del estado a al estado b el cambio de energ´ıa ∆Eab es menor que cero, entonces la probabilidad Wab es 1 (decimos que “aceptamos la movida”). Pero si el cambio ∆Eab > 0, entonces Wab ∝ exp −β∆Eab . En este caso comparamos Wab con un n´ umero r entre [0,1] escogido al azar, y si Wab > r, se acepta la movida, y si es menor se rechaza. En la pr´ actica, si tenemos un sistema de N part´ıculas sometido a una temperatura T , el algoritmo de Metropolis se implementa as´ı: 1.

Seleccione una part´ıcula ni cualquiera y calcule su energ´ıa Ui (r)

2.

D´e a la part´ıcula un desplazamiento aleatorio r0 = r + ∆r y calcule su nueva energ´ıa Ui0 (r).

3.

Si ∆U < 0: acepte la movida y vuelva a 1. Si ∆U > 0, escoja un n´ umero α al azar entre [0,1]. Si α < exp(−β∆U ), acepte y vuelva a 1. Si α > exp(−β∆U ), rechaze la movida, es decir conserva la posici´on de la part´ıcula ni , y vaya a 1.

4.

Luego de barrer la N part´ıculas, guarde la configuraci´on obetnida y vuelva a 1 para comenzar otro ciclo. Despues de un n´ umero razonable de paso, vaya a 5

5.

Calcule la propiedades f´ısicas de inter´es a partir de la configuraciones guardadas. STOP.

Antes de terminar este apartado hacemos las siguientes obeservaciones: a menudo es costumbre usar, en vez de la probabilidad exp(−β∆U ), la probabilidad normalizada 1/(1 + exp(β∆U )). Esto no cambia los resultados del equilibrio, y tiene la ventaja que el m´etodo converje m´as r´apido para altas temperaturas. En t´erminos t´ecnicos, podemos decir que el algoritmo de Metropolis es un proceso de Markov en el cual se construye una caminata aleatoria (random walks) de tal modo que la probabilidad de visitar un punto particular ~rN es proporcional al factor de Boltzmann exp(−βU ). De hecho hay varios maneras de construir tal caminata aleatoria. El algoritmo de Metropolis es s´ olo una de ellas, la m´as usada. Una de la mayores diferencias entre MC y DM es que en MC no es posible calcular propiedades din´ amicas, pues la din´amica que se sigue es ficticia. ¿Entonces cual es la ventaja de MC? La respuesta depende del tipo de problema a tratar. Seg´ un Frenkel [13], siempre que se pueda, deber´ıa preferirse DM. Sin embargo, hay algunos casos (que no son pocos) en los cuales no es posible usar DM: i) sistemas que no poseen din´amica intr´ınseca, como modelos en la red, sistemas de espines (Ising, Heisenberg, modelo de Potts, etc). ii) casos en que las barreras de activaci´on son demasiado altas.

5.

An´ alisis de los resultados

Lo que hemos tratado hasta ahora dice relaci´on con la simulaci´on propiamente tal, esto es, como obtener la trayectoria en el espacio de fases para un sistema de part´ıculas. Ahora veremos como analizar estas trayectorias para obtener de all´ı propiedades f´ısicas macrosc´opicas que puedan ser comparadas con el experimento. Dividiremos estas propiedades en: propiedades termodin´ amicas, propiedades est´ aticas y propiedades din´amicas.

12

5.1.

Evaluaci´ on de las propiedades f´ısicas

Propiedades termodin´ amicas

En el ensemble microcan´ onico la temperatura del sistema se calcula como el promedio de la energ´ıa cin´etica, a trav´es del teorema de equipartici´on hKi =

3 N kB T . 2

(21)

donde kB = 1,38062 × 10−23 J/K es la constante de Boltzmann. El h. . .i se refiere al promedio sobre las N part´ıculas y sobre el ensemble. La presi´on media, p, se calcula a trav´es del teorema del virial pV = N kB T + hW i , (22) donde 1 hW i = 3

*N X

+ ri · Fi

(23)

i=1

es conocido como el virial del sistema, siendo ri la posici´on de la part´ıcula i y Fi la fuerza sobre la part´ıcula i debido a todas las restantes. A partir de las fluctuaciones cuadr´ aticas medias, h(δA)2 i = h(A − hAi)2 i, se obtienen las funciones de respuesta, tales como calor espec´ıfico, m´odulo de compresibilidad adiab´atico y constantes el´ asticas. Por ejemplo, el calor espec´ıfico por part´ıcula, CV , se obtiene de   h(δK)2 i h(δV)2 i 3kB = = N kB T 1 − . (24) hKi hVi 2CV Este tipo de relaciones entre una funci´ on respuesta y la desviaci´on cuadr´atica media de cierta cantidad es un caso especial del teorema de fluctuaci´on–disipaci´on. Su uso en DM ofrece la clara ventaja que mediante una sola corrida de simulaci´on podemos obtener inmediatamente la funci´on de respuesta deseada. Un m´etodo alternativo para obtener, por ejemplo, el calor espec´ıfico, ser´ıa correr el sistema a varias temperaturas diferentes y graficar la energ´ıa total respecto de la temperatura. Su pendiente nos dar´ıa CV , con el evidente gasto de tiempo. Otras propiedades importantes a calcular son la entrop´ıa y las energ´ıas libres. Desgraciadamente esto no es tan directo como lo mostrado arriba y requiere de t´ecnicas m´as elaboradas [28].

5.2.

Propiedades est´ aticas

Las propiedades estructurales est´ aticas de un sistema se pueden describir a trav´es de la funci´on de distribuci´ on de pares (FDP), g(r) y del factor de estructura est´atico S(q). La primera est´ a dada por hn(r, r + ∆r)i V g(r) = , (25) 4πr2 ∆r N donde n(r, r +∆r) indica el n´ umero de part´ıculas que hay en una capa entre r y r +∆r, teniendo como origen una determinada part´ıcula. La funci´on de distribuci´on de pares es proporcional a la probabilidad de encontrar dos part´ıculas separadas por una distancia r + ∆r. Es com´ un tambi´en graficar la funci´ on de distribuci´on radial F DR = 4πρ0 g(r)r2 ; aqu´ı el ´area encerrada por el primer pico es proporcional al n´ umero de coordinaci´on mientras que el cuociente entre la posici´on del primer y segundo pico informa sobre las distancias interat´omicas. Experimentalmente lo que se mide, mediante scattering de neutrones y rayos X, es el factor de estructura est´ atico, S(q). En el caso de l´ıquidos y materiales amorfos, este s´olo depende del

Propiedades est´aticas

13

Figura 3: Funci´ on de distribuci´on de pares y patr´on de difracci´on

m´odulo |q| y puede representarse como una integral sobre g(r): Z S(q) = 1 + 4πρ 0

R

r2 [g(r) − 1]

sen(qr) dr , qr

(26)

donde el valor de R debe escogerse menor que la mitad de la longitud de la caja de simulaci´on. En la Figura 3 se muestra esquem´aticamente la forma t´ıpica que presenta la funci´on de distribuci´ on de pares y el patr´ on de difracci´on para gases, l´ıquidos, amorfos y solidos, en referencia a la distribuci´ on de sus ´ atomos en el espacio real. En el caso de s´ olidos cristalinos una cantidad importante, que permite caracterizar el desorden de las diferentes capas at´ omicas, es el factor de estructura est´atico definido como * S(K, `) =

2 + X N` 1 exp(iK · r ) , j 2 N` j=1

(27)

donde K es un vector de la red rec´ıproca de una capa at´omica, N` es el n´ umero total de part´ıculas en la capa ` y el ´ındice j se refiere a cada una de las part´ıculas en esa capa.

14

Evaluaci´ on de las propiedades f´ısicas

5.3.

Propiedades din´ amicas

Una de la ventajas de DM con respecto a MC reside en la posiblidad que ofrece de calcular propiedades din´ amicas. Como hemos visto, en MC no hay una din´amica genuina y por tanto no se pueden obtener directamente de ´el propiedades que dependan de la evoluci´on temporal del sistema. DM permite calcular ´estas en forma directa. Las propiedades din´ amicas las obtenemos a partir de las funciones de correlaci´on temporal. Estas son muy importantes en simulaci´on computacional pues est´an directamente relacionadas con los coeficientes de transporte y por tanto, con los datos experimentales. El coeficiente de difusi´ on D, por ejemplo, lo podemos calcular a partir de la funci´on de autocorrelaci´ on de velocidades 1 Z(t) = hvi (t) · vi (0)i 3 a trav´es de Z ∞ D=

Z(t)dt .

(28)

0

Este es un ejemplo de la f´ ormula de Green–Kubo, que permite expresar un coeficiente de transporte (macrosc´ opico) como una integral sobre el tiempo de una funci´on de autocorrelaci´on temporal (microsc´ opica). Una manera alternativa de calcularlo es usando la relaci´on de Einstein h|ri (t) − ri (0)|2 i , t→∞ 6t

D = l´ım

(29)

donde h|ri (t)−ri (0)|2 i es la desviaci´ on cuadr´atica media de la posici´on de la part´ıcula i respecto a su posici´ on inicial ri (0). Otra cantidad que es posible calcular a partir de Z(t) es la densidad de estados vibracionales, G(ω), ˜ Z(ω) G(ω) = , Z(t = 0) ˜ donde Z(ω) es la transformada de Fourier de Z(t) y est´a dado por Z ∞ 1 ˜ eiωt Z(t)dt . Z(ω) = 2π −∞

A. A.1.

(30)

Ap´ endice Comentario Bibliogr´ afico

Sobre f´ısica computacional y simulacion en general: Kauffamn y Smarr: libro de divulgaci´on sobre la relevancia que ha tomado y tomar´a la simulaci´ on computacional en pr´ acticamente todas las ´areas de la ciencia. Computing in Science and Engineering, revista bimensual, heredera de Computer in Physics, permite estar al d´ıa del estado del arte en este tema. Gould y Tobochnik: abarca m´etodos num´ericos y t´ecnicas de simulaci´on usadas en f´ısica. Muy pedag´ ogico y con gran cantidad de u ´tiles ejemplos. Altamente recomendable.

Propiedades din´amicas

15

Koonin y Meredith: libro muy utilizado, que pone ´enfasis en lo m´etodos num´ericos m´as usados hoy d´ıa en f´sica, con su respectiva implementaci´on computacional. Sobre MD: Allen y Tildesley: excelente libro para iniciarse en los secretos de la DM y MC, con discusiones precisas y ejemplos de su implementacion computacional, incluidos muchos programas de computaci´ on. Haile: m´ as que recetas, hay un tratamiento detallado de los fundamentos de la din´amica molecular, con muchos ejemplos y problemas. Buen complemento de Allen y Tildesley. Frenkel y Smit: en la misma l´ınea que el de Allen y Tildesley, es un complemento obligado por traer temas m´ as actuales y dar explicaciones muy claras de los distintos m´etodos usados. La discusi´ on sobre equilibrios de fase es notable. Adem´as esta muy bien escrito y llevado lo que hace que la lectura sea muy entretenida. Excelente. Recular y Madden: un paper ”cl´asico”donde se explica de manera legible el metodo de de Car-Parrinello sobre DM ab-initio. Se requiere conocimientos b´asicos de Mec´anica Cu´ antica y DM cl´ asica. Altamente recomendable. Sobre MC Allen y Tildesley: excelente libro para iniciarse en los secretos de la DM y MC, con discusiones precisas y ejemplos de su implementacion computacional, incluidos muchos programas de computaci´ on. Frenkel y Smit: en la misma l´ınea que el de Allen y Tildesley, es un complemento obligado por traer temas m´ as actuales y dar explicaciones muy claras de los distintos m´etodos usados. Excelente. Binder: referencia obligada para cualquier interesado en profundizar sobre el tema.

A.2.

Direcciones en Internet http://www.sissa.it/furio/md : p´agina de Furio Ercolessi, clave en MD. All´ı hay unos excelentes apuntes de MD con ejemplos en Fortran 90 y links muy interesantes. En la SISSA y en el ICTP, en Trieste, Italia, fue donde naci´o el m´etodo de MD ab–initio de Car–Parrinello. La p´ agina contiene los distintos grupos que alli hay: estructura electr´onica, superficies, as´ı como papers y tesis de doctorado donde se ha usado DM. En Inglaterra hay un grupo llamado Computer Simulation of Condensed Phase: http://www.dl.ac.uk/CCP/CCCP5/main.html aqui hay notas sobre MC y MD, software, noticias, links, etc. http://antas.agraria.uniss.it/software : aqu´ı hay un listado de software disponible, tanto libre como comercial, en DM, MC, DM ab-initio, etc. En DM son recomendables los programas MOLDY y DL POLY. http://www.ncsa.uiuc.edu/Apps/CMP/ceperley : David Ceperley es uno de los expertos mundiales en MC y su aplicaci´on en Materia Condensada. En este sitio hay apuntes, ejemplo de software, e interesantes links.

16

A.3.

Evaluaci´ on de las propiedades f´ısicas

Simulaci´ on computacional en Chile

A continuaci´ on damos una lista (ciertamente incompleta e imprecisa) de personas que trabajan en simulaci´ on computacional (o la usan como herramienta auxiliar) en nuestro pa´ıs. Universidad Cat´ olica de Antofagasta: Sergio Curilef utiliza DM para estudiar la estad´ıstica de Tsallis. Universidad de Chile: • Facultad de Ciencias F´ısicas y Matem´aticas: existe un buen grupo formado por Patricio Cordero, Rodrigo Soto, Rosa Jim´enez y J.M Pasini. Han usado DM dirigida por los eventos para estudiar gases y l´ıquidos. Ver http://www.tamarugo.uchile.cl/ • Facultad de Ciencias: Jos´e Rogan y Rodrigo Ferrer est´an usando MC para estudiar sistemas magn´eticos. Jaime Roessler tiene experiencia en MC. Universidad de Santiago: Carlos Esparza ha trabajado desarrollando nuevos m´etodos de DM. Gonzalo Guti´errez usa DM para caracterizar materiales l´ıquidos y amorfos, as´ı como superficies e interfaces de ´solidos. Dora Altbir usa MC para estudiar sistemas magn´eticos. Guillermo Palma y Teodoro Meyer tienen experiencia con MC en la red. P. Universidad Cat´ olica: Miguel Kiwi y Ricardo Ram´ırez usan DM para estudiar s´olidos. Aldo Romero tiene experiencia en DM ab-initio, Car–Parrinello. Universidad del Bio–Bio: Dino Risso investiga en hidrodin´amica usando DM de esferas duras. Universidad de La Frontera: el grupo de f´ısica del s´olido (Eugenio Vogel et al) tiene vasta experiencia en MC para sistemas de espines. Universidad de Magallanes: Mauricio Marin ha desarrollado algoritmos eficientes en DM.

Referencias [1] William J. Kaufmann y Larry L. Smarr, Supercomputing and the Transformation of Science, Scientific American Library, New York, 1993. [2] New dimensions in simulation, Special issue, en Physics World 9, No. 7, p.29–48 (1996). [3] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J.Chem. Phys. 21, 1087 (1953). [4] H. Gould and J. Tobochnik, An Introduction to Computer Simulation Methods: Applications to Physical Systems, second edition, , Addison–Wesley, Reading MA (1996). [5] S. E. Koonin and D. C. Meredith, Computational Physics, Addison–Wesley, Reading MA (1990). [6] K. Binder, Montercarlo Method in Statistical Physics, Springer, Berlin, 1986. [7] D. Ceperley, Apuntes en p´ agina WEB, http://www.ncsa.uiuc.edu/Apps/CMP/ceperley

Referencias

17

[8] D. W. Heerman, Computer Simulation Methods in Theoretical Physics, Springer-Verlag, 1986. [9] M. P. Allen and D. Tildesley, Computer Simulations of Liquids, Clarendon Press, Oxford, 1987. [10] J. P. Hansen, An introduction to Molecular Dynamics with an applicaction to glass transition, en Computer Simulation in Materials Science, Edited by M. Meyer and V. Pontikis, Kluwer Academic Publishers, 1991. [11] J. H. Haile, Molecular Dynamics Simulation, J. Wiley, New York, 1992. [12] D. E. Rappaport, The Art of Molecular Dynamics Simulation, Cambridge Univ. Press, 1996. [13] D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, San Diego, 1996. [14] Simulation of Liquids and Solids, edited by G. Cicotti, D. Frenkel, and I. R. Mc Donald, North-Holland, Amsterdam, 1987. [15] D. Frenkel, Introduction to Computer Simulation, en Simple Molecular Systems at very High Density, Edited by A. Polian, P Loubeyre y N. Boccara, Plenum Press, New York, 1989. [16] S. Nose, Mol. Phys. 50, 255 (1983). [17] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980). [18] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981). [19] F. Lado, J. Chem. Phys. 75, 5461 (1981). [20] D. C. Wallace and G. K. Straub, Phys. Rev. A 27, 2201 (1983). [21] W. H. Press, S. A. Teulkolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes: the art of scientific computing, Second Edition, Cambridge University Press, 1992. [22] M. Gillan, Contemporary Physics 38, No.2, p.115–130 (1997). [23] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). [24] D. K. Reculer and P. A. Madden, Molecular dynamics without effective potentials via the Car–Parrinello approach, Mol. Phys. 70, 921 (1990). [25] G. Galli and A. Pasquarello M. Parrinello, First Principles Molecular Dynamics, en Computer Simulation in Chemical Physics, Edited by M. P. Allen and D. J. Tildesley, Kluwer Academic Publishers, 1993. [26] MRS Bulletin, Vol. 21, No.2, Feb. 1996. [27] L. Verlet, Phys. Rev. 159, 98 (1967). [28] D. Frenkel, Lectures Notes on Free Energy Calculations, en Computer Simulation in Materials Science, Edited by M. Meyer and V. Pontikis, Kluwer Academic Publishers, 1991.