ELEMENTOSDESIMULACIONCOMPUTACIONALDinamicaMolecularyMetododeMonteCarloGonzalo GutierrezDepartamento de Fsica,Universidad de Santiago de Chile.email: ggutierr@lauca.usach.clAbril 2001Santiago, ChileIndice iIndice1. Introduccion 12. ElecciondelModelo 12.1. Ensembles de la Simulacion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.2. El potencial interatomico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.2.1. La forma del potencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.3. Implementacion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33. Inicializacion 53.1. Condiciones iniciales y condiciones de borde. . . . . . . . . . . . . . . . . . . . . 54. Generaciondelasconguraciones 64.1. Evaluacion de la energa potencial y las fuerzas . . . . . . . . . . . . . . . . . . . 64.2. Dinamica Molecular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74.3. Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95. Analisisdelosresultados 115.1. Propiedades termodinamicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125.2. Propiedades estaticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125.3. Propiedades dinamicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14A. Apendice 14A.1. Comentario Bibliograco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14A.2. Direcciones en Internet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15A.3. Simulacion computacional en Chile . . . . . . . . . . . . . . . . . . . . . . . . . . 16ListadeAbreviacionesMC : Monte CarloDM : Dinamica MolecularCBP: Condiciones de Borde PeriodicasFDP: Funcion de Distribucion de paresAcercadeestedocumentoNotas de clase con motivo del curso Simulacion en Materia Condensada, Depto. Fsica, Usach,2001.iic _GonzaloGutierrez,2da.version,Abril2001.Enviarcrticas,comentarios,sugerenciasaggutierr@lauca.usach.clIntroduccion 11. IntroduccionEl desarrollo de los computadores digitales a partir de la decada de los 50, y su aplicacion ala resolucion de problemas cientcos, ha introducido lo que algunos han llamado una tercerametodologaalainvestigacioncientca: lasimulacioncomputacional [1]. Estemetodo, decaracter complementario y muchas veces alternativo a los modos convencionales de hacer ciencia,el experimental y el teorico, ha ejercido un fuerte impacto en practicamente todos los camposdelaciencia(verporejemplo[1, 2]). El objetivodelasimulacioncomputacional esresolverlos modelos teoricos en su total complejidad, mediante la resolucion numerica de las ecuacionesinvolucradas, haciendo uso intensivo (y extensivo) de computadores.En el area de la fsica, la simulacion computacional fue introducida como una herramientaparatratarsistemasdemuchoscuerposacomienzodeladecadadelos50, conel trabajopionero de Metropolis et al. [3]. Mas tarde, auspiciosos resultados iniciales obtenidos en mecanicaestadsticaclasica,enparticularenelestudiodelquidos,dieroncredibilidadalasimulacioncomputacional,extendiendoserapidamentesuusoatemastandiversoscomocromodinamicacuantica, fsicadeudos, relatividadgeneral, fsicadel plasma, materiacondensada, fsicanuclear y ciencia de materiales.Actualmente, gracias al vertiginoso desarrollo de la tecnologa de los computadores, cuya ve-locidad crece aproximadamente un factor 2 cada dieciocho meses, la simulacion computacionalsehaconstitudoenunaherramientadecalculoesencial, tantoparaexperimentalistascomopara teoricos. Mediante un buen modelo computacional no solo se pueden reproducir experimen-tos de laboratorio, sino que ademas, gracias a que se pueden variar libremente los parametrosusados, permiteprobar (ofalsar) modelos teoricos existentes enrangos deparametros im-posibles de alcanzar experimentalmente por ahora, resolviendo as conictos entre explicacionteoricayobservacion.Unpapelfundamentaltambienlojuegahoydalavisualizaciondelosresultados obtenidos. No solo obtenemos datos numericos que pueden ser contrastados con losexperimentos, sino tambien obtenemos una imagen graca del proceso en cuestion.Los dos metodos de simulacion computacional mas usados en fsica actualmente son el dela Dinamica Molecular [8, 9, 10, 11, 12], que es de caracter determinista, y el de Montecarlo,queesdecaracterprobablistico[6].Ambospuedenconsiderarsecomometodosparagenerarconguracionesdiferentesdeunsistemadepartculas1,esdecirpuntosenelespaciodefasescompatibles con las condiciones externas.El metodo de la Dinamica Molecular y el de Montecarlo ha sido empleado con exito parasimular gases, lquidos y solidos [14], ampliandose tanto su uso como el desarrollo de tecnicasespeccasenformaparalelaalavancetecnologicodeloscomputadores.Lossistemasestudi-adosvandesdecientosamilesy ultimamenteinclusoadecenasdemillonesdeatomos. Losaspectos estudiados incluyen propiedades estructurales, termodinamicas, mecanicas y cineticas.Cabe se nalar que estas tecnicas de simulacion tienen aplicaciones mucho mas amplias que lasaqu esbozadas. Enparticularel metododeMCseempleaconexitoenFsicadePartculas(por ejemplo lattice gauge theory), as como para calcular sistemas cuanticos, donde se ocupantecnicas tales como Path integral MC y otros [7].2. ElecciondelModeloElpuntodepartidaparasimularunsistemafsicoesdenirconclaridadelproblemaencuestion: quetipodepropiedadesnosinteresaestudiar, dentrodequerangodeparametros,1Porpartculasentendemosatomosomoleculas2 Elecciondelmodeloconqueprecision.Enfunciondeellodebemosdecidireln umerodepartculasausar,cualesseran las variables de control, que potencial interatomico usar, que tipo de promedios debemoscalcular, en que ensemble conducir la simulacion.2.1. EnsemblesdelaSimulacionLa informacion que genera una corrida de Dinamica Molecular es la posicion y la velocidadde cada partcula del sistema en cada instante de tiempo. Por su parte, en el caso de MC lo queobtenemos es la posicion de las partculas en cada paso de simulacion. Empleando las tecnicastradicionales de la Mecanica Estadstica podemos pasar de esta informacion microscopica a laobtencion de magnitudes macroscopicas que nos permitan conectar con el experimento, va latermodinamica.Supongamos que estamos tratando un sistema puro compuesto por Npartculas, encerradoenunvolumenV yconunaenergajaE. Lasposicionesyvelocidadesdenenunespaciode fases de 6Ndimensiones. Obtener la posicion y la velocidad de cada una de las partculas,en cada instante, signica obtener la trayectoria de un punto del espacio de fase en funciondel tiempo, estoes, (t). Denotemospor /el valorinstantaneodeunciertoobservable. Elpromedio de esta cantidad / esta dado por/)obs = /)tiempo =1obsobs

