Moleculaire Dynamica en Monte Carlo Simulaties Case Study 17 Solid-Liquid Equilibrium of Hard Spheres Joost van Bruggen 0123226 6 juli 2004
1
Inhoudsopgave 1 Thermaliseren
2
2 Waarde van λmax
2
3 Integreren
2
4 het Berekenen
4 Samenvatting
In dit verslag worden de resultaten van het werken aan Case Study 17 uit het boek van Frenkel en Smit [2] besproken. Er wordt tijdens de simulaties gebruik gemaakt van een fcc rooster met daarop verschillende aantallen deeltjes. De thermalisatietijd bedraagt 1000 Monte Carlo stappen. De waarde van λmax is 2000, c heeft de waarde 70 en de parameter ndispl heeft steeds de waarde 1000. Voor het berekenen van de integraal wordt gebruik gemaakt van een 10 punts Gaussische kwadratuur in combinatie met een coordinatentransformatie. Bovendien wordt een extrapolatie gemaakt naar oneindige systeem grootte, omdat de verkregen antwoorden van de grootte af blijken te hangen en we op zoek zijn naar een waarde voor macroscopische systemen.
1
Thermaliseren
Om te kijken hoe lang de thermalisatietijden dienen te zijn, werd voor verschillende aantallen thermalisatiestappen gekeken hoe de uitkomst afhangt van die aantallen. Door steeds de thermalisatietijden te verhogen en te kijken naar de waarde van de gemiddelde kwadratische verplaatsing, werd gezocht naar een limiet waarde. Het verschil tussen de waarde van de gemiddelde kwadratische verplaatsing tussen opeenvolgende geprobeerde waarden voor de thermalisatietijd, werd steeds kleiner. Er is op deze manier gekozen voor een thermalisatietijd van 1000 Monte Carlo stappen. Dit lijkt een goede afweging tussen snelheid en nauwkeurigheid te zijn.
2
Waarde van λmax
De waarde van λmax is op advies van het boek gekozen op 2000. In figuur 1 is te zien dat de met simulaties verkregen curve (zie notebooks op [1] voor details) zeker voor λ > λmax = 2000 voldoet aan de analytische oplossing van een Einstein solid. Hiermee kunnen we dus concluderen dat de gebruikte λmax goed gekozen is.
3
Integreren
Om de resultaten uit de simulatie te verwerken, is, wederom op aanraden van Frenkel en Smit [2], gekozen voor een 10 punts Gaussische kwa-
2
rood : simulatie, groen : Einstein solid 0.014 0.012 0.01 0.008 0.006 0.004 0.002 Λ 0.001
0.01
0.1
1
10
100
1000
Figuur 1: Vergelijking van de met simulaties verkregen curve met de voor het Einstein solid bekende analytische oplossing
dratuur methode. Eveneens is de voorgestelde coordinatentransformatie toegepast.
de coordinatentransformatie en Gaussische kwadratuur De voorgestelde coordinatentransformatie dient er voor om de functie gladder te laten verlopen zodat de kwadratuur methode een nauwkeuriger antwoord oplevert. Deze methode is namelijk geschikt om functies te integreren die benaderd kunnen worden door polynomen van graad tot 2 × p met p het aantal evaluatiepunten in de Gaussische kwadratuur. Volgens de kwadratuur methode geldt:
Z
b
f (x)dx = a
N X
wi f (xi ),
i
met N het aantal evaluatiepunten en f een willekeurige functie die, wil men een nauwkeurig antwoord krijgen, goed benaderd kan worden door een polynoom van graad maximaal 2N .
waarde van c In de coordinatentransformatie komt een parameter c voor die in de uitgevoerde experimenten steeds de waarde 70 had.
3
4
het Berekenen
In deze sectie wordt stap voor stap het gevolgde traject tot het eindantwoord beschreven. Na het vaststellen van waarden voor alle parameters (zoals λmax en c), is er voor systemen van verschillende grootte gesimuleerd. Zie tabel 1 voor de gebruikte waarden. Dit simuleren is gedaan door de bij elke systeem grootte voor 10 verschillende waarden van λ te runnen. Die 10 waarden van λ zijn verkregen door gebruik te maken van Mathematica’s GaussianQuadratureWeights, die ook meteen de gewichten wi uit de hierboven genoemde som oplevert. Zie voor de details weer de notebooks op [1]. Hierna is de som uitgerekend alsmede de fout (in het kwadraat) door kwadratische optelling van de uit de simulatie verkregen fouten. Als laatste zijn de door Frenkel en Smit [2] voorgestelde correcties toegepast. De verkregen resultaten zijn te vinden in tabel 1.
systeem grootte 3*3*6 = 54 3*6*6 = 108 6*6*6 = 216 6*6*12 = 432 9*9*6 = 486
(βF ex + Ln(N ))/N 5.8812 5.8738 5.8996 5.8954 5.9002
fout 0.0005484 0.0005731 0.0003995 0.0001804 0.0002208
Tabel 1: Resultaten van simulaties Omdat de te berekenen excess free energy enigszins afhangt van de grootte van het systeem, zoals blijkt uit de waarden in tabel 1, is een extrapolatie gemaakt naar oneindig grote systemen. De hiermee verkregen waarde voor de excess free energy bedraagt: 5.899(2).
4
Referenties [1] http://www.phys.uu.nl/~vbrug/vakken/mosim [2] Daan Frenkel en Berend Smit, Understanding Molecular Simulation, Academic Press (2002)
5
