Simulari Monte Carlo in 'ansamblul Gibbs'

Simulari Monte Carlo in 'ansamblul Gibbs'

Asa cum a fost argumentat mai inainte, simularile pe calculator sunt foarte asemanatoare cu experimentele din lumea reala, motiv pentru care simularile pe calculator sunt plasate in domeniul stiintelor experimentale.

Aceste asemanari se opresc insa in momentul in care se doreste sa se studieze transformarile de faza, in special cele de speta intai. Intr-un experiment real, o tranzitie de faza de speta intai este relativ usor de localizat: la temperatura si presiunea potrivite se observa ca un sistem initial omogen se va separa in doua faze distincte separate printr-o interfata, proprietatile fizice ale acestor faze putand fi atunci determinate prin masuratori directe.

In cazul simularilor pe calculator, ar trebui sa se determine proprietatile termodinamice ale fazelor individuale, apoi sa se gaseasca punctele in care temperaturile, presiunile si potentialele chimice ale acestora sunt egale. Aceasta metoda, in principiu corecta, sufera insa din cauza necesitatii de a se efectua un mare numar de simulari pe calculator la diversi parametri de lucru, care apoi necesita o examinare atenta pentru a gasi coexistenta celor doua faze.

Situatia este ingreunata si de faptul ca in preajma tranzitiilor de faza sistemul prezinta mari fluctuatii ale parametrilor sai. Un exemplu elocvent in acest sens il reprezinta graficul urmator in care s-a obtinut prin simulare pe un sistem Ising bidimensional magnetizarea ca functie de temperatura.

Explicatia principala a acestor dificultati consta in numarul mic de particule cu care suntem obligati sa lucram, caz in care o mare parte din molecule se afla de fapt la interfata sistemului.

Pentru a se evita situatia in care sa fie nevoie sa se faca un numar foarte mare de simulari pentru a studia coxistenta unui sistem, a fost pusa la punct o metoda speciala de simulare, numita Metoda 'Ansamblului Gibbs', si care permite sa se studieze coexistenta a doua faze intr-o singura simulare.

Metoda 'Ansamblului Gibbs'

Conditia de coexistenta a doua sau mai multe faze este egalitatea parametrilor intensivi ai acestora:

Egalitatea presiunilor:

Egalitatea temperaturilor:

Egalitatea potentialelor chimice:

La prima vedere, s-ar parea ca ansamblul in care ar rebui sa se studieze coexistenta unui astfel de sistem ar trebui sa fie ansamblul 'μPT'. In realitate un astfel de ansamblu nu exista, din motiv ca, prin fixarea tuturor parametrilor intensivi, toti parametrii extensivi sunt fluctuanti si nemarginiti. Conditia ca un ansamblu statistic sa existe este ca acesta sa aiba cel putin un parametru extensiv fixat, cum este cazul ansamblelor canonic (NVT), izoterm izobar (NPT), grand-canonic (μVT) sau microcanonic (ENV).

Metoda de simulare Monte Carlo in cazul 'ansamblului Gibbs' se bazeaza pe observatia ca doua faze aflate in conditii de echilibru realizeaza un schimb de particule una cu cealalta, ceea ce conduce la realizarea egalitatiii potentialelor chimice, iar ceilalti parametrii intensivi sunt egali fie prin fixarea lor de catre   rezervoarele cu care sistemul se afla in contact, fie se incearca egalizarea acestora printr-o alta metoda.

Motivul pentru care un astfel de ansamblu functioneaza este acela ca, la echilibrul de faza, presiunile si temperaturile sunt egale si pot fi determinate prin tehnici standard de simulare, dar valorile absolute ale potentialelor chimice nu se pot de fapt determina, multumindu-ne sa afirmam ca la echilibru diferenta dintre potentialele chimice ale celor doua faze este zero:

Simularile Monte Carlo in cazul 'ansamblului Gibbs' se fac pornind de la un ansamblu canonic avand parametrii N, V si T, dar care este impartit in doua prin intermediul unui piston mobil.