=1/(()) , (1)donde representauntiempodiscreto(pasos dedinamicamolecular)yobssonlospasostotales de la corrida. Suponiendo que el sistema es ergodico, podemos asociar directamente estepromedio con el promedio usual sobre ensemble de la Mecanica Estadstica,/)obs = /)tiempo = /)ens. (2)Enotraspalabras, pormediodel formalismodelasimulacionloquesehaceesgenerarunasucesion de diferentes estados (puntos) del espacio de fases (que se suponen nocorrelacionados),compatibles con las condiciones externas ((N, V, E) es este caso), sobre los cuales se toman lospromedios.La eleccion del ensemble bajo el cual llevar a cabo la simulacion esta dictada fundamental-menteporeltipodeproblemaatratar.Lospromediosestadsticospuedenllevarapeque nasdiferencias en los diferentes ensembles, pero estas desaparecen en el lmite termodinamico, quese alcanza incluso con unos pocos cientos de partculas [9, 13]. Sin embargo la eleccion del en-semble si inuye al momento de calcular las uctuaciones cuadraticas medias de las magnitudestermodinamicas. Estaspermitencalcular, porejemplo, lacapacidadcaloricaoel modulodeelasticidad.El ensemble convencional en Dinamica Molecular es el ensemble microcanonico, (N, V, E),mientras que en Montecarlo es el canonico, (N, V, T), donde T es la temperatura. Tambien se handesarrollado tecnicas para llevar a cabo simulaciones de MC en otros ensembles, as como simu-laciones de DM en el ensemble canonico [16], en el ensemble isotermicoisobarico (N, P, T) [17],apresionytemperaturaconstanteyenel ensembleisoentalpicoisotension(H, t, N)[18], atensionexternatconstante, entreotros. Paraunarevision, conlareferenciaalostrabajosoriginales, de los diferentes metodos existentes remitimos al lector a la Ref. [13].2.2. ElpotencialinteratomicoUnpuntodeimportanciacentraltantoenDinamicaMolecularcomoenMontecarloeslaeleccion del potencial interatomico del sistema a simular. De la delidad con que este representeElpotencialinteratomico 3las interacciones reales entre las partculas dependera la calidad de los resultados: la conclusioninmediataes quemientras mas detalles delainteraccionposeael potencial, mejor seralasimulacion.Lacontrapartidadeestoesquemientrasmayorsealacomplejidadfuncionaldelpotencial, mayor tambien sera el tiempo de computacion requerido. Evidentemente, si lo que sebusca es solo probar ciertos aspectos de un modelo teorico, lo mejor sera emplear un potenciallomassimpleposiblequereproduzcalaesenciadeesemodelo.Diferenteeslasituacionsiloque se desea es simular materiales reales: entonces el potencial debera contener el maximo deinformacion posible de modo de reproducir los resultados no solo cualitativamente, sino tambiencuantitativamente.Sin duda en un solido el mejor metodo para obtener las fuerzas que act uan sobre los atomoses por medio de la Mecanica Cuantica, resolviendo la ecuacion de Schrodinger para un sistema deNpartculas interactuantes (vease [22] para una interesante introduccion). De hecho, ya se handesarrollado metodos para realizar esta tarea desde primeros principios, como el Metodo de CarParrinello, que combina DM con teora del funcional de la densidad [23, 24, 25]. Sin embargo,el costo computacional de esto es alto, pudiendo realizarse en la actualidad simulaciones con alo mas cientos de partculas. Si se quiere ir mas alla, se debe establecer un compromiso entrela calidad del potencial y las posibilidades de calculo. Esto es lo que mantiene viva la vigenciade los llamados potenciales empricos y semiempricos y la b usqueda de nuevos metodos paramejorarlos (vease, por ejemplo, el MRS (Material Research Society) Bulletin del mes de febrerode 1996 dedicado al tema [26]).2.2.1. LaformadelpotencialEn general, la energa potencial 1de un sistema deNatomos puede expresarse en funcionde las coordenadas de los atomos individuales, de la distancia entre dos de ellos, de la posicionrelativa entre tres atomos, etc:1 =N

i=1v1(ri) +N

i=1N

j>iv2(ri, rj) +N

i=1N

j>iN

