Tesztrendszerek és elvégzett számítások

A javasolt módszert négy különböző rendszeren teszteltem DMC és MD eredmények összehasonlításával. A számítások hossza 10⁸ MC lépés volt DMC, illetve 10⁶ időlépés MD szimulációk esetén. A DMC szimulációkat standard NVT sokaságon végeztem. A MD számításokban a rendszer részecskeszáma és térfogata állandó volt, a hőmérsékletet pedig Berendsen-termosztát [136] kontrollálta, amelynek csatolási időtényezője 0,1 – 0,3 ps volt. Az MD időlépések nagysága 1,0, 2,0 illetve 3,0 fs volt, a részecskék mozgásegyenleteinek megoldásra a Verlet-algoritmust [5] használtam.

A DMC szimulációkban az sα, sβ, sγ, …, stb. átlagokat a rendszer konfigurációját minden 1000. lépésben befagyasztva, és egyetlen befagyasztott konfigurációnál a véletlenszerű keresésen alapuló algoritmust egyszer végrehajtva gyűjtöttem. Az átlagokhoz szükséges d_ab konstansok direkt módon meghatározhatók a szimulációkban, de értékei megkaphatók úgy is, mint a gab(r) párkorrelációs függvények első zérustól különböző értékeihez tartozó r távolságok. Megjegyzendő, hogy a véletlenszerű bejáráson alapuló

technikában használt részecskeelléptetéshez és -elforgatáshoz is szükséges egy maximális ellépéshossz illetve egy maximális elforgatási szög definiálása. Megállapítottam azonban, hogy a (h_ij −d_ab) különbség kapott értékei független lesznek e paraméterek megválasztásától, ha garantáljuk, hogy a kereséskor az esetleges szomszédokon nem lépünk túl, vagyis a maximális ellépéshosszt nem állítjuk túlságosan nagyra. Megfelelő választás, ha az elléptetés lehetséges maximumának a részecskék előforduló legkisebb méretparaméterének felét vesszük, illetve az elforgatás maximális szögének az 1-10º értéket adjuk meg. Természetesen az ajánlott értéknél kisebbeket is választhatunk, de ekkor híg rendszerekre a keresési algoritmus futási ideje szükségtelenül hosszú lehet. Az általam végzett számításokban dab, illetve sα, sβ, sγ, …, stb. átlagok értékeit és ezekből rmax értékét a szimuláció kezdeti, az egyensúlyi vagy a stacionárius állapotot létrehozó periódusában határoztam meg. A tapasztalatok azt mutatták, hogy az s_α, s_β, s_γ, …, stb. átlagok értékei néhány ezer befagyasztott konfiguráció mintavételezése után már gyakorlatilag nem változtak.

A számításokban használt potenciállmodellek mindegyike LJ kölcsönhatási centrumokat és parciális töltéseket tartalmazó effektív párpotenciál (ld. függelék). A vegyes kölcsönhatások kiszámítására a Lorentz-Berthelot-szabályt alkalmaztam. A Coulomb-kölcsönhatások kezeléséhez a Wolf-eljárást használtam [14] αw = 2/rc

konvergenciaparaméterrel (rc = Lmin/2, ahol Lmin a szimulációban előforduló legkisebb dobozhossz). Az első három tesztrendszer termodinamikai paramétereit és összetételét az 4.

táblázat foglalja össze. A módszer érvényességét széles sűrűségtartományban teszteltem, különböző összetételek mellett. A rendelkezésre álló és a részecskék által elfoglalt térfogatok arányának kifejezésére bevezettem a következő mennyiséget:

3

ahol, ni az i-edik részecske kölcsönhatási centrumainak száma, L pedig a szimulációs doboz oldalhossza. A vizsgált tesztrendszerekre η kb. 50-66%-a a redukált számsűrűségnek,

3 3 /L