Fie N₁ si N₂ numerele de particule din cele doua compartimente cu , iar V₁ si V₂ volumele acestor compartimente cu

Existenta pistonului mobil asigura ca presiunile in cele doua compartimente sunt egale , iar temperaturile in cele doua compartimente sunt fixate de catre mediul exterior:

Daca se considera ca intre cele doua compartimente poate avea loc un schimb de particule, acest fapt asigura ca la echilibru si potentialele chimice ale sistemului sunt egale:

Partea configurationala a functiei de partitie a unui astfel de sistem se obtine din functia de partitie ansamblului canonic:

Aceasta forma de scriere a functiei de partitie asigura urmatoarele aspecte:

• numarul de particule din fiecare compartiment poate varia intre 0 si N;
• volumul fiecarui compartiment poate varia intre 0 si V;

Datorita faptului ca volumele celor doua compartimente sunt fluctuante, in calculul termenilor de interactie ai functiei de partitie, s-au folosit coordonate scalate, similar cu cele folosite in cazul ansamblului izoterm-izobar, cu diferenta notabila ca doua molecule interactioneaza doar daca fac parte din acelasi compartiment.

Realizarea unei simulari Monte Carlo

in cazul 'ansamblului Gibbs'

Particularitatile fizice ale sistemului studiat conduc la ideea ca in cazul ansamblului Gibbs se pot face urmatoarele tipuri de modificari:

1. Deplasari ale moleculelor.

Aceasta deplasare este identica cu deplasarile din cazul ansamblelor studiate pana acum, cu observatia ca deplasarea unei molecule se poate face doar in cadrul compartimentului din care face parte aceasta (si cu respectarea conditiilor pe frontiera periodice) . Trebuie precizat si faptul ca, datorita necesitatii de a se respecta principiul reversibilitatii microscopice, alegerea compartimentului in care se face deplasarea moleculei trebuie sa fie facuta aleator.

Pentru a vedea care este schema de acceptare a unei astfel de deplasari in cazul ansamblului Gibbs, si in acest caz se pleaca de la schema Metropolis scrisa in densitati de probabilitate:

sistemul trece cu siguranta din starea 'OLD' in starea 'NEW', daca densitatea de probabilitate a starii 'NEW' este mai mare decat densitatea de probabilitate a starii 'OLD', iar daca nu, el poate trece totusi, cu o probabilitate egala cu raportul densitatilor de probabilitate ale celor doua stari.

Din expresia functiei de partitie pentru 'ansamblul Gibbs', densitatea de probabilitate este proportionala cu:

Presupunem ca se deplaseaza o molecula in compartimentul 1, in urma deplasarii ei, setul de coordonate scalate corespunzator acestuia devenind cu densitatea de probabilitate corespunzatoare:

Raportul densitatilor de probabilitate este in acest caz:

identic cu raportul din cazul ansamblului canonic

Prin acelasi rationament ca si in cazul celorlalte ansamble studiate, se gaseste ca raportul densitatilor de probabilitate depinde doar de energia de interactie a particului alese cu celelalte din compartimentul sau, inainte si dupa deplasare:

T Schema dupa care se face acceptarea unei deplasari de molecula este identica cu cea de la ansamblul canonic:

• Se alege la intamplare unul din cele doua compartimente;
• Se alege la intamplare o molecula din acest compartiment;
• se calculeaza marimea
• se genereaza un numar aleator ξ in intervalul (0, 1). Deoarece evenimentul de generare al numarului aleator si evenimentul ca sistemul sa aiba raportul densitatilor de probabilitate egal cu Arg sunt evenimente independente, probabilitatea ca numarul aleator ξ sa fie mai mare decat Arg este egala cu probabilitatea ca ξ sa fie mai mic decat Arg.
• Daca , atunci deplasarea moleculei este acceptata;
• Daca , deplasarea moleculei nu este acceptata.

2. Modificari ale volumului

Existenta pistonului mobil asigura faptul ca presiunile din cele doua compartimente sunt egale. In realitate, acest piston nu exista fizic in sistem, prezenta sa fiind in fapt asigurata prin cerinta ca