k>j>iv3(ri, rj, rk) + (3)dondeelprimerterminov1representalasinteraccionesdeuncuerpo(fuerzaexterna), v2lasinteracciones de dos cuerpos, o de pares,v3interacciones de tres cuerpos y as sucesivamente.El termino de dos cuerpos,v2, solo depende del modulo de la distancia interatomicarij=[ri rj[.Esteterminoesmuyimportantepuessehademostradoque elporssolodescribemuy bien ciertos sistemas fsicos, como es el caso del potencial de LennardJones para los gasesnobles. El resto de los terminosv3(ri, rj, rk), v4(ri, rj, rk, rl) . . . son la llamada interaccion demuchos cuerpos. Estos terminos toman en cuenta los efectos de cluster sobre un atomo causadopor tener mas de un atomo alrededor de el. Por ejemplo, el termino de tres cuerpos v3(ri, rj, rk)es de mucha importancia en el caso de solidos covalentes, por los enlaces direccionales que ellosposeen.Enelcasodemetales,elpotencialsepuedesepararenunterminodedoscuerposyunodemuchoscuerpos,enlaformadeunfuncionalquedependedeladensidadelectronicalocal alrededor del atomo en cuestion.2.3. ImplementacionUna simulacion tpica, tanto de Dinamica Molecular como por medio de Monte Carlo, implicalaelaboraciondeunprogramadecomputacioncuyos elementos centrales seindicanenlaFigura 1. Ellos son:4 ElecciondelmodeloInicializacinGeneracin deConfig.Clculo de fuerzasIntegrac. Ec. NewtonCalculo energa V(r)Mt. Estoctico paramover partculas(Metropolis)Eleccin del potencialCond.de bordeLista vecinos, etcDMMCEleccin del ensembleNmero de partculasCondiciones iniciales:(x, v),T, PResultadosxi(t), vi(t)DMMC xi(t)Anlisis de resultadosClculo de prop. fsicas:Prop. temodinmicasProp. estticasProp.dinmicasFigura 1: Elementos centrales de un programa tpico de simulacion computacional.Condicionesinicialesycondicionesdeborde 51. Inicializacion: una vez realizada la eleccion del modelo, es decir, el ensemble del sistemayel potencial deinteraccion, entreotros, sedebehacerunadescripcionmoleculardelsistemaencuestion, especicandolascondicionesinicialesdelavariablesinvolucradas,tales como posiciones de las partculas, temperatura, volumen, densidad, etc.2. Generacion de las conguraciones: en DM se obtienen integrando numericamente las ecua-cionesclasicasdemovimiento que gobiernanelsistema.En MC las mismas se obtienenmediante un metodo de caracter estocastico.3. Analisis de los resultados. Se trata de evaluar las propiedades fsicas con la informacionadquirida. Estoserealizatomandopromediostemporalessobrelasdiferentescongu-raciones. Losvaloresmediosas obtenidos, consideradosduranteuntiemposuciente-mente largo, corresponden a los promedios termodinamicos suponiendo un comportamien-to ergodico del sistema.Comoseobserva, ladiferenciafundamental entreDMyMCestaenel modulo2, esdecir,enla forma como generar las trayectorias delsistema:mientras en DM se hace por medio dela2daLeydeNewton, locual permitecalcularpropiedadesdinamicas, enMCsesigueunmetodo probabilistico y por tanto no es posible un seguimiento a traves del tiempo del sistema,excluyendo as la posibilidad de calcular propiedades dinamicas2.En lo que sigue, analizaremos cada uno de estos modulos de la manera mas general posible,especicando en cada caso lo que es valido para DM y que es valido para MC.3. Inicializacion3.1. CondicionesinicialesycondicionesdebordeLa especicacion de las condiciones iniciales para la posicion y la velocidad de cada partculapuederealizarseenunavariedaddeformas, dependiendodelascaractersticasdel sistemaasimular.Enelcasodesolidosperfectosescostumbreponerlaspartculasinicialmenteensusposicionesdeequilibrioenlared, tomandocomocajadesimulacionunm ultiplodelaceldaunitaria en cada una de las direccionesx, y, z. Si se trata de un solido con defectos en la redcristalina, por ejemplo un borde de grano, entonces las partculas se ponen inicialmente en posi-ciones que se suponen cercanas al equilibrio y se ocupa alguna tecnica de minimizacion [21], porejemplo el metodo del descenso mas rapido, simulated annealing, quenching molecular dynam-ics u otro, para minimizar la energa del sistema, que al inicio de la simulacion (a temperaturaT= 0 K) corresponde a la energa de cohesion.Los velocidades iniciales se especican generalmente asignandoles a cada partcula una ve-locidad escogida al azar de una distribucion de MaxwellBoltzmann. Estas velocidades inicialespueden ser escaladas para obtener la temperatura deseada. El momentun lineal y angular delsistema se iguala a cero.Para el caso de sistemas lquidos y amorfos (vidrios) se puede proceder de forma similar. Ca-lentando gradualmente, i.e., escalando las velocidades de los atomos, y variando las dimensionesde la caja de simulacion obtenemos la temperatura y densidad deseada.La correcta eleccion de la condiciones de borde es otro aspecto que debe considerarse en lasimulacion. Para simular el bulk es costumbre emplear condiciones de borde periodicas (CBP).En el caso de supercies libres u otros defectos es necesario considerar otro tipo de condicionesde borde.2NotemosqueenciertoscasosparticularesesposibleencontrarunaprescripcionquepermiteasociaralpasodeMCuntiemporeal.6 GeneraciondelasconguracionesFigura 2: Supercelda de simulacion.4. Generaciondelasconguraciones4.1. EvaluaciondelaenergapotencialylasfuerzasLa evaluacion de la energa potencial (en MC) y las fuerzas (en DM) es la parte mas costosa,en terminos de tiempo de computacion, en una simulacion. De hecho, el tiempo que se ocupaen integrar las ecuaciones es casi despreciable al lado de este. Para un sistema de Npartculas,evaluar en forma directa la interaccion de dos cuerpos requiere O(N2) de operaciones, mientrasqueevaluarlapartedetrescuerposrequiereenprincipioO(N3)deoperaciones. Deall lanecesidad de elaborar tecnicas que permitan ahorrar tiempo en esta tarea.TruncamientodelPotencial Supongamos que tenemos un sistema deNpartculas inter-actuando a traves de potencial de pares con CBP y necesitamos evaluar la energa potencial yla fuerza sobre una cierta partculai. La situacion la podemos representar por una supercel-daqueconsisteenlacajadesimulacionrodeadaporsusimagenes, comosemuestraenlaFigura2.LafuerzasobrelapartculaicorrespondealasumasobresusN 1vecinosdelacaja. Pero, debido a las CBP, deberamos tambien sumar sobre sus imagenes. Una manera dehacer esto es mediante la convencion de mnima imagen: a partir de la partculai se construyeunacajaimaginariadeigual dimensionyformaquelacajadesimulacion, ysesumasolosobre las partculas dentro de ella. Claro esta que para usar la convencion de mnima imagense requiere que el alcance del potencial sea menor que la mitad de la longitud de la caja. Esto7 ultimo permite dar una idea de las dimensiones mnimas de la caja de simulacion (y por tantodel n umero de partculas) a emplear en relacion al potencial usado.La contribucion principal a la energa y a las fuerzas sobre una partcula provienen de susvecinos mas cercanos. Por ello es costumbre ocupar potenciales de corto alcance, que general-mente no van mas alla de quintos vecinos, e introducir un corte (cuto), rc, mas alla del cualel potencial es nulo. Para asegurar que las fuerzas y la energa potencial tiendan suavemente acero enr = rcse puede usar un ajuste del tipo1(r) 1(r) 1(r)[r=rc (r rc)d1drr=rc. (4)De este modo, la suma ya no se realiza sobre losN 1 vecinos, sino que queda restringida alos que estan dentro de la esfera de inuencia del potencial, como se muestra en la Figura 2.Lista de vecinos Hasta ahora el tiempo de computacion que hemos ahorrado consiste en quela evaluacion de la energa y la fuerza no se hace paraN 1 partculas, sino para un n umeromucho menor. Sin embargo, para saber cuales partculas son las que estan a distancia mayor delcorte, y por tanto no contribuyen ni a la energa ni a la fuerza, debemos examinar, en cada pasode computacion, la distancia entre todos los pares de partculas. El tiempo de esta operaciones proporcional aN2.Para reducir este tiempo de computacion, Verlet [27] ideo un ingenioso sistema de lista devecinos de cada partcula , que se renueva cada cierto n umero de pasos. El metodo supone quelos vecinos con los cuales interact ua la partcula i, o sea aquellos que estan dentro de la esfera deradiorc (ver Figura 2) no varan mucho entre paso y paso de integracion. A lo mas cada cierton umerodepasos,digamos30porejemplo,algunaspartculasentranyotrassalen,quedandoadistanciasmenoresquers.LoquepropusoVerletfuehacerunalista,paracadapartcula,detodoslosvecinosqueestandentrodesuesferaderadiors, yas envezdeexaminarladistancia de la partculai con todas lasN 1 restantes, se examinan esas distancias solo conlas partculas de su lista. Esta lista se contruye cada cierto n umero de pasos.El ahorro con este metodo es signicativo para sistemas de entre 500 a 5000 partculas, paraloscualeseltiempoporpasodesimulacionbajapracticamentealamitad(ver[9],pag.149).Para sistema de entre 100 a 200 partculas los cambios no son sustanciales, mientras que parasistemasdemasde5000partculassehanideadometodosmasecientes, comoel linkcelllist [9].4.2. DinamicaMolecularUna parte central de todo programa de DM lo constituye el algoritmo de integracion. Lasecuaciones de movimiento de Newton dadas por la Ec.( 6) son ecuaciones diferenciales ordinariasacopladas, nolineales, de segundo orden. Ellas deben ser resueltas numericamente. Dadas lasposiciones y velocidades iniciales a un tiempo inicial t0, la tarea del algoritmo es entregar lasposiciones y velocidades al tiempot0 + t.En lo que sigue revisaremos el formalismo basico de la DM en el ensemble microcanonico. Eneste caso las variables termodinamicas que se mantienen constantes son el n umero de partculasN,elvolumenV ylaenergainternaE.LaLagrangeanavienedadapor L = / 1,donde/ =

Ni=1(mi/2) q2ies la energa cinetica y 1es la energa potencial.La dinamica de este sistema esta gobernada por las ecuaciones de EulerLagrangeLqiddt_L qi_= 0 , i = 1, . . . , N (5)8 Generaciondelasconguracionesdonde qi y qi son las posiciones y las velocidades respectivamente. Esto da origen a las ecuacionesde movimientomiri = ri1, i = 1, . . . , N. (6)Cuando las fuerzas entre las partculas son conservativas la Hamiltoniana H es una constantede movimiento y la energa total se conserva:H = / +1 = E. (7)Cabe hacer notar que generalmente en Dinamica Molecular el momentum tambien es una can-tidad conservada cuando las paredes del recipiente son reemplazadas por condiciones de bordeperiodicas. Esto signica que existe una ligadura adicional, que se reejara en las propiedadestermodinamicas calculadas en la simulacion. Sin embargo, se ha demostrado que tales efectosson despreciables para sistemas de mas de cien partculas [19, 20].Exiten numerosos algoritmos para integrar las ecuaciones de Newton. Todos ellos conviertenlas ecuaciones diferenciales en ecuaciones de diferencias nitas. En DM la eleccion del algoritmoes (nuevamente) un compromiso entre el grado de precision requerido y el costo computacional.LosalgoritmosmasusadossoneldeVerlet[27],elvelocityVerletyelalgoritmodeBeeman.Otros algoritmos muy usados en DM son los del tipo correctorpredictor como el de Gear [9].Revisemosporejemploel algoritmodeVerlet. Paradeducirlo, partimosdel desarrolloenserie de Taylor de r(t),r(t + t) = r(t) +v(t)t + a(t)2(t)2+ (8)donde la aceleracion es a(t) = F(t)/m. Del mismo modo,r(t t) = r(t) v(t)t + a(t)2(t)2+ (9)Sumando ambos desarrollos, obtenemosr(t + t) +r(t t) = 2r(t) +a(t)(t)2+O((t)4) (10)As, la nueva posicion r buscada, en el tiempot + t, viene dada porr(t + t) 2r(t) r(t t)a(t)(t)2(11)El error estimado que contiene la nueva posicion r es del oreden de t4, donde t es el pasodetiempo(timestep)enlasimulaciondeDM. Notesequeparaevaluarlanuevaposicionr(t +t) solo necesitamos conocer la posicion anterior (en t t) y la aceleracion en el tiempot; no se necesita la velocidad. Sin embargo, esta la podemos calcular a partir der(t + t) r(t t) = 2v(t)t +O((t)3) (12)de dondev(t) =r(t + t) r(t t)2t+O((t)2) (13)Como se ve, el error en la velocidad es del orden t2, ademas que no se trata en el mismo nivelque la posicion. Un algoritmo que supera este hecho es el velocity Verlet, donde la posicion y lavelocidad se obtiene al mismo tiempo (t + t).9Notemos que la forma mostrada de conducir la dinamica, esto es, integrando las ecuacionesde movimiento y avanzado con un tiempo subdividido en porcionest es eciente para el casode potenciales continuos, pero si tenemos un potencial tipo esferas duras de radioro, o seaV (r) =_ sir < ro0 sir > roes mejor usar lo que se denomina dinamica dirigida por eventos. Supongamos que tenemos unconjunto de n esferas duras en un caja c ubica y damos a cada esfera una cierta velocidad inicial(temperatura).Entonceselproximoeventoqueocurriraserai)elchoquedeunaesferaconla pared, o ii) el choque entre dos esferas. Dado que sabemos el radio de las esferas, as comola velocidad y posicion inicial, se puede calcular exactamente cual sera el proximo evento y eltiempo te en el cual ocurrira. Entonces en vez de integrar las ecuaciones de movimiento, es masecientemoverdirectamentecadaunadelasesferasdesdesuposicionrialanuevaposicionri=vite. Tantoparael choquedeunaesferacontralaparedoelchoqueentredosesferasentre ellas, se puede suponer una colision elastica y por tanto solo cambiar el sentido y direccionde la velocidad seg un la ley de choques elasticos.En resumen, un programa de DM trabaja de la siguiente forma:1. Se leen los parametros que especican las condiciones de la corrida tales como la temper-aturainicial, el n umerodepartculas, laposiciondelaspartculas(estructurafcc, bcc,etc) la densidad, el paso de tiempo t, tiempo total de simulacion, etc.2. Se inicializa el sistema, esto es, se asignan las posiciones y las velocidades iniciales.3. Se calculan las fuerzas sobre todas las partculas.4. Seintegranlas ecuaciones demovimientodeNewton. Estepasoas comoel anteriorconforman el loop central de la simulacion. Ellos son repetidos hasta haber calculado laevoluciontemporal del sistemaduranteel tiempototal desimulaciondeseado. Sevanguardando las posiciones, velocidades, fuerzas, etc, durante cada paso en un archivo paraluego ser procesadas.5. Despues de haber completado lo anterior, se calculan y se imprimen los diferentes prome-dios relevantes. Entonces STOP.4.3. MonteCarloEl otro metodo usado para generar conguraciones independientes en el espacio de fases esel llamado Metodo de Montecarlo. Este es un metodo de caracter probabilstico o estocastico,que hace uso intensivo de un generador a n umeros aleatorios en su funcionamiento. Precisamentesu nombre es en honor a la ciudad de Montecarlo, famosa por sus casinos, ruletas y juegos deazar, entre otras diversiones.Supongamosquedeseamoscalcularunaciertapropiedad /deunsistemaenelensemblecanonico (NV T). Seg un la mecanica estadstica, el valor promedio de /totalesta dado por/)total =_dpNdrNexp [H(pN, rN)]/(pN, rN)_dpNdrNexp [H(pN, rN)](14)10 GeneraciondelasconguracionesEn general, como la parte de energa cinetica / del Hamiltoniano depende cuadraticamentedelosmomentap,estosgradosdelibertadsepuedenintegraranalticamente,quedandosolouna integral que depende de las coordenadasr/) =_drN/(rN)Z(rN) (15)conZ(rN) =exp[|(r)]_drNexp |(r)(16)Nuestro problema entonces se reduce a que dada la energa potencial |(rN), debemos calcularlafunciondeprobabilidad Z(rN)yluegohacerlaintegral. Paraellopodramosutilizarelmetodo tradicional de cuadraturas. Sin embargo, a poco andar uno se da cuenta que se trata deuna integral multidimensional (en 6Ndimensiones, donde Nes el n umero de partculas) dondelos metodos de cuadraturas son prohibitivos por la lentitud de ellos. Una alternativa a esto eshacer la integral por el Metodo de Montecarlo.Supongamos que por alg un procedimiento podemos generar aleatoriamente puntos del espa-cio de conguraciones de acuerdo a la distribucion de probabilidades Z(rN). Esto signica queen promedio, el n umero de puntos ni generados por unidad de volumen alrededor del punto rNies igual aLZ(rNi), dondeL es el n umero total de puntos generados. O sea/) =1LL