Nσ_m -nak, ahol σm a rendszer egy tipikus LJ méretparamétere. Az 4. táblázatban szereplő redukált hőmérsékletet az alábbi képlet alapján definiálhatjuk:

2

4 0

T e k

T_{∗ = B} πε ε^rσ^m _{, (59)}

ahol, e az elektron töltése, ε0 a vákuum permittivitása és εr a relatív permittivitás.

4. táblázat A számításokban használt tesztrendszerek (AX, PM, WM) paraméterei és összetétele. A PM1-PM8 rendszerek esetén C1, C2 és A jelöli a (+1)-es, (+2)-es és (-1)-es töltésű ionokat.

Rendszer Komponens N Dobozhossz η Hőmérséklet Ar 250

Az első tesztrendszer (AX) az argon és xenon ekvimoláris elegye, a második tesztrendszercsoport (PM1-PM8) pedig különböző sűrűségű és összetételű, folyékony elektrolitok elegye (a potenciálparaméterekhez ld.: függelék). Az utóbbinál valamelyest módosítottam az elektrolitok primitív modelljét [137] úgy, hogy LJ potenciált használtam az eredeti merevgömbi potenciál helyett, részben azért, hogy a javasolt módszert egy komplexebb rendszeren lehessen tesztelni (merevgömbi potenciálmodell esetén d_ab meghatározása jelentősen leegyszerűsödik). A primitív modell merevgömbi átmérőit megtartottam LJ méretparaméterekként, megengedve, hogy a kölcsönhatási centrumok esetlegesen közelebb kerülhessenek, mint ahogyan azt a merevgömbi átmérők engednék, továbbá energiaparaméterként az ε/kB = 200 K-es értéket választottam. Ilyen esetben az oldószert, ha jelen van a vizsgált elektrolitokban, implicit módon, egy ε_r > 1 relatív permittivitással vesszük figyelembe. A PM5-PM8 tesztrendszereknél az εr = 78,65 értéket használtam, és a modellparamétereket úgy választottam meg, hogy a rendszerek híg NaCl-CaCl2 vizes oldatát reprezentálják, amelyben Na⁺(C1)-, Ca²⁺(C2)- és Cl^-(A)-ionok találhatóak [138-^139140141142]. A módszert teszteltem továbbá realisztikus „site-site” kölcsönhatási potenciálokkal modellezett, szobahőmérsékletű és atmoszférikus nyomású víz/metanol elegyre is (WM1-WM2 rendszerek). A vízre illetve a metanolra a korábbi számításokban is alkalmazott potenciálmodelleket ([82] és [83]) használtam.

A módszert, az 1.6. fejezetben bemutatott, membrántranszportot megvalósító számításban is teszteltem, amelyben Lísal és munkatárainak néhány DCV-GCMD szimulációs eredményét reprodukáltuk úgy, hogy a részecsketranszlációs lépésekhez MD helyett DMC módszert használtunk (DCV-GCDMC). A két módszer felcserélésekor feltételezzük, hogy a vizsgált stacionárius rendszer teljesíti az ergodicitási kritériumot (rendszerint minden realisztikus kölcsönhatásokkal leírt stacionárius rendszerről feltételezik, hogy ergodikus [143]). Mivel a DCV-GCMD technikában a transzporttulajdonságokat csak a dinamikus cellában számítjuk, az r_max szimulációs paraméter meghatározása csak a dinamikus cella figyelembevételével történt. A kontrollcellák sűrűsége nagyságrendekkel kisebb a dinamikus cella sűrűségénél, ami azt eredményezi, hogy a legközelebbi szomszédokig megtett s utak, és ennek megfelelően a belőlük számított rmax érték jóval nagyobbak lennének a kontrollcellákban. Ha ezt az rmax értéket használnánk a kontrollcellákban, akkor irreálisan nagy transzlációs lépéseket engedélyeznénk a kontrollcellák és a dinamikus cella