Sa presupunem ca volumul compartimentului 1 s-a modificat cu , astfel incat:

Cerinta de mai sus este asigurata daca in acelasi timp avem:

In starea veche, densitatea de probabilitate a sistemului este:

iar in starea noua este:

Trebuie tinut cont ca, prin modificarea volumelor, distantele dintre particulele ambelor sisteme se modifica, in acelasi mod ca in cazul ansamblului izoterm-izobar, de unde se obtine o modificare a energiei totale de interactie din fiecare compartiment.

Raportul densitatilor de probabilitate este:

Asa cum este usor de inteles, aceasta etapa a simularii, pe langa complicatiile deosebite ce le implica din punct de vedere al programarii, este si extrem de costisitoare din punctul de vedere al timpului de rulare.

O anumita eficientizare se poate realiza in felul urmator. Functia de partitie a sistemului studiat se mai poate scrie sub forma echivalenta:

iar densitatea de probabilitate este proportionala cu:

Se constata ca aceasta forma echivalenta a functiei de partitie nu modifica decat partea corespunzatoare volumelor, criteriile de acceptare pentru deplasarile moleculelor sau schimbul de particule ramanand neschimbate.

Motivul care sta la baza acestei alegeri este acela ca prin noua forma a functiei de partitie, in loc sa se faca o problema de tipul mersului la intamplare in V₁, se face o problema de tip random walk in , care se poate arata ca exploreaza mai eficient toate valorile posibile ale volumelor celor doua compartimente. Aceasta modificare se poate face si in cazul functiei de partitie a ansamblului izoterm izobar, pentru eficientizarea modificarilor de volum ale acestuia.

Raportul densitatilor de probabilitate devine in acest caz:

Din punct de vedere practic, acceptarea modificarii volumelor compartimentelor se realizeaza in modul urmator:

Presupunem ca la un moment dat cele N particule identice ce alcatuiesc sistemul se afla distribuite cate N1 si N2 in cele doua compartimente.

Presupunem ca la acest moment, volumul total al sistemului, VOL, este impartit in doua cutii de forma cubica, de volume VOld1 si VOld2, cu laturile BoxOld1 si BoxOld2. Fie EOld1 si EOld2 energiile totale ale sistemelor in aceasta situatie;

• se modifica aleator volumele sistemelor, dupa schema:

se genereaza numarul aleator

LnDV = Ln ( VOld1 / VOld2 ) + ( RND( ) - 0.5 )*DVMAX

se modifica volumele celor doua compartimente

VNew1 = Vol * Exp ( LnDV ) / ( 1 + Exp ( LnDV ) )

VNew2 = Vol - VNew1

se calculeaza noile laturi ale cutiilor de simulare:

BoxNew1 = (VNew1) ** ( 1 / 3 )

BoxNew2 = (VNew2) ** ( 1 / 3 )

• se rescaleaza pozitiile tuturor moleculelor din cele doua compartimente:
• XNEW( I ) = X( I ) * BoxNew / BoxOld, unde BoxNew si BoxOld sunt laturile cutiei de simulare in care se afla particula respectiva;

• se calculeaza energiile sistemelor din cele doua compartimente pentru noile pozitii ale particulelor, ENew1 si ENew2;

• se calculeaza marimile

Arg1 = - Beta * ( ( ENew1 - EOld1 ) + ( N1 + 1) * Ln ( VNew1 / VOld1 ) / Beta )

Arg2 = - Beta * ( ( ENew2 - EOld2 ) + ( N2 + 1) * Ln ( VNew2 / VOld2 ) / Beta )

Arg = Exp ( Arg1 + Arg2 )

• se genereaza un numar aleator ξ in intervalul (0, 1).
• Daca , atunci modificarea volumului este acceptata;
• Daca , modificarea volumului nu este acceptata;

3. Schimbarea unei molecule dintr-o cutie in alta

Pentru fixarea ideilor, presupunem ca va fi scoasa o molecula din compartimentul 1 si va introdusa in compartimentul 2.