i=1niA(rNi) (17)Claramente generar aleatoriamente los puntos en elespacio de conguraciones no es el mejormetodo, pues como Z(rN) es proporcional al factor de Boltzmann exp(U), los puntos quesolo tienen baja energa contriburan signicativamente, mientras que los de alta energa tendranunbajopesorelativo.Laclaveentoncesestaenidearunmetodoquesologenerepuntosquetenganunpesorelativoalto3. EstofueprecisamenteloqueresolvioMetropolisenlosa noscincuenta, con el algoritmo que lleva su nombre.Acontinuaciondamosexposicionheursticadeestemetodo. Consideremoslatransicionentre estados a y b, donde la probabilidad que el sistema este en el estado a es pa y en el estadob espb. En equilibrio se cumple la condicion del balance detallado, y por tantoWabpa = Wbapb(18)dondeWabeslaprobabilidaddequeocurralatransiciona byWbaeslaprobabilidaddeque ocurra la transicion deb a. Por otro lado, sabemos que en equilibrio la probabilidad deencontrar al sistema en el estadoa(b) es proporcional a exp(Ea(b)), luegopapb= exp(Eab). (19)Reemplazando en ( 18) quedaWbaWab= exp(Eab) (20)As, vemos que tambien la probabilidad de transicion entre dos estados posibles esta controladapor el factor de Boltzmann. Esto nos da la clave de como escoger la conguracion (o sea, como3Estoesloquesellamaimportancesampling.11hacer el sampling) para tener siempre un peso relativo alto: si al ir del estadoa al estadob elcambiodeenergaEabesmenorquecero,entonceslaprobabilidadWabes1(decimosqueaceptamoslamovida).PerosielcambioEab> 0,entoncesWab exp Eab.Enestecaso comparamosWabcon un n umeror entre [0,1] escogido al azar, y siWab> r, se acepta lamovida, y si es menor se rechaza.Enlapractica,sitenemosunsistemadeNpartculassometidoaunatemperaturaT,elalgoritmo de Metropolis se implementa as:1. Seleccione una partculanicualquiera y calcule su energaUi(r)2. Dealapartculaundesplazamientoaleatorior