találkozásánál. A számításokban ezért a kontrollcellákban is a dinamikus cellára jellemző rmax értéket használtam (a kontrollcellákban a részecskepozíciók pillanatnyi értékeit elsősorban az eltávolítási/beillesztési lépések döntik el).

A memrántranszport jellemzésére szolgáló tulajdonságok közül a komponensek relatív permeabilitására vonatkozó eredményeket vizsgáltam. Tetszőleges α komponens permeabilitása:

lm

p K j

α/

α = ∆ α , (60)

ahol lm a membrán vastagsága, ∆p_α =y_α,cc1p_cc1−y_α,cc2p_cc2, y és p a gázfázisbeli móltört illetve nyomás, cc1 és cc2 pedig a kontrolcellák indexe (p_cc1 > p_cc2). A kifejezésben szereplő ja jelöli az α komponensnek a membrán egységnyi felületére vonatkoztatott áramsűrűségét, amely meghatározható a részecskék egy kiválasztott referenciasíkon való áthaladásának regisztrálásával. Az áramsűrűség ekkor megfelel a referenciasíkon balról-jobbra illetve jobbról-balra áthaladó részecskék számának adott időegységre vonatkoztatott különbségével. Számos gyakorlati szempontból fontos dinamikus jellemző (pl.: egyensúlyi oldatokban a komponensek relatív mozgékonysága) számításakor az eltelt idő mértékegysége irreleváns, ha a keresett tulajdonság ugyan időfüggő összetevőkből számítható, de önmagában időfüggetlen. Jelen esetben ilyen tulajdonság a permeabilitások hányadosaként értelmezett dinamikus szelektivitás is.

A tesztszámításokban Lísal és munkatárainak ekvimoláris elegyekre vonatkozó eredményeit határoztam meg. A számításokban az alkalmazott Powles-membrán [144]

pórusátmérője 0,4104 nm, a membrán számsűrűsége 16,635 nm^-3, továbbá a hőmérséklet 1100 K és a cc1 kontrollcellában alkalmazott állandó nyomás 50 bar volt (a cc2 kontrollcella nyomása szimulációnként más-más értékekre volt beállítva). A komponensek illetve a membrán potenciálmodelljének paramétereit a függelék tartalmazza. Az általam végzett számításokban a cellák geometriai elrendezése (2. ábra) és méretei megfeleltek a Lísal és munkatársai által alkalmazottaknak [33]. A teljes szimulációs doboz az áramlással párhuzamos z-irányban 38,1 nm, y- és x-irányban pedig 9,525 nm széles volt. A szimulációs dobozt z irányban lágy taszító fal határolta, és a „minimum image” konvenciót a teljes szimulációs dobozra vettem figyelembe. Az átmeneti tartomány kiterjedése a kontroll- és a dinamikus cella között zérus volt, továbbá a kontrollcellák és a falak között

néhány részecske szélességnyi (0,762 nm) átmeneti tartomány lett definiálva. A 1,818 nm szélességű membrán z-irányban a szimulációs doboz középen helyezkedett el.

5.1.2. Eredmények

A 5. táblázat foglalja össze a MD referenciaszámítások és a javasolt módszert használó DMC szimulációk eredményeit. A számadatok összehasonlításából látható, hogy a bemutatott módszerrel kapott MSD-arányok jó egyezést mutatnak a referenciaértékekkel gyakorlatilag az összes tesztrendszer esetén (az MD számításokban a diffúziós állandók meghatározásának hibája rendszerint 5-15%, és emiatt a táblázatban szereplő hányadosok statisztikus bizonytalansága helyenként még ennél is nagyobb). A véletlenszerű bejáráson alapuló módszerrel meghatározott rmax értékek változása jól mutatja, hogy a rendszer sűrűségének csökkenésével nő a részecskék átlagos távolsága. A módszer érzékenységét az is jól mutatja, hogy az azonos sűrűségű PM rendszerek esetén ott kisebb rmax értéke - ennek megfelelően a legközelebbi szomszédokig ütközés nélkül megtehető úthossz -, ahol a kétszeres pozitív töltésű ionok vannak többségben: a kölcsönhatások intenzitása ekkor lokálisan kisebb sűrűséget eredményezhet (ld. pl. PM5 és PM6). Ez tehát nem mond ellent annak az előzetes feltevésnek, miszerint e tesztrendszereknél rmax számított értékének változása inkább a sűrűségre és nem annyira az összetétel változására érzékeny.