Densitatea de probabilitate a starii vechi este:

Se foloseste aceasta forma a densitatii de probabilitate, deoarece in modificarile volumului s-a folosit mersul la intamplare in in locul celui in volum.

Densitatea de probabilitate a starii noi este:

unde este energia compartimentului 1, care acum ar urma sa aiba particule, iar este energia compartimentului 2, care acum ar urma sa aiba particule.

Raportul densitatilor de probabilitate este:

Se constata ca aparent si in acest caz este necesar sa se calculeze energia totala a ambelor compartimente inainte si dupa schimbarea cutiei in care se afla particula, operatie extrem de costisitoare. In aceasta situatie se poate face un rationament analog cu cel de la ansamblul grand-canonic:

Energia compartimentului cu N₁ particule poate fi scrisa ca energia sistemului cu sistemului cu particule, , plus energia de interactie a particulei ce urmeaza sa paraseasca sistemul, (pentru simplitate, am presupus ca urmeaza sa paraseasca sistemul a N₁ - a particula):

Aceasta particula devine a - a particula in compartimentul 2. Energia compartimentului 2 poate fi scrisa ca energia sistemului cu sistemului cu particule, , plus energia de interactie a particulei ce urmeaza sa intre in sistem, :

Cu aceste consideratii, raportul densitatilor de probabilitate devine:

depinzand de diferenta dintre energia pe care o va avea particula in noua cutie si energia pe care aceasta o avea in vechea cutie.

Realizarea practica a schimbului de particule intre cele doua cutii se realizeaza in modul urmator:

se alege la intamplare compartimentul din care va fi scoasa o molecula;

se alege la intamplare o molecula din acest compartiment si se calculeaza energia ei EOld;

se alege aleator pozitia in care urmeaza ca particula sa fie introdusa in noul compartiment;

se calculeaza marimea:

Arg = NBox ( Out ) * Vol ( In ) / ( ( NBox ( In ) + 1 ) * Vol ( Out ) ) * Exp ( Beta * ( EOld   - ENew ) ),

unde EOld este energia particulei ce urmeaza sa fie scoasa din compartimentul 'Out', iar ENew este energia pe care ar urma sa o aiba aceasta in compartimentul 'In'.

• se genereaza un numar aleator ξ in intervalul (0, 1).
• Daca , atunci schimbarea particulei este acceptata;
• Daca , schimbarea particulei nu este acceptata;

Program de simulare Monte Carlo

in ansamblul 'Gibbs'

Pentru exemplificare, consideram cazul unui lichid monoatomic a carui simulare Monte Carlo se realizeaza in ansamblul grand-canonic.

Parametrii sistemului supus simularii sunt:

Integer, Parameter :: NPart = 100 !Numarul total de particule (Constant).

Real, Parameter :: TStar = 1.2 !Temperatura

Real, Parameter :: BetaMC = 1. / TStar

Real DRMax !Deplasarea maxima in MC conventional. Initializat in Initializare

Real DVMax !Variatia maxima a volumului in MC de volum. Initializat in Initializare

Diferenta notabila dintre acest program si cele folosite pana acum este aceea ca lucram cu doua cutii de simulare, motiv pentru care marimile legate de dimensiunile fizice ale sistemelor si numerelor de particule sunt definite ca matrici (indicele 1 corespunde cutiei 1, iar indicele 2 corespunde cutiei 2):

Real, Dimension(1:2) :: Box !Laturile cutiilor de simulare (Fluctuante)

Real, Dimension(1:2) :: Vol !Volumele cutiilor de simulare (Fluctuante)

Real VTotal  !Volumul total al sistemului: Vol(1) + Vol(2) (Constant)

Real, Dimension(1:2) :: Rho !Densitatile in cele doua cutii (fluctuante)

Integer, Dimension(1:2) :: NBox !Numerele de particule in fiecare cutie (fluctuante) N(1) + N(2) = NPart