=r + rycalculesunuevaenergaU

i(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 posicion de la partculani, y vaya a 1.4. Luegodebarrer laNpartculas, guardelaconguracionobetnidayvuelvaa1paracomenzar otro ciclo. Despues de un n umero razonable de paso, vaya a 55. Calcule la propiedades fsicas de interes a partir de la conguraciones guardadas. STOP.Antes de terminar este apartado hacemos las siguientes obeservaciones:amenudoes costumbreusar, envezdelaprobabilidadexp(U), laprobabilidadnormalizada1/(1 + exp(U)).Estonocambialosresultadosdelequilibrio,ytienelaventaja que el metodo converje mas rapido para altas temperaturas.Enterminostecnicos, podemosdecirqueel algoritmodeMetropolisesunprocesodeMarkov en el cual se construye una caminata aleatoria (randomwalks) de tal modo quela probabilidad de visitar un punto particular rNes proporcional al factor de Boltzmannexp(U). De hecho hay varios maneras de construir tal caminata aleatoria. El algoritmode Metropolis es solo una de ellas, la mas usada.Unadelamayores diferencias entreMCyDMes queenMCnoes posiblecalcularpropiedadesdinamicas, puesladinamicaquesesigueescticia. Entoncescual eslaventaja de MC? La respuesta depende del tipo de problema a tratar. Seg un Frenkel [13],siempre que se pueda, debera preferirse DM. Sin embargo, hay algunos casos (que no sonpocos) en los cuales no es posible usar DM: i) sistemas que no poseen dinamica intrnseca,como modelos en la red, sistemas de espines (Ising, Heisenberg, modelo de Potts, etc). ii)casos en que las barreras de activacion son demasiado altas.5. AnalisisdelosresultadosLoquehemostratadohastaahoradicerelacionconlasimulacionpropiamentetal, estoes, comoobtenerlatrayectoriaenel espaciodefasesparaunsistemadepartculas. Ahoraveremos como analizar estas trayectorias para obtener de all propiedades fsicas macroscopicasque puedan ser comparadas con el experimento. Dividiremos estas propiedades en: propiedadestermodinamicas, propiedades estaticas y propiedades dinamicas.12 Evaluaciondelaspropiedadesfsicas5.1. PropiedadestermodinamicasEn el ensemble microcanonico la temperatura del sistema se calcula como el promedio de laenerga cinetica, a traves del teorema de equiparticion/) =32NkBT. (21)dondekB= 1,380621023J/K es la constante de Boltzmann. El . . .) se reere al promediosobre las Npartculas y sobre el ensemble. La presion media, p, se calcula a traves del teoremadel virialpV= NkBT +W) , (22)dondeW) =13_N