A rendszerek sűrűségének csökkenésével, azaz a szabad térfogat növekedésével a d_ab konstansokkal való korrekció, illetve a molekulák alakja egyre kisebb mértékben befolyásolja az s átlagok számított értékeit, egészen híg rendszerekre pedig a hatásuk gyakorlatilag elhanyagolható. Ezért a több kölcsönhatási centrummal rendelkező molekuláris modellekkel csak folyadékfázisban végeztem tesztszimulációkat (víz/metanol elegyek különböző összetétellel, ld. 4. táblázat). Megjegyzendő továbbá, hogy híg rendszerek esetén a kapott r_max értékek használata teljesen nem zárhatja ki a szomszédokon való „átugrások” eseményeit, bár kis értéken tartja számukat. Ugyanakkor a sűrűbb rendszerekre az ilyen események valószínűsége - a kapott kis rmax értékek miatt - gyakorlatilag zérus.

5. táblázat A javasolt és a saját szerkezeti függvény számításán alapuló módszerrel kapott dinamikus tulajdonságokból képzett hányadosok összehasonlítása a megfelelő MD számításokban

kapott értékekkel, valamint az rmax értékére kapott számadatok. A zárójelben feltüntetett számok jelentik az arányok utolsó számjegyekben jelentkező statisztikus bizonytalanságát.

MSDα / MSDß rmax/nm q/nm^-1

1javasolt módszer

2a saját szerkezeti függvény számításán alapuló módszer a tömeggel való korrekció nélkül

3a saját szerkezeti függvény számításán alapuló módszer a tömeggel való korrekcióval

A kapott eredményeket összevetettem a Berthier és Kob által javasolt eljárás [42]

alapján meghatározott MSD-arányokkal is. Minden egyes tesztrendszer esetén az rmax vs. τ görbe minimumát a 0 < rmax < Lmin/2 tartományban kerestem (τ jelöli a relaxációs időt: ld.

1.7.2. fejezet). A táblázatban a parciális szerkezeti függvény első diffrakciós csúcsának pozícióját is megadtam (q, az α vagy ß főkomponensre). Az argon/xenon ekvimoláris elegy esetén a minimális relaxációs időhöz tartozó rmax paramétert külön-külön mindkét

komponensre meghatároztam, noha a kapott q_Ar = 16,7 nm^-1 illetve a q_Xe = 17,4 nm^-1 értékekhez tartozó rmax paraméterek szimulációs hibán belül megegyeztek. A saját szerkezeti függvény számításán alapuló módszer esetén a kapott arányok a komponensek tömegének négyzetgyökével való korrekcióval és nélküle is fel vannak tüntetve (az eredeti módszer nem tartalmazza ezt a korrekciót). Jól látható, hogy a legtöbb tesztrendszernél a saját szerkezeti függvény meghatározásán alapuló módszer szolgáltatta eredmények MD referenciaszámításokkal való egyezése néha az utólagos korrekció nélkül, néha annak alkalmazásával lesz kielégítő (ld. pl. WM1 és PM3). Megjegyzendő, hogy ha a parciális szerkezeti függvény első diffrakciós csúcsa helyett más csúcsok pozíciója alapján választjuk ki qα értékét, akkor ugyan egyes esetekben javítható az egyezés, de általában nem. A számítások azt mutatták továbbá, hogy a minimális relaxációs időhöz tartozó r_max értékek qα-függése a sűrűbb rendszerek esetén a legkisebb. A legjellegzetesebb példát erre a víz/metanol tesztrendszerek esetén tapasztaltam, ahol rmax értéke gyakorlatilag független volt q_α megválasztásától.