Deoarece numarul de particule din cele doua cutii fluctueaza in timpul simularii, pentru a se evita ca in cursul acesteia sa se redimensionexe matricile ce reprezinta particulele din cele 2 cutii (ca in cazul ansamblului Grand-Canonic), o solutie este ca in momentul in care se defineste structura ce reprezinta o molecula, sa se mai foloseasca o componenta ce va indica in ce cutie se afla particula respectiva:

Type Vector

Real X, Y, Z

Integer Cutie !Cutia in care se afla molecula respectiva

End Type

Moleculele ce alcatuiesc sistemul sunt o singura matrice de dimensiune fixata, apartenenta unei molecule la una din cutiile de simulare se stabileste prin intermedul componentei Cutie a vectorului.

Type (Vector), Dimension(1:NPart) :: Molecula !Moleculele ce alcatuiesc sistemul

Integer, Parameter :: NEquil = 10**3 !Nr de cicli MC de echilibrare.

Integer, Parameter :: NProd = 2*10**3 !Nr de cicli MC de productie.

Integer, Parameter :: NSamp = 10**2 ! Frecventa de mediere

Integer, Parameter :: NVol = NPart/10 !Nr de incercari de schimbare a volumului

Integer, Parameter :: NSwap = NPart/10 !Nr de incercari de schimbare a unei molecule dintr-o cutie in alta

Marimea CycleLength se alege astfel incat lungimea unui ciclu sa fie egala cu numarul de particule ce afla in sistem, NPart, plus un numar de modificari ale volumelor, NVol, plus un schimb de particule intre sisteme, NSwap. Aceasta alegere a lungimii unui ciclu are la baza necesitatea asigurarii principiului reversibilitatii microscopice: intr-un ciclu, trebuie ca fiecare molecula sa incerce in medie sa se deplaseze o data si sa se realizeze un numar predefinit de schimburi de particule intre compartimente si de modificari ale volumelor, alegerea acestor miscari facandu-se aleator.

Integer, Parameter :: CycleLength = NPart + NVol + NSwap !Lungimea unui ciclu

Integer :: NumarMedii = 0 !Nr de medii efectuate

Etapa de initializare a programului urmareste determinarea volumelor sistemului, laturile cutiilor de simulare si a deplasarii maxime in decursul unei deplasari Monte Carlo conventionale. In aceasta etapa plaseaza si particulele din cele doua cutii in retele cubice obisnuite:

Subroutine InitializareSistem()

Integer i

NBox = NPart/2 !Initial in cele doua cutii nr de particule este acelasi

Rho = 0.5 !Initial in cele doua cutii densitatea este aceeasi