i=1ri Fi_(23)es conocido como el virial del sistema, siendo ri la posicion de la partcula i y Fi la fuerza sobrela partculai debido a todas las restantes.Apartirdelasuctuacionescuadraticasmedias, (/)2)= (A A))2),seobtienenlasfuncionesderespuesta, talescomocalorespecco, modulodecompresibilidadadiabaticoyconstantes elasticas. Por ejemplo, el calor especco por partcula,CV , se obtiene de(/)2)/)= (1)2)1)= NkBT_1 3kB2CV_. (24)Este tipo de relaciones entre una funcion respuesta y la desviacion cuadratica media de ciertacantidadesuncasoespecial del teoremadeuctuaciondisipacion. SuusoenDMofrecelaclaraventajaquemedianteunasolacorridadesimulacionpodemosobtenerinmediatamentelafuncionderespuestadeseada. Unmetodoalternativoparaobtener, porejemplo, el calorespecco, seracorrerel sistemaavariastemperaturasdiferentesygracarlaenergatotalrespecto de la temperatura. Su pendiente nos daraCV , con el evidente gasto de tiempo.Otras propiedades importantes a calcular son la entropa y las energas libres. Desgraciada-mente esto no es tan directo como lo mostrado arriba y requiere de tecnicas mas elaboradas [28].5.2. PropiedadesestaticasLaspropiedadesestructuralesestaticasdeunsistemasepuedendescribiratravesdelafuncion de distribucion de pares (FDP),g(r) y del factor de estructura estaticoS(q).La primera esta dada porg(r) = n(r, r + r))4r2rVN, (25)donde n(r, r+r) indica el n umero de partculas que hay en una capa entre r y r+r, teniendocomoorigenunadeterminadapartcula. Lafunciondedistribuciondeparesesproporcionalalaprobabilidaddeencontrardospartculasseparadasporunadistanciar + r.Escom untambiengracarlafunciondedistribucionradial FDR = 40g(r)r2;aquelareaencerradapor el primer pico es proporcional al n umero de coordinacion mientras que el cuociente entrela posicion del primer y segundo pico informa sobre las distancias interatomicas.Experimentalmente lo que se mide, mediante scattering de neutrones y rayos X, es el factorde estructura estatico,S(q). En el caso de lquidos y materiales amorfos, este solo depende delPropiedadesestaticas 13Figura 3: Funcion de distribucion de pares y patron de difraccionmodulo [q[ y puede representarse como una integral sobreg(r):S(q) = 1 + 4_R0r2[g(r) 1]sen(qr)qrdr , (26)donde el valor de R debe escogerse menor que la mitad de la longitud de la caja de simulacion.En la Figura 3 se muestra esquematicamente la forma tpica que presenta la funcion de dis-tribucion de pares y el patron de difraccion para gases, lquidos, amorfos y solidos, en referenciaa la distribucion de sus atomos en el espacio real.En el caso de solidos cristalinos una cantidad importante, que permite caracterizar el des-orden de las diferentes capas atomicas, es el factor de estructura estatico denido comoS(K, ) =_1N2

N

j=1exp(iK rj)2_, (27)dondeKes unvector delaredrecprocadeunacapaatomica, N

es el n umerototal departculas en la capa y el ndicejse reere a cada una de las partculas en esa capa.14 Evaluaciondelaspropiedadesfsicas5.3. PropiedadesdinamicasUna de la ventajas de DM con respecto a MC reside en la posiblidad que ofrece de calcularpropiedadesdinamicas.Comohemosvisto,enMCnohayunadinamicagenuinayportantonosepuedenobtenerdirectamentede elpropiedadesquedependandelaevoluciontemporaldel sistema. DM permite calcular estas en forma directa.Las propiedades dinamicas las obtenemos a partir de las funciones de correlacion temporal.Estas son muy importantes en simulacion computacional pues estan directamente relacionadascon los coecientes de transporte y por tanto, con los datos experimentales.El coecientededifusionD, porejemplo, lopodemoscalcularapartirdelafunciondeautocorrelacion de velocidadesZ(t) =13vi(t)vi(0))a traves deD =_0Z(t)dt . (28)Estees unejemplodelaformuladeGreenKubo, quepermiteexpresar uncoecientedetransporte (macroscopico) como una integral sobre el tiempo de una funcion de autocorrelaciontemporal (microscopica).Una manera alternativa de calcularlo es usando la relacion de EinsteinD =lmt[ri(t) ri(0)[2)6t, (29)donde [ri(t)ri(0)[2) es la desviacion cuadratica media de la posicion de la partcula i respectoa su posicion inicial ri(0).Otra cantidad que es posible calcular a partir de Z(t) es la densidad de estados vibracionales,G(),G() =Z()Z(t = 0),dondeZ() es la transformada de Fourier deZ(t) y esta dado porZ() =12_eitZ(t)dt . (30)A. ApendiceA.1. ComentarioBibliogracoSobre fsica computacional y simulacion en general:KauamnySmarr:librodedivulgacionsobrelarelevanciaquehatomadoytomaralasimulacion computacional en practicamente todas las areas de la ciencia.ComputinginScience andEngineering, revistabimensual, herederade Computer inPhysics, permite estar al da del estado del arte en este tema.Gould y Tobochnik: abarca metodos numericos y tecnicas de simulacion usadas en fsica.Muy pedagogico y con gran cantidad de utiles ejemplos. Altamente recomendable.Propiedadesdinamicas 15Koonin y Meredith: libro muy utilizado, que pone enfasis en lo metodos numericos masusados hoy da en fsica, con su respectiva implementacion computacional.Sobre MD:AllenyTildesley: excelentelibroparainiciarseenlos secretos delaDMyMC, condiscusionesprecisasyejemplosdesuimplementacioncomputacional, incluidosmuchosprogramas de computacion.Haile: mas que recetas, hay un tratamiento detallado de los fundamentos de la dinamicamolecular, con muchos ejemplos y problemas. Buen complemento de Allen y Tildesley.Frenkel y Smit: en la misma lnea que el de Allen y Tildesley, es un complemento obligadoportraertemasmasactualesydarexplicacionesmuyclarasdelosdistintosmetodosusados. La discusion sobre equilibrios de fase es notable. Ademas esta muy bien escrito yllevado lo que hace que la lectura sea muy entretenida. Excelente.RecularyMadden:unpaperclasicodondeseexplicademaneralegibleelmetododedeCar-ParrinellosobreDMab-initio. Serequiereconocimientos basicos deMecanicaCuantica y DM clasica. Altamente recomendable.Sobre MCAllenyTildesley: excelentelibroparainiciarseenlos secretos delaDMyMC, condiscusionesprecisasyejemplosdesuimplementacioncomputacional, incluidosmuchosprogramas de computacion.Frenkel y Smit: en la misma lnea que el de Allen y Tildesley, es un complemento obligadoportraertemasmasactualesydarexplicacionesmuyclarasdelosdistintosmetodosusados. Excelente.Binder: referencia obligada para cualquier interesado en profundizar sobre el tema.A.2. DireccionesenInternethttp://www.sissa.it/furio/md: paginadeFurioErcolessi, claveenMD. All hayunosexcelentesapuntesdeMDconejemplosenFortran90ylinksmuyinteresantes. EnlaSISSAyenel ICTP, enTrieste, Italia, fuedondenacioel metododeMDabinitiodeCarParrinello. La pagina contiene los distintos grupos que alli hay: estructura electronica,supercies, 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, tantolibrecomocomercial, enDM, MC, DMab-initio, etc. EnDMsonrecomendables losprogramas MOLDY y DL POLY.http://www.ncsa.uiuc.edu/Apps/CMP/ceperley : David Ceperley es uno de los expertosmundialesenMCysuaplicacionenMateriaCondensada. Enestesitiohayapuntes,ejemplo de software, e interesantes links.16 EvaluaciondelaspropiedadesfsicasA.3. SimulacioncomputacionalenChileA continuacion damos una lista (ciertamente incompleta e imprecisa) de personas que tra-bajan en simulacion computacional (o la usan como herramienta auxiliar) en nuestro pas.Universidad Catolica de Antofagasta: Sergio Curilef utiliza DM para estudiar la estadsticade Tsallis.Universidad de Chile:Facultad de Ciencias Fsicas y Matematicas: existe un buen grupo formado por Patri-cio Cordero, Rodrigo Soto, Rosa Jimenez y J.M Pasini. Han usado DM dirigida porlos eventos para estudiar gases y lquidos. Verhttp://www.tamarugo.uchile.cl/Facultad de Ciencias: Jose Rogan y Rodrigo Ferrer estan usando MC para estudiarsistemas magneticos. Jaime Roessler tiene experiencia en MC.Universidad de Santiago: Carlos Esparza ha trabajado desarrollando nuevos metodos deDM. Gonzalo Gutierrez usa DM para caracterizar materiales lquidos y amorfos, as comosupercies e interfaces de solidos. Dora Altbir usa MC para estudiar sistemas magneticos.Guillermo Palma y Teodoro Meyer tienen experiencia con MC en la red.P. Universidad Catolica: Miguel Kiwi y Ricardo Ramrez usan DM para estudiar solidos.Aldo Romero tiene experiencia en DM ab-initio, CarParrinello.Universidaddel BioBio: DinoRissoinvestigaenhidrodinamicausandoDMdeesferasduras.Universidad de La Frontera: el grupo de fsica del solido (Eugenio Vogel et al) tiene vastaexperiencia en MC para sistemas de espines.Universidad de Magallanes: Mauricio Marin ha desarrollado algoritmos ecientes en DM.Referencias[1] William J. Kaufmann y Larry L. Smarr,SupercomputingandtheTransformationofSci-ence, Scientic American Library, New York, 1993.[2] New dimensions in simulation, Special issue, en Physics World 9, No. 7, p.2948 (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: Applica-tions to Physical Systems, second edition, , AddisonWesley, Reading MA (1996).[5] S. E. Koonin and D. C. Meredith, Computational Physics, AddisonWesley, Reading MA(1990).[6] K. Binder,MontercarloMethodinStatisticalPhysics, Springer, Berlin, 1986.[7] D. Ceperley, Apuntes en pagina WEB, http://www.ncsa.uiuc.edu/Apps/CMP/ceperleyReferencias 17[8] D.W. Heerman, ComputerSimulationMethodsinTheoretical Physics,Springer-Verlag,1986.[9] M. P. Allen and D. Tildesley,ComputerSimulationsofLiquids, Clarendon Press, Oxford,1987.[10] J. P. Hansen, An introduction to Molecular Dynamics with an applicaction to glass transi-tion, enComputerSimulationinMaterialsScience, Edited by M. Meyer and V. Pontikis,Kluwer Academic Publishers, 1991.[11] J. H. Haile,MolecularDynamicsSimulation, J. Wiley, New York, 1992.[12] D. E. Rappaport, TheArt of Molecular Dynamics Simulation, CambridgeUniv. Press,1996.[13] D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, San Diego,1996.[14] SimulationofLiquidsandSolids, edited by G. Cicotti, D. Frenkel, and I. R. Mc Donald,North-Holland, Amsterdam, 1987.[15] D. Frenkel, IntroductiontoComputerSimulation, enSimpleMolecularSystemsatveryHighDensity, 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 artof scientic computing, Second Edition, Cambridge University Press, 1992.[22] M. Gillan, Contemporary Physics 38, No.2, p.115130 (1997).[23] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).[24] D. K. Reculer and P. A. Madden, Molecular dynamics without eective potentials via theCarParrinello approach, Mol. Phys. 70, 921 (1990).[25] G. Galli and A. Pasquarello M. Parrinello, First Principles Molecular Dynamics, enCom-puterSimulationinChemicalPhysics, Edited by M. P. Allen and D. J. Tildesley, KluwerAcademic Publishers, 1993.[26] MRS Bulletin, Vol. 21, No.2, Feb. 1996.[27] L. Verlet, Phys. Rev. 159, 98 (1967).[28] D. Frenkel, LecturesNotesonFreeEnergyCalculations, enComputerSimulationinMa-terialsScience, Edited by M. Meyer and V. Pontikis, Kluwer Academic Publishers, 1991.