Utolsó tesztként a nyomásgradienssel kontrollált stacionárius membrántranszportra vonatkozó dinamikus tulajdonságokat hasonlítottam össze referenciaszámítások eredményeivel, ahol a komponensek permeabilitásainak arányai fejezik ki a membrán szelektivitását. Lísal és munkatársai azt találták, hogy az általuk tanulmányozott rendszerekben a kapott permeabilitások szimulációs hibán belül nem mutatnak nyomásgradienstől való függést. Hasonló viselkedés volt megfigyelhető a DMC számításokban is a vizsgált (p_cc1-p_cc2) = 5-20 bar tartományban. A 5. táblázat utolsó sorai mutatják ezen eredmények összehasonlítását, ahol megint csak elmondható, hogy a javasolt módszer jó egyezést mutat az MD referenciaszimulációkkal. Érdekességként megemlítendő, hogy ennél a tesztrendszernél a saját szerkezeti függvény számításán alapuló módszer a relaxációs időre monoton csökkenő tendenciát adott az rmax

függvényében (a membránt alkotó részecskék távolságánál, néhány tized nm-nél jóval nagyobb r_max értékkel ugyanakkor a membrántranszport dinamikájának vizsgálata megkérdőjelezhető precizitású lett volna). A komponensek részecskéinek relaxációs ideje nagymértékben függ a részecskéket körülvevő környezettől, és ahol a környezetet reprezentáló részecskék nagy többsége szabályos rácsot alkot (kötött pozíciójú

membránrészecskék), ott a relaxációs idő nagy valószínűséggel mindig monoton csökken rmax értékének növekedésével.

Sok esetben a DMC szimulációk rmax értékeinek ±20%-os változtatása a kapott optimumhoz képest nem befolyásolja jelentős mértékben a számított MSD-arányokat.

Valószínűleg ez az oka annak, hogy a saját szerkezeti függvény számításán alapuló módszerrel meghatározott MSD-arányok is helyenként jó egyezést mutatnak a referenciaértékekkel. Általánosságban azonban az nem jelenthető ki, hogy az MSD-arányok rmax értékétől való függése kicsi. A 25. ábrán a PM rendszer közepes/kis sűrűségtartományára jellemző MSD_α/MSD_ß vs. r_max görbék láthatók. Az ábráról leolvasható, hogy a számított MSD-arány igen jelentős mértékben eltérhet a referenciaadattól, ha pl. kétszer nagyobb r_max értéket használunk a javasolt helyett.

A hagyományos MC szimulációkban az rmax paramétert a megkísérelt és elfogadott szimulációs lépések arányával szokás beállítani, amelyet gyakran kb. 50%-nak választanak közepesen sűrű vagy híg rendszerek esetén, és ha elérhető, kb. 30-40%-nak sűrű rendszerek esetén. Érdemes megjegyezni, hogy a WM és PM rendszerekre a kapott r_max értékek nagyjából a 30%, 50% és 70%-os elfogadási arányt eredményezik a szimulációs lépésekre sorrendben az η > 0,5, 0,1 < η < 0,5 és η < 0,1 „sűrűségek” esetén. Ugyanakkor nehéz egyértelmű kapcsolatot kimutatni a rendszer sűrűsége és az elvárt elfogadási valószínűség között, ugyanis néhány tesztrendszerre (AX, DCV) ez a szabály nem volt alkalmazható.

25. ábra Az A/C1 és a C2/C1 MSD-arányok rmax értékétől való függése a PM1 tesztrendszer esetén (az MD számítások eredményeinek statisztikus bizonytalansága a szimbólumok méretének közelítőleg kétszerese). A háromszög helye jelöli a véletlenszerű bejáráson alapuló módszerrel

kapott rmax értéket.