Vol = Real(NBox) / Rho !Volumele celor 2 cutii de simulare (initial identice

VTotal = Vol(1) + Vol(2) !Volumul total al sistemului

Box = Vol ** (1.0 / 3.0) !Laturile cutiilor de simulare

DRMax = Box(1) / 20.0 !Initial cele doua cutii sunt egale

DVMax = Vol(1) / 20.0 !Initial cele doua cutii sunt egale

Call Lattice(1) !Plasam in latice moleculele din cutia 1

!In cutia 2 pozitiile sunt identice, dar Cutie = 2

Do i = 1, NPart/2

Molecula(NBox(1)+i) = Molecula(i)

Molecula(NBox(1)+i)%Cutie = 2

End Do

End Subroutine InitializareSistem

Rutina Lattice este identica cu rutinele de plasare a moleculelor in cutia de simulare din cazul ansamblelor canonic, izoterm-izobar sau grand canonic, cu diferenta ca aceasta trebuie sa precizeze si cutia in care se afla particulele plasate.

Programul de simulare Monte Carlo este urmatorul:

Program Gibbs_MC

Call InitializareSistem()

!Echilibrare

Do ICycle = 1,NEquil

Do IMC = 1,CycleLength

Call GibbsMonteCarlo()

End Do

End Do

!Productie

Do ICycle = 1,NProd

Do IMC = 1,CycleLength

Call GibbsMonteCarlo()

End Do

If ( Mod (ICycle, NSamp) == 0) Then

Call Sample() ! Acumulare medii

End If

End Do

End Program Gibbs_MC

Subrutina GibbsMonteCarlo este:

Subroutine GibbsMonteCarlo()

Real Decid !Folosit la luarea deciziei pentru ce tip de miscare se opteaza

Decid = Int(Rnunf() * CycleLength) + 1

If (Decid .LE. NPart) Then

Call MC_Move() !Deplasare a unei molecule in cutia ei

ElseIf (Decid .LE. (NPart +NVol)) Then

Call MC_Vol() !Modificarea volumelor cutiilor

Else

Call MC_Swap() !Schimb de particule intre cutii

End If

End Subroutine GibbsMonteCarlo

Subrutina MC_Move este o rutina clasica de deplasare a unei   molecule dintr-un anumit compartiment. Singura diferenta comparativ cu rutinele folosite pana acum este la nivelul precizarii corespunzatoare a parametrilor compartimentului:

Subroutine MC_Move()

Integer i !Indexul unei particule alese la intamplare

Real EOld, ENew !Energiile inainte si dupa deplasare

Type(Vector) CurrentParticle !Particula curenta

Type(Vector) TestParticle !Particula test

i = Int(RNUNF()*NPart) + 1 !Aleg o molecula la intamplare

CurrentParticle = Molecula(i)

EOld = Energie_I(CurrentParticle, i) !Calculeaza energia initiala a particulei ce se deplaseaza

!NOILE POZITII ALE PARTICULEI I

TestParticle%X = CurrentParticle%X + (RNUNF() - 0.5) * DRMax

TestParticle%Y = CurrentParticle%Y + (RNUNF() - 0.5) * DRMax

TestParticle%Z = CurrentParticle%Z + (RNUNF() - 0.5) * DRMax

TestParticle%Cutie = CurrentParticle%Cutie !!!Particula test se afla in aceeasi cutie cu particula curenta

!Consideratii periodice

TestParticle%X = TestParticle%X - ANINT ( TestParticle%X / Box(TestParticle%Cutie) ) * Box(TestParticle%Cutie)

TestParticle%Y = TestParticle%Y - ANINT ( TestParticle%Y / Box(TestParticle%Cutie) ) * Box(TestParticle%Cutie)

TestParticle%Z = TestParticle%Z - ANINT ( TestParticle%Z / Box(TestParticle%Cutie) ) * Box(TestParticle%Cutie)

ENew = Energie_I(TestParticle, i) !Calculeaza energia noua a particulei ce se deplaseaza

If (RNUNF() .LT. EXP(- BetaMC * (ENew - EOld))) THEN

Molecula(i) = TestParticle !Criteriul Metropolis clasic 

EndIf

End Subroutine MC_Move

Energia particulei I inainte si dupa deplasare se face in rutina EnergieI, foarte asemanatoare cu cele folosite pana in prezent. Aceasta rutina trebuie sa primeasca parametrii particulei a carei energie se doreste sa fie determinata, precum si indicele ei, pentru a evita calculul energiei de interactie cu ea insasi.

Real Function Energie_I(Arg, i)

!Calculeaza energia de interactie a moleculei I cu toate celelalte din cutia sa

Type(Vector), Intent(IN) :: Arg

Integer, Intent(IN) :: i !Indexul particulei I

Type(Vector) :: CurrentParticle !Particula curenta

Real RXIJ, RYIJ, RZIJ, RIJSQ, E

Integer j

E = 0.0

Do j = 1, NPart !Cicleaza pe toate particulele din sistemul dublu

If (Arg%Cutie .EQ. Molecula(j)%Cutie) Then !Se asigura ca lucram in aceeasi cutie

If (I .NE. J) Then !Se evita pe ea insasi

CurrentParticle = Molecula(j)

RXIJ = Arg%X - CurrentParticle%X

RYIJ = Arg%Y - CurrentParticle%Y

RZIJ = Arg%Z - CurrentParticle%Z

!Consideratii periodice

RXIJ = RXIJ - ANINT ( RXIJ / Box(Arg%Cutie) ) * Box(Arg%Cutie)

RYIJ = RYIJ - ANINT ( RYIJ / Box(Arg%Cutie) ) * Box(Arg%Cutie)

RZIJ = RZIJ - ANINT ( RZIJ / Box(Arg%Cutie) ) * Box(Arg%Cutie)

RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ !Patratul distantei dintre molecula I si molecula J

E = E + LJ(RIJSQ) !Energia de interactie depinde de tipul efectiv al potentialului

EndIf

End If

End Do

Energie_I = E

End Function Energie_I

Subrutina MC_Vol este subrutina ce implementeaza modificarea volumului celor doua compartimente, rutina ce prezinta anumite asemanari cu cea corespunzatoare din cazul ansamblului izoterm-izobar:

Subroutine MC_Vol()

Real, Dimension(2) :: EOld, ENew !Energiile particulelor din cele doua compartimente. Sunt matrici cu dimensiunea 2, corespunzand fiecarui compartiment

Real, Dimension(2) :: Scala !Factorii de scala pt cele doua cutii (matrici)

Real LnDV !Folosit la random walk in ln(V)

Type(Vector), Dimension(1:Npart) :: MoleculaOld !Vechile pozitii ale moleculelor

Real, Dimension(2) :: BoxOLd, VolOld !Vechile volume si laturi ale cutiilor (matrici)

Real Arg1, Arg2

Integer i, CeCutie

!Pune la pastrare vechii parametrii, in caz de REJECT

MoleculaOld = Molecula

BoxOld = Box

VolOld = Vol

!Energiile initiale in fiecare cutie

EOld(1) = EnergieTotala(1)

EOld(2) = EnergieTotala(2)

!Volumele noi

LnDV = Log(Vol(1)/Vol(2))+ (RNUNF()-0.5)*DVMax

Vol(1) = VTotal*Exp(LnDV)/(1+Exp(LnDV))

Vol(2) = VTotal - Vol(1)

If ((Vol(1) <= 0.0) .OR. (Vol(2) <= 0.0)) return !Una din cutii are volum negativ

Box = Vol**(1.0/3.0) !Noile laturi ale cutiilor de simulare

Scala = Box/BoxOLd !Factorii de scala pentru pozitii

Do i = 1, NPart !Rescaleaza pozitiile tuturor moleculelor, tinand cont de cutia in care se afla fiecare

CeCutie = Molecula(i)%Cutie

Molecula(i)%X = Molecula(i)%X * Scala(CeCutie)

Molecula(i)%Y = Molecula(i)%Y * Scala(CeCutie)

Molecula(i)%Z = Molecula(i)%Z * Scala(CeCutie)

End Do

!Energiile noi in fiecare cutie

ENew(1) = EnergieTotala(1)

ENew(2) = EnergieTotala(2)

Arg1 = - BetaMC*((ENew(1)-EOld(1)) + (NBox(1)+1) * Log(Vol(1)/VolOld(1)) / BetaMC)

Arg2 = - BetaMC*((ENew(2)-EOld(2)) + (NBox(2)+1) * Log(Vol(2)/VolOld(2)) / BetaMC)

!Criteriul Metropolis in volume

If (RNUNF() .GT. Exp(Arg1+Arg2)) Then

!In caz de REJECT se revine la parametrii initiali

Molecula = MoleculaOld

Box = BoxOld

Vol = VolOld

End If

End Subroutine MC_Vol

Functia EnergieTotala calculeaza energia sistemului dintr-o anumita cutie.

Real Function EnergieTotala(CeCutie)

!Calculeaza energia totala a sistemului din cutia CeCutie

Integer, Intent(IN) :: CeCutie !Cutia in care se face calculul energiei de interactie

Integer i,j, NPart_1

Real RXIJ, RYIJ, RZIJ, RIJSQ, E

Type(Vector) :: MoleculaI, MoleculaJ

E = 0.0

NPart_1 = NPart - 1

Do i = 1, Npart_1

If (CeCutie .EQ. Molecula(i)%Cutie) Then !Se asigura ca Molecula(i) se afla in cutia CeCutie

MoleculaI = Molecula(i)

Do j = i + 1, NPart

If (CeCutie .EQ. Molecula(j)%Cutie) Then !Se asigura ca Molecula(j) se afla in cutia CeCutie

MoleculaJ = Molecula(j)

RXIJ = MoleculaI%X - MoleculaJ%X

RYIJ = MoleculaI%Y - MoleculaJ%Y

RZIJ = MoleculaI%Z - MoleculaJ%Z

RXIJ = RXIJ - ANINT ( RXIJ / Box(CeCutie) ) * Box(CeCutie)

RYIJ = RYIJ - ANINT ( RYIJ / Box(CeCutie) ) * Box(CeCutie)

RZIJ = RZIJ - ANINT ( RZIJ / Box(CeCutie) ) * Box(CeCutie)

RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ !Patratul distantei dintre molecula I si molecula J

E = E + LJ(RIJSQ)

End If

End Do

End If

End Do

EnergieTotala = E

End Function EnergieTotala

Subrutina MC_Swap este ea care realizeaza schimbul de particule intre cutii:

Subroutine MC_Swap()

!Rutina de schimbat o particula intre cutii

Integer In, Out, i, Flag

Real ENew, EOld, Arg

Type(Vector) TestParticle

!Aleg cutia

If (RNUNF() .LT. 0.5) Then

In = 1; Out = 2

Else

In = 2; Out =1

End If

If (NBox(Out) == 0) return !Cutia de plecare este goala

!Aleg o molecula din cutia Out, care va pleca

Flag = 0

Do While(Flag .NE. Out)

i = Int(RNUNF()*NPart) + 1 !Aleg o molecula la intamplare

Flag = Molecula(i)%Cutie

End Do

EOld = Energie_I(Molecula(i),i) !Energia particulei in cutia veche

!In noua cutie, particula se insereaza la intamplare

TestParticle%X = (RNUNF() - 0.5) * Box(In)

TestParticle%Y = (RNUNF() - 0.5) * Box(In)

TestParticle%Z = (RNUNF() - 0.5) * Box(In)

TestParticle%Cutie = In !Se precizeaza cutia de inserare

ENew = Energie_Ghost(TestParticle) !Energia particulei in noua cutie

Arg = NBox(Out) * Vol(In) / ((Nbox(In) + 1) * Vol(Out)) * Exp(BetaMc*(EOld - ENew))

If (RNUNF() .LT. Arg) Then

!ACCEPTARE

Molecula(i)=TestParticle

NBox(Out) = NBox(Out) - 1

NBox(In) = NBox(In) + 1

End If

End Subroutine MC_Swap

Energia particulei de inserat in noua sa cutie se calculeaza cu functia Energie_Ghost:

Real Function Energie_Ghost(Arg)

!Calculeaza energia de interactie a moleculei ce ar trebui inserata in cutia CeCutie

Type(Vector), Intent(IN) :: Arg

Type(Vector) :: MoleculaJ

Integer CeCutie

Real RXIJ, RYIJ, RZIJ, RIJSQ, E

Integer j

E = 0.0

CeCutie = Arg%Cutie

Do j = 1, NPart

MoleculaJ = Molecula(j)

If (CeCutie .EQ. MoleculaJ%Cutie) Then !Se asigura ca lucrez in aceeasi cutie

RXIJ = Arg%X - MoleculaJ%X

RYIJ = Arg%Y - MoleculaJ%Y

RZIJ = Arg%Z - MoleculaJ%Z

RXIJ = RXIJ - ANINT ( RXIJ / Box(CeCutie) ) * Box(CeCutie)

RYIJ = RYIJ - ANINT ( RYIJ / Box(CeCutie) ) * Box(CeCutie)

RZIJ = RZIJ - ANINT ( RZIJ / Box(CeCutie) ) * Box(CeCutie)

RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ !Patratul distantei dintre molecula Ghost si molecula J

E = E + LJ(RIJSQ)

End If

End Do

Energie_Ghost = E

End Function Energie_Ghost