Univerzita Jana Evangelisty Purkyně v Ústí nad Labem Přírodovědecká fakulta
PARALELNÍ PROGRAMOVÁNÍ S APLIKACEMI Martin Lísal
2007
Studijní opora je určena studentům, kteří jsou zběhlí v programování v jazyce FORTRAN či v jazyce C nebo C++. Nejdříve je proveden úvod do problematiky paralelního programování. Jsou vysvětleny důvody paralelizace a její efektivita. Dále je uveden přehled úloh a architektury počítačů z hlediska paralelizace a přehled prostředků pro paralelizaci numerických výpočtů. Vlastní část studijní opory se týká úvodu do jedné z nejpoužívanějších technik paralelizace, tzv. „Message Passing Interfaceÿ (MPI) knihovny. Nejdříve jsou na jednoduchých programech v jazyce FORTRAN 90 demonstrovány základní MPI příkazy. Poté následuje přehled nejpoužívanějších MPI příkazů doplněný ilustračními programy. V závěrečné části jsou uvedeny konkrétní příklady paralelizace tří počítačových metod: paralelizace metody paralelního temperingu pro pohyb částice v jednorozměrném silovém poli a při modelování polymerů, paralelizace molekulární dynamiky a paralelizace Monte Carlo metody na příkladu Lennardovy-Jonesovy tekutiny.
Recenzovali: Ing. Zdeněk Wágner, CSc. Ing. Alexandr Babič, Ph.D. c Martin Lísal, 2007
ISBN 978-80-7044-902-8
Obsah 1 ÚVOD 1.1 Co je to paralelní počítání? . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Proč potřebujeme paralelní počítače? . . . . . . . . . . . . . . . . . . . . 1.3 Problémy spojené s vývojem paralelního počítání . . . . . . . . . . . . . 1.4 Kdy se vyplatí paralelizace? . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Závisí způsob paralelizace výpočtů na architektuře paralelních počítačů? 1.6 Rozdělení paralelních úloh z hlediska jejich spolupráce během výpočtu . . 1.7 SPMD úlohy a strategie paralelizace . . . . . . . . . . . . . . . . . . . .
. . . . . . .
5 5 5 6 6 7 10 11
2 MPI 2.1 Co je to MPI? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Vytvoření prostředí pro paralelní počítání . . . . . . . . . . . . . . . . . . 2.3 Program typu „Hello World” . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Argumenty příkazů MPI SEND a MPI RECV . . . . . . . . . . . . . . . . 2.5 Více o MPI SEND a MPI RECV . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Numerická integrace lichoběžníkovou metodou . . . . . . . . . . . . 2.5.2 Paralelní program pro numerickou integraci lichoběžníkovou metodou 2.6 Vstup a výstup v paralelních programech . . . . . . . . . . . . . . . . . . . 2.6.1 Vstup z terminálu s použitím MPI SEND a MPI RECV . . . . . . 2.6.2 Vstup z terminálu s použitím MPI BCAST . . . . . . . . . . . . . . 2.6.3 Vstup z terminálu s použitím MPI PACK a MPI UNPACK . . . . 2.6.4 Vstup ze souboru . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.5 Srovnání jednotlivých metod . . . . . . . . . . . . . . . . . . . . . . 2.7 Příkazy MPI REDUCE a MPI ALLREDUCE . . . . . . . . . . . . . . . . 2.8 Často používané MPI příkazy . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Příkazy pro vytvoření a správu paralelního prostředí . . . . . . . . 2.8.2 Příkazy pro kolektivní komunikaci . . . . . . . . . . . . . . . . . . . 2.8.3 Příkazy pro operace na proměnných distribuovaných na jednotlivých procesech . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13 13 13 14 16 17 17 18 20 21 23 24 26 27 27 30 30 31
3 APLIKACE 3.1 Paralelní „tempering” . . . . . . . . . . . . . . . . . . . . 3.1.1 Částice v jednorozměrném silovém poli . . . . . . . 3.1.2 Monte Carlo metoda . . . . . . . . . . . . . . . . . 3.1.3 Metoda paralelního „temperingu” . . . . . . . . . . 3.2 Paralelní „tempering” při modelování polymerů . . . . . . 3.2.1 Úvod . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Konformace polymeru . . . . . . . . . . . . . . . . 3.2.3 Model polymeru . . . . . . . . . . . . . . . . . . . . 3.2.4 Rosenblutův Monte Carlo algoritmus . . . . . . . . 3.2.5 Rosenblutovo Monte Carlo a paralelní „tempering”
43 43 43 44 49 57 57 58 59 60 61
3
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
39
3.3
3.4
Paralelní molekulární dynamika . . . 3.3.1 Úvod . . . . . . . . . . . . . . 3.3.2 Paralelizace . . . . . . . . . . Paralelní Monte Carlo . . . . . . . . 3.4.1 Úvod . . . . . . . . . . . . . . 3.4.2 Hybridní Monte Carlo metoda 3.4.3 Paralelizace . . . . . . . . . .
4
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
82 82 82 93 93 94 95
1
ÚVOD
1.1
Co je to paralelní počítání?
Paralelní počítání je počítání na paralelních počítačích či jinak řečeno využití více než jednoho procesoru při výpočtu (viz obr. 1).
(a) Jednoprocesorový poþítaþ
memory CPU
(b) Paralelní poþítaþ typu výpoþtový cluster
memory
memory
memory
CPU
CPU
CPU
Obrázek 1: Schema (a) jednoprocesorového a (b) paralelního počítače typu výpočtový cluster; paměť (memory), procesor (CPU).
1.2
Proč potřebujeme paralelní počítače?
• Od konce 2. světové války dochází k bouřlivému rozvoji počítačových simulací (computational science). Počítačové simulace řeší problémy, které teoretická věda (theoretical science) nedokáže vyřešit analyticky či problémy, které jsou obtížné nebo nebezpečné pro experimentální vědu (experimental science). Jako příklad uveďme proudění okolo trupu letadla, které teoretická věda dokáže popsat parciálními diferenciálními rovnicemi; ty ale nedokáže analyticky řešit. Dalším příkladem je určení vlastností plazmy. Jelikož plazma existuje při vysokých teplotách, je experimentální měření většiny jejích vlastností obtížné či přímo nemožné. Počítačové simulace vyžadují náročné numerické výpočty, které se v současnosti neobejdou bez nasazení paralelizace. 5
• Vývoj procesorů naráží (či v budoucnu narazí) na svůj materiálový limit t.j. vývoj procesorů nepůjde zvyšovat do nekonečna. (To pull a bigger wagon, it is easier to add more oxen than to grow a gigantic ox.) • Cena nejvýkonějších (advanced) procesorů na trhu obvykle roste rychleji než jejich výkon. (Large oxen are expensive.)
1.3
Problémy spojené s vývojem paralelního počítání
• Hardware: v současné době lze konstruovat paralelní počítače mající společnou (sdílenou) paměť pro maximálně asi 100 procesorů. Trendem je tedy spojení jednoprocesorových či několikaprocesorových počítačů do clusteru pomocí rychlé komunikační sítě tzv. switche. Pozornost v oblasti hardwaru se proto převážně zaměřuje na vývoj switchů, které umožňují rychlost meziprocesorové komunikace srovnatelnou s rychlostí komunikace mezi procesorem a vnitřní pamětí počítače. • Algoritmy: bohužel se ukazuje, že některé algoritmy používané na jednoprocesorových počítačích nelze dobře paralelizovat. To vede k nutnosti vytvářet zcela nové paralelní algoritmy. • Software: existují prostředky pro automatickou paralelizaci jakými jsou např. paralelní kompilátory [High Performance Fortran (HPF) Standard] či paralelní numerické knihovny. Platí (a pravděpodobně dlouho bude platit), že nejlepší paralelizace je ta, kterou si uděláme sami.
1.4
Kdy se vyplatí paralelizace?
Zrychlení výpočtu paralelizací lze odhadnout pomocí Amdahlova pravidla: 100 − P +
P n
(1)
kde P je část programu, kterou lze paralelizovat a n je počet použitých procesorů. Amdahlovo pravidlo nám říká, kolik procent výpočtového času bude trvat paralelizovaný program oproti stejnému programu spuštěnému na jednoprocesorovém počítači. Uveďme si jednoduchý ilustrativní příklad. Představme si, že máme sečíst 1 000 000 000 prvků vektoru a. Načtení prvků trvá 8 % celkového výpočtového času, vlastní výpočet součtu trvá 90 % celkového výpočtového času a tisk trvá 2 % celkového výpočtového času. Je zřejmé, že jen výpočet součtu lze paralelizovat t.j. P = 90 % (viz. obr. 2). V případě použití n = 10 procesorů dostaneme z Amdahlova pravidla výsledek 19 % t.j. zrychlení výpočtu paralelizací přibližně 5krát. Z Amdahlova pravidla a uvedeného příkladu plynou 6
dva závěry: (i) neplatí přímá úměra mezi počtem procesorů a urychlením výpočtu („užití 10 procesorů neznamená urychlení výpočtu 10krát”) a (ii) pro limitní případ n → ∞ je urychlení výpočtu konečné a rovno 100/(100 − P ).
(a) Sériový program
(b) Paralelní program
a.out
S1
P
S2
S1
P
S1
P
S2
S1
P
S2
S1
P
S2
S2
Obrázek 2: Schéma (a) sériového a (b) paralelního programu (běžícího na 4 procesorech). S1 a S2 jsou části programu, které nelze paralelizovat a P je paralelizovatelná část programu.
1.5
Závisí způsob paralelizace výpočtů na architektuře paralelních počítačů?
Způsob paralelizace výpočtů bohužel závisí na architektuře paralelních počítačů. V současnosti se nejvíce používají dva typy architektur: (i) paralelní počítače se sdílenou pamětí (shared memory) a (ii) paralelní počítače s distribuovanou pamětí (distributed memory); viz obr. 3. 7
(a) Paralelní poþítaþ se sdílenou pamČtí (shared memory)
(b) Paralelní poþítaþ s distribuovanou pamČtí (distributed memory)
switch ~ komunikace (shared) memory
CPU
CPU
CPU
CPU
memory
memory
memory
memory
CPU
CPU
CPU
CPU
Obrázek 3: Schema paralelního počítače (a) se sdílenou pamětí (shared memory) a (b) s distribuovanou pamětí (distributed memory).
Z hlediska způsobu paralelizace je výhodnější architektura se sdílenou pamětí, neboť sdílení a výměna informací mezi procesory během výpočtu se děje přes sdílenou (společnou) paměť. Paralelizace na „shared-memory” počítačích se provádí pomocí tzv. OpenMP. OpenMP je seznam direktiv, které se vkládají na místa, které chceme paralelizovat. Jinými slovy říkáme počítači „kde a co” paralelizovat. Uveďme si jako příklad paralelizaci cyklu (do) pomocí OpenMP: !$OMP PARALLEL DO do i = 2, n b(i) = (a(i) + a(i-1)) / 2.0 end do !$OMP END PARALLEL DO
Direktiva před začátek cyklu, !$OMP PARALLEL DO, říká: „paralelizuj cyklus”, zatímco direktiva na konci cyklu, !$OMP END PARALLEL DO, říká: „ukonči paralelizaci cyklu”.1 Paralelní počítače se sdílenou pamětí jsou oproti paralelním počítačům s distribuovanou pamětí dražší a počet procesorů sdílející jednu paměť je v současnosti omezen maximálně asi na 100. Proto se běžněji setkáváme s „distributed-memory” počítači. „Distributedmemory” počítače mají fyzicky oddělené paměti a komunikace (sdílení dat) mezi procesory 1
Jelikož direktivy OpenMP začínají ”!” a ten se ve FORTRANu90 interpretuje jako začátek komentářového řádku, lze programy s direktivami OpenMP spouštět bez problému na jednoprocesorových počítačích.
8
zajišťuje switch. Paralelizace na „distributed-memory” počítačích je obtížnější než paralelizace na „shared-memory” počítačích, neboť musíme zajistit a pracovat s předáváním informací (message passing) mezi procesory během výpočtu. Na paralelních počítačích s distribuovanou pamětí se používají dva systémy: (i) Parallel Virtual Machine (PVM) vhodný pro heterogenní clustery a (ii) Message Passing Interface (MPI) vhodný pro homogenní clustery. PVM a MPI je knihovna příkazů pro FORTRAN, C a C++, jež umožňují předávání informací mezi procesory během výpočtu.
þas komunikace (overhead)
Komunikace představuje časovou ztrátu, o kterou se sníží zrychlení výpočtu paralelizací dané Amdahlovým pravidlem. Čas komunikace (overhead) jako funkce množství předávané informace je znázorněn na obr. 4. Z obr. 4 je patrno, že čas komunikace závisí na latenci (času, který potřebuje procesor a switch k „navázání komunikace”) a dále, že „overhead” roste s množstvím předávané informace mezi procesory. Z obr. 4 plynou následující zásady, které bychom měli dodržovat při paralelizaci na počítačích s distribuovanou pamětí: (i) přenášet co nejméně informací a (ii) komunikovat co nejméně.
latence
množství pĜenášené informace
Obrázek 4: Čas komunikace (overhead) versus množství předávané informace mezi procesory.
9
1.6
Rozdělení paralelních úloh z hlediska jejich spolupráce během výpočtu
Podle spolupráce během výpočtu můžeme rozdělit paralelní úlohy na • MPMD (Multiple Program Multiple Data) úlohy
• SPMD (Single Program Multiple Data) úlohy
MPMD úlohy lze dále rozdělit na typ
• „Master/Worker” (obr. 5a)
• „Coupled Multiple Analysis” (obr. 5b)
Úlohy typu „Master/Worker” jsou charakteristické existencí řídícího procesu (master, a.out v obr. 5a), který řídí ostatní procesy (worker, b.out v obr. 5a) a nechává si od nich počítat mezivýsledky. Úlohy typu „Coupled Multiple Analysis” jsou charakteristické použitím nezávislých procesů (a.out, b.out a c.out v obr. 5b), které řeší různé aspekty určitého problému. Např. při návrhu trupu letadla by program a.out řešil proudění, b.out by řešil pevnostní analýzu a c.out analýzu teplotní. SPMD úlohy (obr. 5c) používají stejnou kopii programu (single program) na všech procesorech. Procesy ale během výpočtu zpracovávají různá data (multiple data). SPMD úlohy jsou typické pro numerické výpočty a my se SPMD úlohami budeme výhradně zabývat v rámci tohoto předmětu. 10
(a) MPMD úloha typu "Master/Worker" a.out b.out
a.out
(b) MPMD úloha typu "Coupled Multiple Analysis" b.out c.out
(c) SPMD úloha a.out
Obrázek 5: Paralelní úlohy z hlediska jejich spolupráce během výpočtu. (a) MPMD (Multiple Program Multiple Data) úloha typu „Master/Worker”, (b) MPMD úloha typu „Coupled Multiple Analysis” a (c) SPMD (Single Program Multiple Data) úloha.
1.7
SPMD úlohy a strategie paralelizace
Uveďme si jednoduchý příklad paralelizace programu pro součet prvků vektoru, který ozřejmí pojem SPMD úlohy a naznačí strategii paralelizace SPMD úloh. Představme si, že máme vektor a mající 9 000 000 000 prvků. Úkol je sestavit paralelní program pro souP čet prvků vektoru a ( 9i=1000 000 000 ai ), přičemž máme k dispozici 3 procesory. Rozdělme si úlohu pomyslně na 3 části: (i) načtení vektoru a, S1 , (ii) vlastní výpočet součtu, P , a (iii) tisk součtu, S2 . Je zřejmé, že jen vlastní výpočet součtu P lze paralelizovat.2 Paralelizaci provedeme tak, že každý procesor spočte 1/3 součtu a tu pak pošle na jeden procesor, který provede celkový součet, jež následně vytiskne. Paralelní program pro součet prvků vektoru zapsaný v jazyce FORTRAN může vypadat následovně (v místech programu, jež by měla obsahovat MPI paralelní příkazy, jsou komentáře označeny jako „. . .”): program soucet implicit none integer :: i,my_rank,p,n_loc,i_begin,i_end real :: a(9000000000), & sum_loc,sum_tot ! 2
Načtení dat se nedoporučuje paralelizovat, neboť způsob načítání dat se liší systém od systému.
11
open(10,file=’input.dat’) do i=1,9000000000 read(10,*) a(i) end do close(10) ! "MPI: zjisti počet procesů ’p’, které se podílí na výpočtu" "MPI: zjisti moje pořadí ’my_rank’ ! n_loc=9000000000/p ! Předpokládáme i_begin=my_rank*n_loc+1 i_end=(my_rank+1)*n_loc sum_loc=0.0 do i=i_begin,i_end sum_loc=sum_loc+a(i) end do ! "MPI: zajisti sečtení jednotlivých a pošli ’sum_tot’ např. na proces ! if(my_rank == 0) then write(*,*) sum_tot end if stop "soucet: Konec vypoctu!" ! end program soucet
ve výpočtu; my_rank=0,1,...,p-1" dělitelnost beze zbytku :-)
’sum_loc’ do ’sum_tot’ mající ’my_rank=0’"
12
2
MPI
2.1
Co je to MPI?
MPI (Message Passing Interface) je knihovna funkcí a podprogramů, která umožňuje: • vytvořit prostředí pro paralelní počítání; • komunikaci mezi procesy během výpočtu např. posílání informací z procesu na proces; • kolektivní operace např. sečtení mezivýsledků z procesů a uložení celkové součtu na určitý proces; • topologizaci procesů t.j. např. seskupení určitého počtu procesů za účelem jejich spolupráce v rámci výpočtů. MPI lze použít ve spojení s programovacími jazyky FORTRAN, C a C++.
2.2
Vytvoření prostředí pro paralelní počítání
Vytvořením prostředí pro paralelní počítání rozumíme: • inicializaci a ukončení používání MPI; • definování celkového počtu procesů nprocs, které se účastní výpočtu; • definování pořadí procesu myrank. Demonstrujme si vytvoření paralelního prostředí na jednoduchém programu, který vypíše na obrazovku pořadí jednotlivých procesů.3 program par_env implicit none include ’mpif.h’ ! Budeme pouzivat MPI knihovnu ! integer :: ierr,nprocs,myrank ! call MPI_INIT(ierr) ! INITializuje MPI ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! Default COMMunikator a celkovy pocet (SIZE) procesu ’nprocs’ ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! Poradi (RANK) ’myrank’; myrank = 0, 1, ..., nprocs-1 ! print*," Poradi procesu =",myrank ! call MPI_FINALIZE(ierr) ! Ukonci (FINALIZE) MPI end program par_env 3 Pokud nebude řečeno jinak, pracujeme s SPMD modelem paralelního programu: na všech uzlech se rozběhne stejná kopie programu (SP, single program), každá kopie programu pracuje s různými daty (MD, multiple data).
13
Komentář k programu: • Include ’mpif.h’ říká, že budeme používat MPI knihovnu. • MPI INIT(ierr) inicializuje MPI; ierr je „chybová” proměnná typu INTEGER (ierr = 0 když je vše v pořádku). • MPI COMM SIZE(MPI COMM WORLD,nprocs,ierr) definuje celkový počet procesů nprocs, které se účastní výpočtu. • MPI COMM RANK(MPI COMM WORLD,myrank,ierr) definuje pořadí procesu myrank. • MPI FINALIZE(ierr) ukončí MPI. • V celém programu používáme předdefinovaný (default) tzv. „komunikátor” MPI COMM WORLD – označení pro skupinu spolupracujících procesů.
2.3
Program typu „Hello World”
Při učení nového programovacího jazyka se pravděpodobně každý setkal s programem, který vypíše na obrazovku „pozdrav” - program typu „Hello World”. Napišme si takovou verzi programu pomocí MPI. MPI „Hello World” program pošle z procesů majících pořadí myrank = 1, 2, . . . , nprocs −1 informaci o svém pořadí na proces s pořadím 0 a proces s pořadím 0 vytiskne tuto informaci na obrazovku. program hello_world implicit none include ’mpif.h’ ! Budeme pouzivat MPI knihovnu ! integer :: STATUS(MPI_STATUS_SIZE) ! Pomocna promenna pouzita v MPI_RECV integer :: ierr,nprocs,myrank,tag,isource,act_proc ! call MPI_INIT(ierr) ! INITializujeme MPI ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! Default COMMunikator a celkový pocet (SIZE) procesu ’nprocs’ ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! Poradi (RANK) ’myrank’; myrank = 0, 1, ..., nprocs-1 ! tag=1 if(myrank /= 0) then ! Ne_Mastr procesy act_proc=myrank call MPI_SEND(act_proc,1,MPI_INTEGER,0,tag,MPI_COMM_WORLD,ierr) ! Ne_Mastr procesy SEND info ’act_proc’ else ! Mastr proces do isource=1,nprocs-1 call MPI_RECV(act_proc,1,MPI_INTEGER,isource,tag,MPI_COMM_WORLD, & STATUS,ierr) ! Mastr proces RECV info ’act_proc’
14
print*," Hello World =",act_proc enddo endif ! call MPI_FINALIZE(ierr) ! Ukonci (FINALIZE) MPI end program hello_world
Komentář k programu: • MPI SEND(. . . ) posílá informace. • MPI RECV(. . . ) zajistí obdržení informace. • 1. argument v MPI SEND/MPI RECV je posílaná/obdržená informace. • 2. argument v MPI SEND/MPI RECV je velikost posílané/obdržené informace (1 pro skalár, počet prvků vektoru či matice, atd.). • 3. argument v MPI SEND/MPI RECV značí MPI typ, který přibližně odpovídá typům v jazycích FORTRAN, C a C++. • 4. argument v MPI SEND značí, na jaký proces se má zpráva poslat (v našem případě na proces 0). • 4. argument v MPI RECV značí, z jakých procesů zpráva přichází (v našem případě z procesů 1, 2, . . . , nprocs −1). • 5. argument v MPI SEND/MPI RECV je proměnná typu INTEGER (0 – 32767), která umožňuje dodatečné (pomocné) označení zprávy. • Proměnná STATUS v MPI RECV zpětně „informuje” MPI SEND o úspěšném doručení zprávy.4 • MPI RECV je až na STATUS zrcadlovým obrazem MPI SEND. • MPI SEND a MPI RECV jsou vzájemně synchronizované příkazy t.j. běh programu pokračuje, až když se oba příkazy dokončí.
4
Ostatní argumenty v MPI SEND/MPI RECV byly již vysvětleny.
15
2.4
Argumenty příkazů MPI SEND a MPI RECV
V programu typu „Hello World” jsme použili příkaz MPI SEND pro posílání informace z uzlu na uzel a příkaz MPI RECV pro přijímání informací z uzlů. Popišme si detailně argumenty těchto příkazů. Tyto dva příkazy lze obecně zapsat jako: call MPI_SEND(send_message,count,MPI_data_type,dest,tag,comm, & ierr)
call MPI_RECV(recv_message,count,MPI_data_type,source,tag,comm, & status,ierr) kde • send message je posílaná informace, jejíž typ specifikuje MPI data type a jež může být skalár, vektor či matice. • recv message je přijímaná informace, jejíž typ specifikuje MPI data type a jež může být skalár, vektor či matice. • count je typu INTEGER a vyjadřuje velikost posílané resp. přijímané informace; např. 1 pro skalár či 10 pro matici (2,5). • MPI data type: MPI data type MPI INTEGER MPI REAL MPI DOUBLE PRECISION MPI COMPLEX MPI LOGICAL MPI CHARACTER MPI PACKED MPI BYTE
FORTRAN INTEGER REAL DOUBLE PRECISION COMPLEX LOGICAL CHARACTER
• dest je typu INTEGER a vyjadřuje pořadí procesu, na který se send message posílá. • source je typu INTEGER a vyjadřuje pořadí procesu, z kterého se recv message přijímá. • tag je proměnná typu INTEGER, která může nabývat hodnot od 0 do 32767 a která umožňuje dodatečné (pomocné) označení posílané resp. přijímané informace. • comm je označení pro skupinu spolupracujících procesů („komunikátor”); MPI COMM WORLD je předdefinovaný (default) „komunikátor”. 16
• status je vektor typu INTEGER deklarovaný jako integer :: status(MPI_STATUS_SIZE)
status zpětně „informuje” MPI SEND o úspěšném doručení zprávy. • ierr je „chybová” proměnná typu INTEGER; ierr = 0 když je vše v pořádku.
2.5
Více o MPI SEND a MPI RECV
Ukažme si další použití MPI SEND a MPI RECV na paralelní verzi programu pro numerický výpočet určitého integrálu lichoběžníkovou metodou. 2.5.1
Numerická integrace lichoběžníkovou metodou
Úkolem je spočíst integrál funkce f (x): I=
Z b
f (x) dx
(2)
a
Krok 1: rozdělme interval ha, bi na n stejných dílů: i = 0, 1, 2, . . . , n − 1, n
xi = a + ih x0 = a xn = b
(3) (4) (5)
kde
b−a n Krok 2: aproximujme funkci f (x) mezi body xi a xi+1 přímkou. Krok 3: body xi , xi+1 , f (xi ) a f (xi+1 ) tvoří lichoběžník, jehož obsah je h=
(6)
f (xi ) + f (xi+1 ) 2
(7)
Ai = h
Krok 4: sečtením obsahů všech lichoběžníků, dostaneme přibližnou hodnotu integrálu (2): X f (a) + f (b) n−1 I ≈ h + f (xi ) 2 i=1 "
"
≈
#
f (a) f (b) + f (a + h) + f (a + 2h) + . . . + f (b − h) + 2 2 17
(8) #
(9)
Sériový program pro numerickou integraci lichoběžníkovou metodou může vypadat následovně: program ser_lich_int implicit none ! integer :: i,n real :: f,a,b,h,x,int_sum,Int ! print *,"Dolni mez intervalu a:" read *,a print *,"Horni mez intervalu b > a:" read *,b print *,"Zadejte deleni intervalu n > 0:" read *,n ! h=(b-a)/REAL(n) ! "krok" int_sum=(f(a)+f(b))/2.0 ! [f(a)+f(b)]/2 do i=1,n-1 x=a+REAL(i)*h ! x=a+i*h int_sum=int_sum+f(x) ! pomocna suma end do Int=h*int_sum ! vysledny integral ! print *,"Integral =",Int end program ser_lich_int ! function f(x) implicit none ! real :: f,x f=x*x end function f
2.5.2
Paralelní program pro numerickou integraci lichoběžníkovou metodou
Paralelní verze programu pro numerickou integraci lichoběžníkovou metodou 1. rozdělí interval ha, bi na podintervaly hai , bi i dle celkového počtu procesů nprocs, 2. pro jednotlivé podintervaly spočtou jednotlivé procesy integrály Z bi
f (x) dx
ai
pomocí lichoběžníkové metody 3. a tyto integrály se pošlou na proces mající pořadí 0, kde se sečtou a proces mající pořadí 0 součet (celkový integrál) vytiskne. Paralelní program pro numerickou integraci lichoběžníkovou metodou může vypadat následovně:5 6 5 Abychom se vyhnuli problematice načítání vstupních dat v paralelních programech, zadejme vstupní hodnoty a, b a n přímo do programu. 6 Dále předpokládáme, že n je beze zbytku dělitelné nprocs.
18
! ! Trapezoidal Rule, MPI version 1. ! ! a, b, and n are in DATA statement ! ! Algorithm: ! 1. Each process calculates "its" interval of integration. ! 2. Each process estimates the integral of f(x) over its interval ! using the trapezoidal rule. ! 3a. Each process /= 0 sends its integral to 0. ! 3b. Process 0 sums the calculations received from the individual ! processes and prints the result. ! program par_lich_int1 implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank source, & ! process sending integer dest, & ! process receiving integer tag,ierr integer :: status(MPI_STATUS_SIZE) integer :: n, & ! # of trapezoids local_n ! # of trapezoids for a processor ! real :: a, & ! left endpoint b, & ! right endpoint h, & ! trapezoid base length local_a, & ! left endpoint for a processor local_b, & ! right endpoint for a processor integral, & ! integral over a processor total ! total integral ! data a,b,n /0.0,1.0,1024/ ! data for integration, n must be even data dest,tag /0,50/ ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! if(mod(n,nprocs) /= 0) stop "par_lich_int1: Wrong n/nprocs ratio!" h=(b-a)/REAL(n) ! h is the same for all processes local_n=n/nprocs ! # of trapezoids ! ! length of interval of integration for a process ! local_a=a+REAL(myrank)*REAL(local_n)*h local_b=local_a+REAL(local_n)*h ! ! calculate integral on a process ! call Trap(local_a,local_b,local_n,h,integral) ! ! add up the integrals calculated by each process ! if(myrank /= 0) then call MPI_SEND(integral,1,MPI_REAL,dest,tag, & MPI_COMM_WORLD,ierr) else total=integral do source=1,nprocs-1 call MPI_RECV(integral,1,MPI_REAL,source,tag, & MPI_COMM_WORLD,status,ierr) total=total+integral enddo endif ! ! print the results ! if(myrank == 0) then
19
write(*,’(1x,a,e13.5)’) " Integral =",total endif ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program par_lich_int1 ! ! Subroutine trapezoid for a processor ! subroutine Trap(local_a,local_b,local_n,h,integral) implicit none ! integer :: local_n,i real :: local_a,local_b,h,integral,x,f ! integral=(f(local_a)+f(local_b))/2.0 x=local_a do i=1,local_n-1 x=x+h integral=integral+f(x) enddo integral=integral*h end subroutine Trap ! ! Function for integration ! function f(x) implicit none ! real :: f,x f=x*x end function f
2.6
Vstup a výstup v paralelních programech
MPI standard neobsahuje speciální příkazy pro načítání dat či jejich výstup na terminál nebo do souboru. Vstup a výstup v paralelních programech se proto musí řešit pomocí příkazů programovacího jazyka read, write či print a MPI příkazů pro komunikaci mezi uzly. Je nutno si uvědomit, že např. způsob načtení dat použitý v sériovém programu pro numerickou integraci lichoběžníkovou metodou: print *,"Dolni mez intervalu a:" read *,a print *,"Horni mez intervalu b > a:" read *,b print *,"Zadejte deleni intervalu n > 0:" read *,n
může v paralelním programu fungovat různě v závislosti na systému a architektuře počítače. Při podobném použití příkazu read v paralelních programech může dojít k následujícím situacím: • Pokud jen jeden proces „umí” číst z terminálu, data se načtou jen na tento proces a ne na ostatní procesy. 20
• Pokud „umí” všechny procesy číst, pak může nastat případ zadávat data tolikrát, kolik procesorů používáme.7 Situace se dále komplikuje, pokud načítáme data ze souboru a např. více procesů musí číst současně data z jednoho souboru. Obdobné problémy a nejednoznačnosti jsou spojeny s výstupem na terminál či do souborů. Systémy a architektury počítačů vždy minimálně umožňují následující: • Alespoň jeden proces dokáže číst z terminálu a alespoň jeden proces dokáže zapisovat na terminál. • Každý proces dokáže číst ze souboru či do souboru zapisovat, pokud ostatní procesy z tohoto souboru nečtou či do něj nezapisují. • Procesy dokáží číst současně z různých souborů či do různých souborů zapisovat. Na základě těchto „minimálních” schopností systémů a architektur počítačů se doporučuje používat následující dvě zásady při řešení vstupu a výstupu do paralelních programů přes terminál: 1. Načíst data jen jedním procesem (např. procesem s myrank = 0) a následně rozeslat data z tohoto procesu na ostatní procesy. 2. Poslání výsledků výpočtů z jednotlivých procesů na jeden určitý proces (např. proces s myrank = 0) a tisk výsledků z tohoto procesu na terminál. V následujícím si ukažme různá řešení vstupů dat do paralelních programů z terminálu a souboru na příkladu paralelního programu pro numerickou integraci lichoběžníkovou metodou. Abychom se vyhnuli opakovanému psaní paralelního programu pro numerickou integraci lichoběžníkovou metodou, nahradíme příkaz data a,b,n /0.0,1.0,1024/ ! data for integration, n must be even
podprogramem, který provede načtení dat. Zbytek programu se nezmění.
2.6.1
Vstup z terminálu s použitím MPI SEND a MPI RECV
V první verzi vstupu z terminálu načte proces s pořadím myrank = 0 data a tato data pošle s použitím příkazů MPI SEND a MPI RECV na ostatní procesy. Příslušný podprogram může vypadat následovně: ! ! Algorithm: ! 1. Process 0 reads data. 7
Dále není zaručeno, že čtení proběhne podle pořadí procesů.
21
! 2. Process 0 sends data to other processes. ! subroutine Get_Data1(a,b,n,myrank,nprocs) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank nprocs, & ! # of processes source, & ! process sending integer dest, & ! process receiving integer tag,ierr integer :: n ! # of trapezoids integer :: status(MPI_STATUS_SIZE) ! real :: a, & ! left endpoint b ! right endpoint ! if(myrank == 0) then print *," a:" read*,a print *," b:" read*,b print *," n:" read*,n ! do dest=1,nprocs-1 tag=0 call MPI_SEND(a,1,MPI_REAL,dest,tag,MPI_COMM_WORLD,ierr) tag=1 call MPI_SEND(b,1,MPI_REAL,dest,tag,MPI_COMM_WORLD,ierr) tag=2 call MPI_SEND(n,1,MPI_INTEGER,dest,tag,MPI_COMM_WORLD,ierr) enddo else tag=0 call MPI_RECV(a,1,MPI_REAL,0,tag,MPI_COMM_WORLD, & status,ierr) tag=1 call MPI_RECV(b,1,MPI_REAL,0,tag,MPI_COMM_WORLD, & status,ierr) tag=2 call MPI_RECV(n,1,MPI_INTEGER,0,tag,MPI_COMM_WORLD, & status,ierr) endif end subroutine Get_Data1
Tato verze programu volá několikrát příkazy MPI SEND a MPI RECV. Někdy proto bývá výhodnější načtená data nejdříve uložit do pomocného vektoru a ten poslat s použitím příkazů MPI SEND a MPI RECV na ostatní procesy. V takovém případě se sníží latentní čas komunikace. Příslušný podprogram může vypadat následovně: ! ! Algorithm: ! 1. Process 0 reads data and save them to a vector. ! 2. Process 0 sends data as vector to other processes. ! subroutine Get_Data11(a,b,n,myrank,nprocs) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank
22
nprocs, & ! # of processes source, & ! process sending integer dest, & ! process receiving integer tag,ierr integer :: n ! # of trapezoids integer :: status(MPI_STATUS_SIZE) ! real :: a, & ! left endpoint b, & ! right endpoint temp_vec(3) ! temporary vector ! if(myrank == 0) then print *," a:" read*,a print *," b:" read*,b print *," n:" read*,n ! temp_vec=(/a,b,REAL(n)/) ! do dest=1,nprocs-1 tag=0 call MPI_SEND(temp_vec,3,MPI_REAL,dest,tag,MPI_COMM_WORLD, & ierr) enddo else tag=0 call MPI_RECV(temp_vec,3,MPI_REAL,0,tag,MPI_COMM_WORLD,status, ierr) a=temp_vec(1) b=temp_vec(2) n=INT(temp_vec(3)) endif end subroutine Get_Data11
2.6.2
&
Vstup z terminálu s použitím MPI BCAST
V druhé verzi vstupu z terminálu načte opět proces s pořadím myrank = 0 data a tato data se rozešlou s použitím příkazu MPI BCAST8 na ostatní procesy. Příslušný podprogram může vypadat následovně: ! ! Algorithm: ! 1. Process 0 reads data. ! 2. Process 0 broadcasts data to other processes. ! subroutine Get_Data2(a,b,n,myrank) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank ierr integer :: n ! # of trapezoids ! real :: a, & ! left endpoint b ! right endpoint 8
Příkaz MPI BCAST si můžeme představit jako „sloučení” příkazů MPI SEND a MPI RECV do jednoho.
23
! if(myrank == 0) then print *," a:" read*,a print *," b:" read*,b print *," n:" read*,n endif call MPI_BCAST(a,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(b,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) end subroutine Get_Data2
Příkaz MPI BCAST lze obecně zapsat jako:
call MPI_BCAST(broadcast_message,count,MPI_data_type,root,comm,ierr)
kde význam všech argumentů kromě root je zřejmý z předchozích probíraných MPI příkazů. root je typu INTEGER a vyjadřuje pořadí procesu, z kterého se broadcast message rozesílá.
2.6.3
Vstup z terminálu s použitím MPI PACK a MPI UNPACK
V třetí verzi vstupu z terminálu načte opět proces s pořadím myrank = 0 data. Tato data před rozesláním „zabalí” do pomocné proměnné buffer s použitím příkazu MPI PACK. Proměnná buffer se pak rozešle na ostatní procesy s pomocí příkazu MPI BCAST, kde se zpětně „rozbalí” s použitím příkazu MPI UNPACK.9 Příslušný podprogram může vypadat následovně: ! ! Use of Pack/Unpack ! subroutine Get_Data3(a,b,n,myrank) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank position, & ! in the buffer ierr integer :: n ! # of trapezoids ! real :: a, & ! left endpoint b ! right endpoint ! character(len=100) :: buffer ! if(myrank == 0) then print *," a:" read*,a print *," b:" 9
Užití MPI PACK a MPI UNPACK je výhodné při posílání velkého množství dat různých typů.
24
read*,b print *," n:" read*,n ! ! Pack the data into buffer. ! Beginning of buffer ’position=0’ ! position=0 ! Position is in/out call MPI_PACK(a,1,MPI_REAL,buffer,100,position, & MPI_COMM_WORLD,ierr) ! Position has been incremented; ! the first free location in buffer call MPI_PACK(b,1,MPI_REAL,buffer,100,position, & MPI_COMM_WORLD,ierr) ! Position has been incremented again call MPI_PACK(n,1,MPI_INTEGER,buffer,100,position, & MPI_COMM_WORLD,ierr) ! Position has been incremented again ! ! Broadcast contents of buffer ! call MPI_BCAST(buffer,100,MPI_PACKED,0,MPI_COMM_WORLD, & ierr) else call MPI_BCAST(buffer,100,MPI_PACKED,0,MPI_COMM_WORLD, & ierr) ! ! Unpack the data from buffer ! position=0 call MPI_UNPACK(buffer,100,position,a,1,MPI_REAL, & MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,100,position,b,1,MPI_REAL, & MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,100,position,n,1,MPI_INTEGER, & MPI_COMM_WORLD,ierr) endif end subroutine Get_Data3
Příkazy MPI PACK a MPI UNPACK lze obecně zapsat jako:
call MPI_PACK(pack_data,count,MPI_data_type,buffer,size,position, & comm,ierr)
call MPI_UNPACK(buffer,size,position,unpack_data,count,MPI_data_type, & comm,ierr)
kde význam všech argumentů kromě buffer, size a position je zřejmý z předchozích probíraných MPI příkazů. Všiměte si, že proměnné „zabalené” pomocí příkazu MPI PACK (buffer v našem příkladu) jsou v MPI příkazech typu MPI PACKED. • buffer je pomocná proměnná, jež musí být deklarovaná jako typ CHARACTER. 25
• size je typu INTEGER a vyjadřuje velikost („délku”) pomocné proměnné buffer.10 • position je typu INTEGER a vyjadřuje pozici v buffer, od které se pack data ukládají či unpack data jsou uloženy. Před prvním použitím MPI PACK či MPI UNPACK je nutno position vynulovat a pak je position automaticky aktualizována voláním MPI PACK či MPI UNPACK.
2.6.4
Vstup ze souboru
Předpokládejme, že vstupní data a, b a n jsou v souboru input.dat. Aby došlo ke správnému načtení dat jednotlivými procesy, je nutno zajistit načítání dat z tohoto souboru postupně. To lze docílit použitím příkazu MPI BARRIER. Příslušný podprogram může vypadat následovně: ! ! Reading from file ! subroutine Get_Data4(a,b,n,myrank,nprocs) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank nprocs, & ! # of processes irank,ierr integer :: n ! # of trapezoids ! real :: a, & ! left endpoint b ! right endpoint ! do irank=0,nprocs-1 if(irank == myrank) then open(10,file="input.dat",status="old") read(10,*) a read(10,*) b read(10,*) n close(10) endif call MPI_BARRIER(MPI_COMM_WORLD,ierr) enddo end subroutine Get_Data4
Příkaz MPI BARRIER funguje tak, že program pokračuje až poté, když všechny procesy zavolají tento příkaz. Obecný tvar příkazu MPI BARRIER je: call MPI_BARRIER(comm,ierr) kde význam argumentů je zřejmý z předchozích probíraných MPI příkazů. 10
size musí být větší než součet velikostí („délek”) „balených” dat.
26
2.6.5
Srovnání jednotlivých metod
Použití výše zmíněných metod závisí na množství a typech rozesílaných dat: • V případě rozesílání malého množství dat různých typů, několikanásobné použití MPI SEND a MPI RECV či MPI BCAST je nejjednodušším řešením. • V případě rozesílání většího množství dat podobných typů (např. REAL, INTEGER), uložení dat do pomocného vektoru a použití MPI SEND a MPI RECV či MPI BCAST se jeví jako nejvýhodnější řešení. • V případě rozesílání velkého množství dat podobných typů (např. REAL, INTEGER) či většího množství dat různých typů (např. REAL, INTEGER, COMPLEX, CHARACTER) je užití MPI PACK a MPI UNPACK nejlepším řešením.
2.7
Příkazy MPI REDUCE a MPI ALLREDUCE
Příkazy MPI REDUCE a MPI ALLREDUCE umožňují provést určité typy operací (např. sčítání či násobení) na proměnných („mezivýsledcích”), které se nachází na jednotlivých procesech. Zároveň umožňují výsledek těchto operací buď poslat na určitý proces (MPI REDUCE) či rozeslat na všechny procesy (MPI ALLREDUCE). Ukažme si opět použití těchto příkazů v paralelním programu pro numerickou integraci lichoběžníkovou metodou. Připomeňme si, že paralelní program pro numerickou integraci lichoběžníkovou metodou: (i) rozdělil celkový interval na podintervaly podle celkového počtu procesů; (ii) jednotlivé procesy spočetly integrály pro jednotlivé podintervaly pomocí lichoběžníkové metody a (iii) tyto integrály se poslaly pomocí MPI SEND a MPI RECV na proces mající pořadí 0, kde se sečetly a výsledek se vytiskl. Bod (iii) nyní nahradíme příkazem MPI REDUCE a modifikovaná verze paralelního programu pro numerickou integraci lichoběžníkovou metodou může vypadat následovně: ! ! Trapezoidal Rule, MPI version 6. ! ! Algorithm: ! 1. Each process calculates "its" interval of integration. ! 2. Each process estimates the integral of f(x) over its interval ! using the trapezoidal rule. ! 3a. MPI_REDUCE sums the integrals on 0. ! 3b. Process 0 prints the result. ! program par_lich_int6 implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr integer :: n, & ! # of trapezoids local_n ! # of trapezoids for a processor !
27
real :: a, & ! left endpoint b, & ! right endpoint h, & ! trapezoid base length local_a, & ! left endpoint for a processor local_b, & ! right endpoint for a processor integral, & ! integral over a processor total ! total integral ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! get data ! call Get_Data2(a,b,n,myrank) ! if(mod(n,nprocs) /= 0) stop "par_lich_int6: Wrong n/nprocs ratio!" h=(b-a)/REAL(n) ! h is the same for all processes local_n=n/nprocs ! # of trapezoids ! ! length of interval of integration for a process ! local_a=a+REAL(myrank)*REAL(local_n)*h local_b=local_a+REAL(local_n)*h ! ! calculate integral on a process ! call Trap(local_a,local_b,local_n,h,integral) ! ! add up the integrals calculated by each process using MPI_REDUCE ! call MPI_REDUCE(integral,total,1,MPI_REAL,MPI_SUM,0, & MPI_COMM_WORLD,ierr) ! ! print the results ! if(myrank == 0) then write(*,’(1x,a,e13.5)’) " Integral =",total endif ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program par_lich_int6 ! ! Subroutine trapezoid for a processor ! subroutine Trap(local_a,local_b,local_n,h,integral) implicit none ! integer :: local_n,i real :: local_a,local_b,h,integral,x,f ! integral=(f(local_a)+f(local_b))/2.0 x=local_a do i=1,local_n-1 x=x+h
28
integral=integral+f(x) enddo integral=integral*h end subroutine Trap ! ! Algorithm: ! 1. Process 0 reads data. ! 2. Process 0 broadcasts data to other processes. ! subroutine Get_Data2(a,b,n,myrank) implicit none include ’mpif.h’ ! preprocessor directive ! integer :: myrank, & ! my process rank ierr integer :: n ! # of trapezoids ! real :: a, & ! left endpoint b ! right endpoint ! if(myrank == 0) then print *," a:" read*,a print *," b:" read*,b print *," n:" read*,n endif call MPI_BCAST(a,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(b,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) end subroutine Get_Data2 ! ! Function for integration ! function f(x) implicit none ! real :: f,x f=x*x end function f
Argument MPI SUM v příkazu MPI REDUCE říká, že se má provést součet proměnných integral z jednotlivých procesů do proměnné total. Šestý argument 0 v příkazu MPI REDUCE říká, že se výsledek total (součet proměnných integral ) má poslat na proces s pořadím 0. Obecný tvar příkazu MPI REDUCE lze obecně zapsat jako:
call MPI_REDUCE(operand,result,count,MPI_data_type,operation,root, & comm,ierr)
kde význam argumentů kromě operand, result, operation a root je zřejmý z předchozích probíraných MPI příkazů. • operand je proměnná, s kterou MPI REDUCE provádí určité typy operací. 29
• result je proměnná typu MPI data type, která obsahuje výsledek určitého typu operace na proměnných operand. • operation říká, jaký typ operace se provádí:11 operation Význam P MPI SUM součet, Q MPI PROD násobení, MPI MAX maximum MPI MIN minimum MPI MAXLOC maximum a pořadí procesu, kde se maximum nachází MPI MINLOC minimum a pořadí procesu, kde se minimum nachází • root je typu INTEGER a vyjadřuje pořadí procesu, na který se result posílá. Kdybychom použili v předchozím programu příkaz MPI ALLREDUCE místo příkazu MPI REDUCE, pak by se total rozeslal na všechny procesy. Je tedy zřejmé, že syntaxe příkazu MPI ALLREDUCE je totožná se syntaxí příkazu MPI REDUCE s tím, že neobsahuje argument root.
2.8
Často používané MPI příkazy
Uveďme si pro přehlednost probrané MPI příkazy a doplňme si je o nové. 2.8.1
Příkazy pro vytvoření a správu paralelního prostředí
A. Inicializace MPI
call MPI_INIT(ierr)
B. Ukončení MPI
call MPI_FINALIZE(ierr)
C. Definování komunikátoru a určení celkového počtu procesů
call MPI_COMM_SIZE(comm,nprocs,ierr)
11
Uvádíme jen nejpoužívanější operace.
30
D. Definování pořadí jednotlivých procesů
call MPI_COMM_RANK(comm,myrank,ierr)
E. Zastavení paralelního programu
call MPI_ABORT(comm,errorcode,ierr)
Příkaz MPI ABORT je paralelní verzí příkazu programovacího jazyka STOP. Použije-li jakýkoliv proces příkaz MPI ABORT (např. dojde-li k chybě během čtení na procesu), ukončí se běh paralelního programu na všech procesech.
2.8.2
Příkazy pro kolektivní komunikaci
A. Synchronizace či blokace paralelního programu
call MPI_BARRIER(comm,ierr)
B. Rozeslání informace na jednotlivé procesy
call MPI_BCAST(broadcast_message,count,MPI_data_type,root,comm,ierr)
C. Distribuce informací stejné velikosti na jednotlivé procesy
call MPI_SCATTER(sendbuf,sendcount,send_MPI_data_type, & recvbuf,recvcount,recv_MPI_data_type,root,comm,ierr)
Funkce příkazu je vysvětlena na obr. 6 a použití MPI SCATTER demonstruje následující program. program scatter implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr integer :: i ! real :: sendmsg(4), & ! send message recvmsg ! recieve message ! ! start up MPI !
31
call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatter: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! if(myrank == 0) then do i=1,nprocs sendmsg(i)=real(i) enddo endif ! call MPI_SCATTER(sendmsg,1,MPI_REAL,recvmsg,1,MPI_REAL,0, MPI_COMM_WORLD,ierr) ! print*," myrank =",myrank," recvmsg =",recvmsg ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program scatter
&
Obrázek 6: Schematický popis funkce příkazu MPI SCATTER.
D. Shromáždění informací stejné velikosti z jednotlivých procesů
call MPI_GATHER(sendbuf,sendcount,send_MPI_data_type, & recvbuf,recvcount,recv_MPI_data_type,root,comm,ierr)
Funkce příkazu je vysvětlena na obr. 7 a použití MPI GATHER demonstruje následující program. 32
program gather implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr ! real :: sendmsg, & ! send message recvmsg(4) ! recieve message ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gather: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! sendmsg=real(myrank+1) ! call MPI_GATHER(sendmsg,1,MPI_REAL,recvmsg,1,MPI_REAL,0, MPI_COMM_WORLD,ierr) ! if(myrank == 0) then print*," recvmsg:",recvmsg endif ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program gather
&
processes in comm need to call this routine.
Obrázek 7: Schematický popis funkce příkazu MPI GATHER.
E. Distribuce informací nestejné velikosti na jednotlivé procesy
33
call MPI_SCATTERV(sendbuf,sendcounts,displs,send_MPI_data_type, recvbuf,recvcount, recv_MPI_data_type, root,comm,ierr)
& &
Funkce příkazu je vysvětlena na obr. 8 a použití MPI SCATTERV demonstruje následující program.
program scatterv implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr integer :: send_count(0:3), & ! send count recv_count, & ! receive count displ(0:3) ! displacement ! real :: sendmsg(10), & ! send message recvmsg(4) ! recieve message ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatterv: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! if(myrank == 0) then sendmsg=(/1.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,4.0,4.0/) send_count=(/1,2,3,4/) displ=(/0,1,3,6/) endif recv_count=myrank+1 ! call MPI_SCATTERV(sendmsg,send_count,displ,MPI_REAL, recvmsg,recv_count, MPI_REAL, 0,MPI_COMM_WORLD,ierr) ! print*," myrank =",myrank," recvmsg:",recvmsg ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program scatterv
34
& &
Obrázek 8: Schematický popis funkce příkazu MPI SCATTERV. F. Shromáždění informací nestejné velikosti z jednotlivých procesů
call MPI_GATHERV(sendbuf,sendcount, send_MPI_data_type, recvbuf,recvcounts,displs,recv_MPI_data_type, root,comm,ierr)
& &
Funkce příkazu je vysvětlena na obr. 9 a použití MPI GATHERV demonstruje následující program. program gatherv implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,i integer :: send_count, & ! send count recv_count(0:3), & ! receive count displ(0:3) ! displacement ! real :: sendmsg(4), & ! send message recvmsg(10) ! recieve message ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gatherv: nprocs > 4!" ! ! get my process rank !
35
call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! do i=1,myrank+1 sendmsg(i)=myrank+1 enddo send_count=myrank+1 ! recv_count=(/1,2,3,4/) displ=(/0,1,3,6/) ! call MPI_GATHERV(sendmsg,send_count, MPI_REAL, recvmsg,recv_count,displ,MPI_REAL, 0,MPI_COMM_WORLD,ierr) ! if(myrank == 0) then print*," recvmsg:",recvmsg endif ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program gatherv
& &
Obrázek 9: Schematický popis funkce příkazu MPI GATHERV.
Příkazy MPI GATHER a MPI GATHERV mají ještě „ALL” verze MPI ALLGATHER a MPI ALLGATHERV, které shromažďují informace na všechny procesy a tudíž neobsahují argument root. G. Posílání informací stejné velikosti z jednotlivých procesů na všechny procesy
call MPI_ALLTOALL(sendbuf,sendcount,send_MPI_data_type, recvbuf,recvcount,recv_MPI_data_type, comm,ierr)
36
& &
Funkce příkazu je vysvětlena na obr. 10 a použití MPI ALLTOALL demonstruje následující program.
Obrázek 10: Schematický popis funkce příkazu MPI ALLTOALL. program alltoall implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,i ! real :: sendmsg(4), & ! send message recvmsg(4) ! recieve message ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gather: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! do i=1,nprocs sendmsg(i)=REAL(i+nprocs*myrank) enddo print*," myrank =",myrank," sendmsg:",sendmsg ! call MPI_ALLTOALL(sendmsg,1,MPI_REAL, & recvmsg,1,MPI_REAL, & MPI_COMM_WORLD,ierr) ! print*," myrank =",myrank," recvmsg:",recvmsg ! ! shut down MPI call MPI_FINALIZE(ierr) end program alltoall
37
H. Posílání informací nestejné velikosti z jednotlivých procesů na všechny procesy
call MPI_ALLTOALLV(sendbuf,sendcounts,sdispls,send_MPI_data_type, recvbuf,recvcounts,rdispls,recv_MPI_data_type, comm,ierr)
& &
Funkce příkazu je vysvětlena na obr. 11 a použití MPI ALLTOALLV demonstruje následující program. program alltoallv implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,i integer :: scnt(0:3), & ! send counts sdsp(0:3), & ! send displacements rcnt(0:3), & ! recv counts rdsp(0:3) ! recv displacements real :: sendmsg(10),recvmsg(16) ! sendmsg=(/1.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,4.0,4.0/) scnt=(/1,2,3,4/) sdsp=(/0,1,3,6/) ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gather: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! sendmsg ! do i=1,10 sendmsg(i)=sendmsg(i)+REAL(nprocs*myrank) enddo print*," myrank =",myrank," sendmsg:",sendmsg ! ! rcnt and rdsp ! do i=0,nprocs-1 rcnt(i)=myrank+1 rdsp(i)=i*(myrank+1) enddo ! call MPI_ALLTOALLV(sendmsg,scnt,sdsp,MPI_REAL, & recvmsg,rcnt,rdsp,MPI_REAL, & MPI_COMM_WORLD,ierr) !
38
print*," myrank =",myrank," recvmsg:",recvmsg ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program alltoallv
Obrázek 11: Schematický popis funkce příkazu MPI ALLTOALLV.
2.8.3
Příkazy pro operace na proměnných distribuovaných na jednotlivých procesech
A. Příkaz „Reduce”
call MPI_REDUCE(operand,result,count,MPI_data_type,operation,root, & comm,ierr)
Příkaz MPI REDUCE má ještě „ALL” verzi, které neobsahuje argument root. B. Příkaz „Scan”
call MPI_SCAN(sendbuf,recvbuf,count,MPI_data_type,operation,comm,ierr)
Funkce příkazu je vysvětlena na obr. 12 a použití MPI SCAN demonstruje následující program. 39
program scan implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr real :: sendmsg,recvmsg ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gather: nprocs > 4!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! sendmsg ! sendmsg=REAL(myrank+1) ! call MPI_SCAN(sendmsg,recvmsg,1,MPI_REAL,MPI_SUM, MPI_COMM_WORLD,ierr) ! print*," myrank =",myrank," recvmsg:",recvmsg ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program scan
&
Obrázek 12: Schematický popis funkce příkazu MPI SCAN. C. Příkaz „Reduce Scatter”
call MPI_MPI_REDUCE_SCATTER(sendbuf,recvbuf,recvcounts,MPI_data_type, operation,comm,ierr)
40
&
Funkce příkazu je vysvětlena na obr. 13 a použití MPI REDUCE SCATTER demonstruje následující program.
Obrázek 13: Schematický popis funkce příkazu MPI REDUCE SCATTER. program reduce_scatter implicit none include ’mpif.h’ ! preprocessor directive ! integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,i integer :: rcnt(0:3) ! recv counts ! real :: sendmsg(10), & ! send message recvmsg(4) ! recieve message ! rcnt=(/1,2,3,4/) ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatter: nprocs > 4!" !
41
! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! do i=1,10 sendmsg(i)=real(i+myrank*10) enddo print*," myrank =",myrank," sendmsg:",sendmsg ! call MPI_REDUCE_SCATTER(sendmsg,recvmsg,rcnt,MPI_REAL,MPI_SUM, MPI_COMM_WORLD,ierr) ! print*," myrank =",myrank," recvmsg =",recvmsg ! ! shut down MPI ! call MPI_FINALIZE(ierr) end program reduce_scatter
42
&
3
APLIKACE
3.1
Paralelní „tempering”
Paralelní „tempering” je způsob jak zefektivnit, urychlit či zlepšit vzorkování např. v Monte Carlo (MC) metodě.12 Ukážeme si princip paralelního „temperingu” na jednoduchém případu částice v jednorozměrném silovém poli.
3.1.1
Částice v jednorozměrném silovém poli
Představme si, že máme částici v jednorozměrném silovém poli. Částice má teplotu T a silové pole je charakterizované potenciálem U (x):
U (x) = = = = = = =
∞ 1 · [1 + sin (2πx)] 2 · [1 + sin (2πx)] 3 · [1 + sin (2πx)] 4 · [1 + sin (2πx)] 5 · [1 + sin (2πx)] ∞
x < −2 − 2 ≤ x ≤ −1.25 − 1.25 ≤ x ≤ −0.25 − 0.25 ≤ x ≤ 0.75 0.75 ≤ x ≤ 1.75 1.75 ≤ x ≤ 2 x>2
(10)
Průběh U (x) je graficky zobrazen na obr. 14. Potenciál U (x) je charakterizován čtyřmi lokálními minimy oddělenými bariérami vzrůstající velikosti ve směru kladné osy x.
12
Význam anglického výrazu „tempering” je zřejmý z dalšího výkladu.
43
10
8
U(x)
6
4
2
0 -2
-1
0
1
2
x
Obrázek 14: Průběh potenciálu U (x).
3.1.2
Monte Carlo metoda
„Pohyb” částice v intervalu x ∈ h−2; 2i (vzorkování fázového prostoru) budeme realizovat pomocí MC metody. MC metoda používá následující algoritmus: 1. Nechť se částice nachází v poloze xo a má zde hodnotu potenciálu U (xo ). 2. Vygenerujme novou pozici xn příkazem xn = xo + (2ζ − 1) ∆xmax
(11)
kde ζ je náhodné číslo z intervalu (0; 1) a ∆xmax je maximální posunutí ve směru x. Současně spočtěme hodnotu potenciálu v xn , U (xn ). 3. Spočtěme faktor exp {−β [U (xn ) − U (xo )]}
(12)
exp {−β [U (xn ) − U (xo )]} ≥ ζ
(13)
a proveďme test zda V rovnicích (12) a (13) β = 1/(kB T ) a kB je Boltzmannova konstanta. Pokud je podmínka (13) splněna, pak nahradíme „starou” pozici xo „novou” pozicí xn . Pokud podmínka (13) není splněna, „pohyb” se nepřijme a pozice xo se nezmění. 44
Výše popsaný postup se opakuje dostatečně dlouho než dojde k náležitému „provzorkování” fázového prostoru. MC procesem spočteme pravděpodobnost nalezení částice v místě x, P (x), pro různé teploty T . P (x) určíme následujícím způsobem: 1. Nadefinujeme velikost BIN u pro histogram, ∆x. 2. Příkazem BIN = INT ((xo − xmin )/∆x) + 1
(14)
spočteme pořadí BIN u v histogramu. V rovnici (14) je BIN deklarován jako INTEGER, xmin je minimální hodnota x, kterou částice může dosáhnout (v našem případě xmin = −2) a INT je vnitřní FORTRANská funkce, jež převádí čísla typu REAL na čísla typu INTEGER. 3. Akumulujeme histogram příkazem histogram(BIN ) = histogram(BIN ) + 1
(15)
4. Po skončení MC procesu normalizujeme histogram celkovým počtem MC pokusů. MC program pro částici v jednorozměrném silovém poli může vypadat následovně: ! ! A Single Particle on a 1D Potential with Multiple Local Minima ! program SerPT implicit none integer,parameter :: maxhist=1000 real(8),parameter :: x_min=-2.0d0,x_max=2.0d0 integer :: idum,ncyclus,nprint, & nacp_loc,ntri_loc,nc,ih, & histogram(maxhist) real(8) :: temp,beta, & dxmax,x,Ux, & dx,rat_loc ! ! Auxiliary Variables ! dxmax=0.1d0 dx=0.025d0 ! ! Init idum ! idum=-471 ! ! Input ! print*," # MC Cyclus:" read(*,*) ncyclus print*," Print Frequency:" read(*,*) nprint ! ! Iput Values (Hard Wired) !
45
x=-1.2d0 ! Initial x (1st valley) temp=2.00d0 ! kB.T beta=1.0d0/temp ! 1/(kB.T) ! ! Zero Accumulators ! nacp_loc=0 ntri_loc=0 histogram=0 ! ! Initial Potential ! call potential(x,Ux) ! ! Begin Loop Over Cyclus ! do nc=1,ncyclus ! ! Local Updates ! call update_loc(beta,dxmax,x,x_min,x_max,Ux, & dx,histogram,idum,nacp_loc) ntri_loc=ntri_loc+1 ! ! Output Every ’nprint’ Steps ! if(mod(nc,nprint) == 0) then rat_loc=dble(nacp_loc)/dble(ntri_loc) print*," # MC Cylus =",nc," Rat_Loc =",rat_loc print*," x =",x," U(x) =",Ux endif enddo ! ! Final Output ! x=x_min+(dx/2.0d0) do ih=1,maxhist write(10,"(1x,e13.5,1x,e13.5)") x,dble(histogram(ih))/dble(ntri_loc) x=x+dx if(x > x_max) exit enddo stop " SerPT: End of Calc!" end program SerPT ! ! Local Updates ! subroutine update_loc(beta,dxmax,x,x_min,x_max,Ux, & dx,histogram,idum,nacp_loc) implicit none integer :: idum,nacp_loc,bin, & histogram(*) real(8) :: random real(8) :: beta,dxmax,x,x_min,x_max,Ux,dx, & x_new,Ux_new,DelUx,beta_DelUx ! ! Accumulate Histogram ! bin=INT((x-x_min)/dx)+1 histogram(bin)=histogram(bin)+1 ! ! Move Particle ! x_new=x+(2.0d0*random(idum)-1.0d0)*dxmax ! ! Check Boundaries
46
! if(x_new < x_min) return if(x_new > x_max) return ! ! Calculate a New Value of Potential ! call potential(x_new,Ux_new) ! ! Check for Acceptance ! DelUx=Ux_new-Ux beta_DelUx=beta*DelUx if(dexp((-1.0d0)*beta_DelUx) > random(idum)) THEN ! ! Bookkeeping ! Ux=Ux_new x=x_new nacp_loc=nacp_loc+1 endif end subroutine update_loc ! ! Potential ! subroutine potential(x,Ux) implicit none real(8) :: pi,x,Ux ! pi=dacos(-1.0d0) ! ! Calculate Potential at x ! if(x <= -1.25d0) then Ux=1.0d0*(1.0d0+sin(2.0d0*pi*x)) else if(x <= -0.25d0) then Ux=2.0d0*(1.0d0+sin(2.0d0*pi*x)) else if(x <= 0.75d0) then Ux=3.0d0*(1.0d0+sin(2.0d0*pi*x)) else if(x <= 1.75d0) then Ux=4.0d0*(1.0d0+sin(2.0d0*pi*x)) else Ux=5.0d0*(1.0d0+sin(2.0d0*pi*x)) endif end subroutine potential ! ! RANDOM GENERATOR (Numerical Recipes in Fortran90) ! FUNCTION random(idum) ! ! "Minimal" random number generator of Park and Miller combined with a Marsaglia shift ! sequence. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint ! values). This fully portable, scalar generator has the "traditional" (not Fortran 90) calling ! sequence with a random deviate as the returned function value: call with idum a negative ! integer to initialize; thereafter, do not alter idum except to reinitialize. The period of this ! generator is about 3.1 x 10^18. ! IMPLICIT NONE ! ! Declaration of Variables ! INTEGER,PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B),INTENT(INOUT) :: idum REAL :: ran INTEGER(K4B),PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836 REAL,SAVE :: am
47
INTEGER(K4B),SAVE :: ix=-1,iy=-1,k ! real(8) :: random ! ! Initialize ! if((idum <= 0).or.(iy < 0)) then am=nearest(1.0,-1.0)/IM iy=ior(ieor(888889999,abs(idum)),1) ix=ieor(777755555,abs(idum)) ! ! Set idum Positive ! idum=abs(idum)+1 endif ! ! Marsaglia Shift Sequence with Period 2^32 - 1 ! ix=ieor(ix,ishft(ix,13)) ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) ! ! Park-Miller Sequence by Schrage’s Method, Period 2^31 - 2 ! k=iy/IQ iy=IA*(iy-k*IQ)-IR*k if(iy < 0) iy=iy+IM ! ! Combine the Two Generators with Masking to Ensure Nonzero Value ! ran=am*ior(iand(IM,ieor(ix,iy)),1) ! random=dble(ran) END FUNCTION random
Částice by měla navštívit všechna čtyři lokální minima přibližně stejně často, neboť se nacházejí v nulové potenciálové hladině. Obr. 15 ukazuje P (x) získané z 1 · 108 MC pokusů. V MC procesu bylo použito ∆xmax = 0.1 a ∆x = 0.025, MC proces začal s částicí v xini = −1.2 a byl proveden pro tři teploty: kB T = 2, kB T = 0.3 a kB T = 0.05. Z obr. 15 je vidět, že částice „navštěvuje” všechna čtyři lokální minima pouze pro nejvyšší teplotu kB T = 2. Pro nižší teploty částice vzhledem k bariérám není schopna dosáhnout všechna čtyři lokální minima a „navštěvuje” jen např. 1. a 2. lokální minimum v případě kB T = 0.3 či zůstává jen v 1. lokálním minimum v případě kB T = 0.05. 48
0.3
k B .T=2.0 k B .T=0.3 k B .T=0.05
P(x)
0.2
0.1
0.0 -2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
x
Obrázek 15: Pravděpodobnost nalezení částice v místě x, P (x), pro teploty: kB T = 2, kB T = 0.3 a kB T = 0.05; kB je Boltzmannova konstanta. P (x) byla získána z 1 · 108 Monte Carlo pokusů, v MC procesu bylo použito ∆xmax = 0.1 a ∆x = 0.025 a Monte Carlo proces začal s částicí v xini = −1.2.
3.1.3
Metoda paralelního „temperingu”
Použitím principu paralelního „temperingu” lze zefektivnit MC vzorkování při nízkých teplotách. Uvažujeme, že provádíme paralelně MC procesy při různých teplotách. Jednotlivé MC procesy C1 , C2 , . . ., Cn seřadíme tak, že T1 > T2 > ... > Tn (β1 < β2 < ... < βn ) Základní myšlenkou paralelního „temperingu” je provádět vedle standardních MC „pohybů” („lokálních” pohybů) vzájemnou výměnu „studenějších” a „teplejších” konfigurací („globální” pohyby). To se provede následujícím způsobem: 1. V kopiích {Ci }ni=1 probíhá standardní MC proces („lokální” pohyby). 2. Po určitém počtu MC kroků se proces zastaví a provede se pokus o „globální” pohyb: • vybere se náhodně jedna konfigurace Ci z {Ci }n−1 i=1 ; 49
• když exp (∆) ≥ ζ
(16)
provede se výměna konfigurací Ci ↔ Ci+1 . V rovnici (16) je ∆ = ∆β∆U ∆β = βi+1 − βi ∆U = Ui+1 − Ui
(17) (18) (19)
Paralelizace „globálních” pohybů pomocí MPI lze provést následujícím způsobem: 1. Proces s pořadím myrank == 0 vygeneruje náhodně jednu konfiguraci Ci z {Ci }n−1 i=1 příkazem i = IN T (ζ(n − 1)) + 1 (20) kde n odpovídá počtu paralelních procesů nprocs. 2. Proces s pořadím myrank == 0 rozešle i na ostatní procesy příkazem MPI BCAST. (MPI BCAST zároveň provede synchronizaci všech MC procesů.) 3. Procesy s pořadím myrank /= i − 1 a myrank /= i pokračují dále v „lokálních” pohybech. 4. Proces s pořadím myrank == i pošle na proces s myrank == i − 1 hodnotu U (Ci+1 ) pomocí příkazů MPI SEND a MPI RECV. 5. Na procesu s pořadím myrank == i − 1 se provede test exp (∆) ≥ ζ a výsledek testu se pošle na proces s pořadím myrank == i opět pomocí příkazů MPI SEND a MPI RECV. 6. Pokud je podmínka v předchozí rovnici splněna, provede se výměna konfigurací Ci ↔ Ci+1 t.j. x(Ci ) ↔ x(Ci+1 ) a U (Ci ) ↔ U (Ci+1 ) mezi procesy s pořadím myrank == i a s pořadím myrank == i − 1 pomocí příkazů MPI SEND a MPI RECV. Program pro paralelní „tempering” částice v jednorozměrném silovém poli může vypadat následovně: ! ! A Single Particle on a 1D Potential with Multiple Local Minima ! program ParPT implicit none ! include ’mpif.h’ ! preprocessor directive integer :: nprocs, & ! # of processes myrank, & ! my process rank
50
ierr ! integer,parameter :: maxhist=1000,maxtemp=10,nptemp=50 real(8),parameter :: x_min=-2.0d0,x_max=2.0d0 integer :: idum,ncyclus,nprint,nc,ih,notemp,iunit, & nacp_loc,ntri_loc,nacp_pt,ntri_pt, & histogram(maxhist) real(8) :: temp,beta, & dxmax,x,Ux, & dx,rat_loc,rat_glo, & tempar(maxtemp) character(len=7) :: cnfile character(len=7),dimension(maxtemp) :: outfil ! ! File Names ! outfil(1) ="f01.out";outfil(2) ="f02.out";outfil(3) ="f03.out" outfil(4) ="f04.out";outfil(5) ="f05.out";outfil(6) ="f06.out" outfil(7) ="f07.out";outfil(8) ="f08.out";outfil(9) ="f09.out" outfil(10)="f10.out" ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > maxtemp) stop "ParPT: # of processes > max # of temperatures!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! Define File Unit and Open Ooutput Files ! iunit=10+myrank*5 cnfile=outfil(myrank+1) open(unit=iunit,file=cnfile,status="unknown") ! ! Input ! if(myrank == 0) then print*," # MC Cyclus:" read(*,*) ncyclus print*," Print Frequency:" read(*,*) nprint endif ! ! broadcast input ! call MPI_BCAST(ncyclus,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nprint,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ! ! Temperatures kB.T; T1 > T2 > ... > Tn ! notemp=4 tempar(1)=2.00d0 tempar(2)=0.30d0 tempar(3)=0.10d0 tempar(4)=0.05d0 if(nprocs /= notemp) stop "ParPT: # of processes /= # of temperatures!" ! ! Auxiliary Variables !
51
dxmax=0.1d0 dx=0.025d0 ! ! Assign Actual Temperature ! temp=tempar(myrank+1) ! ! Init idum ! idum=-471 ! ! Iput Values (Hard Wired) ! x=-1.2d0 ! Initial x (1st valley) beta=1.0d0/temp ! 1/(kB.T) ! ! Zero Accumulators ! nacp_loc=0 nacp_pt=0 ntri_loc=0 ntri_pt=0 histogram=0 ! ! Initial Potential ! call potential(x,Ux) ! ! Begin Loop Over Cyclus ! do nc=1,ncyclus ! ! Local Updates ! call update_loc(beta,dxmax,x,x_min,x_max,Ux, & dx,histogram,idum,nacp_loc) ntri_loc=ntri_loc+1 ! ! Parallel Tempering After ’nptemp’ Cyclus ! if(mod(nc,nptemp) == 0) then call update_global(nprocs,myrank, & idum,beta,x,Ux,tempar,nacp_pt,ntri_pt) endif ! ! Output Every ’nprint’ Steps ! if(mod(nc,nprint) == 0) then rat_loc=dble(nacp_loc)/dble(ntri_loc) if(ntri_pt == 0) ntri_pt=1 rat_glo=dble(nacp_pt)/dble(ntri_pt) write(iunit,fmt="(1x,a,i15)") " # MC Cylus =",nc write(iunit,fmt="(1x,2(a,e13.5))") " Rat_Loc =",rat_loc, & " Rat_glo =",rat_glo write(iunit,fmt="(1x,2(a,e13.5))") " x =",x," U(x) =",Ux endif enddo ! ! Final Output ! x=x_min+(dx/2.0d0) do ih=1,maxhist write(iunit,"(1x,e13.5,1x,e13.5)") x,dble(histogram(ih))/dble(ntri_loc) x=x+dx if(x > x_max) exit
52
enddo close (unit=iunit) ! ! shut down MPI ! call MPI_FINALIZE(ierr) ! stop " ParPT: End of Calc!" end program ParPT ! ! Local Updates ! subroutine update_loc(beta,dxmax,x,x_min,x_max,Ux, dx,histogram,idum,nacp_loc) implicit none integer :: idum,nacp_loc,bin, & histogram(*) real(8) :: random real(8) :: beta,dxmax,x,x_min,x_max,Ux,dx, & x_new,Ux_new,DelUx,beta_DelUx ! ! Accumulate Histogram ! bin=INT((x-x_min)/dx)+1 histogram(bin)=histogram(bin)+1 ! ! Move Particle ! x_new=x+(2.0d0*random(idum)-1.0d0)*dxmax ! ! Check Boundaries ! if(x_new < x_min) return if(x_new > x_max) return ! ! Calculate a New Value of Potential ! call potential(x_new,Ux_new) ! ! Check for Acceptance ! DelUx=Ux_new-Ux beta_DelUx=beta*DelUx if(dexp((-1.0d0)*beta_DelUx) > random(idum)) THEN ! ! Bookkeeping ! Ux=Ux_new x=x_new nacp_loc=nacp_loc+1 endif end subroutine update_loc ! ! Potential ! subroutine potential(x,Ux) implicit none real(8) :: pi,x,Ux ! pi=dacos(-1.0d0) ! ! Calculate Potential at x ! if(x <= -1.25d0) then Ux=1.0d0*(1.0d0+sin(2.0d0*pi*x))
&
53
else if(x <= -0.25d0) then Ux=2.0d0*(1.0d0+sin(2.0d0*pi*x)) else if(x <= 0.75d0) then Ux=3.0d0*(1.0d0+sin(2.0d0*pi*x)) else if(x <= 1.75d0) then Ux=4.0d0*(1.0d0+sin(2.0d0*pi*x)) else Ux=5.0d0*(1.0d0+sin(2.0d0*pi*x)) endif end subroutine potential ! ! Global Updates ! subroutine update_global(nprocs,myrank, & idum,beta,x,Ux,tempar,nacp_pt,ntri_pt) implicit none ! include ’mpif.h’ ! preprocessor directive integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr integer :: status(MPI_STATUS_SIZE) integer :: idum,nacp_pt,ntri_pt,ii,ip1,itag,iacpt real(8) :: random real(8) :: beta,x,Ux,Ux1, & Delta_beta,Delta_U,acp_pt, & tempar(*) real(8) :: vectrn0(2),vectrn1(2) ! ! Choose ’i’ Randomly on "mastr" and then Broadcast ! if(myrank == 0) THEN ii=INT(dble(nprocs-1)*random(idum))+1 endif call MPI_BCAST(ii,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ii=ii-1 ! synchronize with ’myrank’ ! ! ’i+1’ ! ip1=ii+1 if((myrank /= ii).and.(myrank /= ip1)) return ! ! Send Potential from ’i+1’ on ’i’ ! itag=10 if(myrank == ip1) then call MPI_SEND(Ux,1,MPI_REAL8,ii,itag,MPI_COMM_WORLD,ierr) else if(myrank == ii) then call MPI_RECV(Ux1,1,MPI_REAL8,ip1,itag,MPI_COMM_WORLD,status,ierr) endif ! ! Check for Acceptance on ’i’: ! min[1,exp(DeltaBeta.DeltaU)] (Delta = ’i+1’ - ’i’) ! if(myrank == ii) then Delta_beta=(1.0d0/tempar(ip1+1))-beta Delta_U=Ux1-Ux acp_pt=dexp(Delta_beta*Delta_U) if(acp_pt > random(idum)) then iacpt=1 nacp_pt=nacp_pt+1 else iacpt=0 endif ntri_pt=ntri_pt+1
54
endif ! ! Send ’iacpt’ from ’i’ on ’i+1’ ! itag=20 if(myrank == ii) then call MPI_SEND(iacpt,1,MPI_INTEGER,IP1,itag,MPI_COMM_WORLD,ierr) else if(myrank == IP1) THEN call MPI_RECV(iacpt,1,MPI_INTEGER,ii,itag,MPI_COMM_WORLD,status,ierr) endif if(iacpt == 0) return ! ! Exchange Replicas between ’i’ and ’i+1’ ! ! Send Replica from ’i’ on ’i+1’ ! itag=30 if(myrank == ii) then vectrn0(1)=x vectrn0(2)=Ux call MPI_SEND(vectrn0,2,MPI_REAL8,ip1,itag,MPI_COMM_WORLD,ierr) else if(myrank == ip1) then call MPI_RECV(vectrn0,2,MPI_REAL8,ii,itag,MPI_COMM_WORLD,status,ierr) endif ! ! Send Replica from ’i+1’ on ’i’ ! itag=40 if(myrank == ip1) then vectrn1(1)=x vectrn1(2)=Ux call MPI_SEND(vectrn1,2,MPI_REAL8,ii,itag,MPI_COMM_WORLD,ierr) else if(myrank == ii) then call MPI_RECV(vectrn1,2,MPI_REAL8,ip1,itag,MPI_COMM_WORLD,status,ierr) endif ! ! Update Information ! if(myrank == ii) then ! ! on ’i’ ! x=vectrn1(1) Ux=vectrn1(2) else if(myrank == ip1) then ! ! on ’i+1’ ! x=vectrn0(1) Ux=vectrn0(2) endif end subroutine update_global ! ! RANDOM GENERATOR (Numerical Recipes in Fortran90) ! FUNCTION random(idum) ! ! "Minimal" random number generator of Park and Miller combined with a Marsaglia shift ! sequence. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint ! values). This fully portable, scalar generator has the "traditional" (not Fortran 90) calling ! sequence with a random deviate as the returned function value: call with idum a negative ! integer to initialize; thereafter, do not alter idum except to reinitialize. The period of this ! generator is about 3.1 x 10^18. ! IMPLICIT NONE
55
! ! Declaration of Variables ! INTEGER,PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B),INTENT(INOUT) :: idum REAL :: ran INTEGER(K4B),PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836 REAL,SAVE :: am INTEGER(K4B),SAVE :: ix=-1,iy=-1,k ! real(8) :: random ! ! Initialize ! if((idum <= 0).or.(iy < 0)) then am=nearest(1.0,-1.0)/IM iy=ior(ieor(888889999,abs(idum)),1) ix=ieor(777755555,abs(idum)) ! ! Set idum Positive ! idum=abs(idum)+1 endif ! ! Marsaglia Shift Sequence with Period 2^32 - 1 ! ix=ieor(ix,ishft(ix,13)) ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) ! ! Park-Miller Sequence by Schrage’s Method, Period 2^31 - 2 ! k=iy/IQ iy=IA*(iy-k*IQ)-IR*k if(iy < 0) iy=iy+IM ! ! Combine the Two Generators with Masking to Ensure Nonzero Value ! ran=am*ior(iand(IM,ieor(ix,iy)),1) ! random=dble(ran) END FUNCTION random
Obr. 16 ukazuje P (x) získané z 1 · 108 MC pokusů v každé kopii s pokusem o výměnu konfigurací každých 50 MC pokusů. V MC procesech bylo opět použito ∆xmax = 0.1 a ∆x = 0.025 a MC procesy začaly s částicí v xini = −1.2. MC procesy probíhaly při čtyřech teplotách: kB T = 2, kB T = 0.3, kB T = 0.1 a kB T = 0.05. Z obr. 16 je patrné, že částice „navštěvuje” všechna čtyři lokální minima přibližně stejně často pro všechny teploty. 56
Obrázek 16: Pravděpodobnost nalezení částice v místě x, P (x), pro teploty: kB T = 2, kB T = 0.3, kB T = 0.1 a kB T = 0.05; kB je Boltzmannova konstanta. P (x) byla získána z 1 · 108 Monte Carlo pokusů v každé kopii s pokusem o výměnu konfigurací každých 50 MC pokusů. V MC procesech bylo použito ∆xmax = 0.1 a ∆x = 0.025 a Monte Carlo procesy začaly s částicí v xini = −1.2.
3.2 3.2.1
Paralelní „tempering” při modelování polymerů Úvod
I „nechemici” se setkali s pojmem polymery. Nejrozšířenějším polymerem je polyethylen [CH2 − CH2 ]n . Jelikož je „velikost” vodíku výrazně menší než „velikost” uhlíku, lze polyethylen schematicky zapsat jako CH3 − ... − CH2 − CH2 − CH2 − ... − CH3
(21)
Z tohoto zápisu je vidět, že polymer je „řetízková” molekula tvořená tzv. segmenty (monomery). V případě polyethylenu segmenty odpovídají chemickým skupinám CH3 a CH2 . Podle tvaru řetízku se polymery dále dělí na lineární, rozvětvené atd. 57
3.2.2
Konformace polymeru
Jednou z nejdůležitějších informací o polymerech je znalost jejich konformačního chování. Pod pojmem konformační chování máme většinou na mysli vliv vnějšího prostředí (teploty, typu rozpouštědla či hustoty rozpouštědla) na tvar (konformaci) polymerů. Na obr. 17 jsou konformace modelového polymeru mající 20 segmentů při dvou teplotách, nízké a vysoké. Polymer je ve vakuu (bez vlivu rozpouštědla) a při nízké teplotě zaujme globulární konformace zatímco při vysoké teplotě má polymer tvar spirály („coilu”). Teplota, při které dochází ke změně konformace ze spirály do globule, se nazývá „coil-to-globule” přechodová teplota TC−G . TC−G je důležitou charakteristikou polymerů. Pro naše účely budeme konformaci polymerů charakterizovat dvěma parametry. Prvním je čtverec vzdálenosti konců polymeru, R2 , a druhým je čtverec poloměru setrvačnosti, Rg2 . R2 je definován jako
R2 = (r1 − rM )2 = (x1 − xM )2 + (y1 − yM )2 + (z1 − zM )2
(22)
kde r1 = (x1 , y1 , z1 ) a rM = (xM , yM , zM ) jsou polohové vektory resp. 1. a posledního (M ) segmentu polymeru. Rg2 je definován jako
Rg2 =
1
M X
mpolymeru
i=1
mi (ri − rcom )2
(23)
kde mpolymeru = M i=1 mi je hmotnost polymeru daná jako součet hmotností jednotlivých segmentů mi , ri = (xi , yi , zi ) je polohový vektor segmentů i a rcom = (xcom , ycom , zcom ) je polohový vektor středu hmotnosti polymeru. P
R2 a Rg2 zprůměrované přes velký soubor konformací, < R2 > a < Rg2 >, definují tzv. < parametr <=
hR2 i D
6 Rg2
E
(24)
Z teorie polymerů plyne, že pro globulární konformace je < < 1 a pro spirálové konformace je < ≥ 1.05. Polymer mající < ' 1 se nazývá „ideální”. TC−G je pak definována jako teplota, kde < = 1 58
Nízká teplota
Vysoká teplota
Obrázek 17: Konformace modelového polymeru mající 20 segmentů při nízké a vysoké teplotě. Polymer je ve vakuu a při nízké teplotě polymer zaujme globulární konformace zatímco při vysoké teplotě má polymer tvar spirály („coilu”).
3.2.3
Model polymeru
Budeme pracovat s velmi jednoduchým modelem polymeru. Segmenty polymeru budou reprezentovány tuhými koulemi o průměru σ, jejichž středy jsou vzdáleny od sebe opět o σ (dotýkající se tuhé koule). Segmenty vzdálené od sebe o více než jeden segment budou interagovat pomocí potenciálu pravoúhlé jámy uSW (r; σ, λ) = +∞ = − = 0
r<σ σ < r < λσ r > λσ
kde je hloubka (energie) pravoúhlé jámy a λ charakterizuje délku dosahu atrakce. 59
(25)
3.2.4
Rosenblutův Monte Carlo algoritmus
Budeme se zabývat jedním polymerem ve vakuu. „Pohyb” polymeru (vzorkování fázového prostoru) lze efektivně provést pomocí Rosenblutova MC algoritmu. Ten začíná náhodným výběrem segmentu i z intervalu < 1, M −1 > a následným spočtením Rosenblutových faktorů W (n) a W (o) pro novou a starou konformaci segmentů i + 1, i + 2, ..., M . Nová konformace se přijme, když W (n) >ζ (26) W (o) A. Výpočet W (n) W (n) počítáme po jednotlivých segmentech, přičemž začínáme od segmentu i + 1. Pro segment i + 1: 1. Vygenerujme k orientací, {bj }kj=1 , a spočteme {exp [−βui+1 (bj )]}kj=1
(27)
kde ui+1 (bj ) je energie na segment i + 1 od segmentů 1, 2, ..., i. 2. Vybereme z k orientací jednu, n, s pravděpodobností pi+1 (n) = kde wi+1 (n) = gram.
Pk
j=1
exp [−βui+1 (bn )] wi+1 (n)
(28)
exp [−βui+1 (bj )]. Algoritmus pro výběr n-té orientace viz. pro-
Tento postup zopakujeme pro segmenty i + 2, i + 3, ..., M a W (n) určíme jako W (n) =
M Y
wj (n) = wi+1 (n) · wi+2 (n) · ... · wM (n)
(29)
j=i+1
B. Výpočet W (o) W (o) počítáme obdobně jako W (n). Pro segment i + 1: 1. Vygenerujme k − 1 orientací,13 {bj }k−1 j=1 , a spočteme k−1 {exp [−βui+1 (bj )]}j=1 13
Máme k dispozici původní orientaci bo .
60
(30)
2. Spočteme wi+1 (o) = exp [−βui+1 (bo )] +
k−1 X
exp [−βui+1 (bj )]
(31)
j=1
Tento postup zopakujeme pro segmenty i + 2, i + 3, ..., M a W (o) určíme jako W (o) =
M Y
wj (o) = wi+1 (o) · wi+2 (o) · ... · wM (o)
(32)
j=i+1
3.2.5
Rosenblutovo Monte Carlo a paralelní „tempering”
Použitím paralelního „temperingu” lze dále zefektivnit Rosenblutovo MC a simulovat konformace polymerů i při nízkých teplotách. Uvažujeme opět, že provádíme paralelně Rosenblutovy MC procesy při různých teplotách. Jednotlivé Rosenblutovy MC procesy C1 , C2 , . . ., Cn seřadíme tak, že T1 > T2 > ... > Tn (β1 < β2 < ... < βn ) a budeme uvažovat, že počet procesů n je liché číslo. Dále bychom mohli pokračovat jako při paralelním „temperingu” částice v jednorozměrném silovém poli. Uveďme si však modifikaci paralelního „temperingu”, která ještě více zefektivní vzorkování konfiguračního prostoru. V této verzi nebudeme vybírat kopie pro výměnu náhodně, ale provedeme výměny pro všechny páry kopií následujícím způsobem: C1 ↔ C2 C2 ↔ C3
C3 ↔ C4 C4 ↔ C5
C5 ↔ C6 C6 ↔ C7
... Cn−2 ... Cn−1
↔ ↔
Cn−1 Cn
(33) (34)
přičemž výměnu párů kopií dle rovnic (33) a (34) postupně střídáme. Výměna páru konfigurací Ci ↔ Ci+1 se opět provede na základě kritéria exp (∆) ≥ ζ
(35)
Program Před vytvořením programu si všimněme, že: • Ve výpočtu nemá cenu pokračovat, pokud wj (n) během výpočtu W (n) nabyde hodnoty nula. • Energie polymeru v nové konfiguraci, utot , je utot = utot (o)1,2,...,i +
M X j=i+1
61
uj (bn )
(36)
Program pro Rosenblutovo MC a paralelní „tempering” polymeru ve vakuu pak může vypadat následovně: ! ! SW POLYMER CHAIN IN VACUUM WITH ’1,3 > INTRA’ ! Parallel Tempering ! module system_data implicit none ! integer,parameter :: NSITE=10 integer,parameter :: KMAX=10 integer,parameter :: KMAXQ=20 integer,parameter :: NBLOCK=500000 integer,parameter :: NBL=200 integer,parameter :: nptemp=50 ! # o Steps Between PT integer,parameter :: maxtemp=10 integer :: idum=-471 ! Seed for Random Generator ! real(8),parameter :: RLAMSW=1.5D0 end module system_data ! PROGRAM PT_SWCHAIN_VACCUM ! ! System Data from Module ! use system_data ! implicit none ! include ’mpif.h’ ! preprocessor directive integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr ! integer :: IMODE,NTRIAL_1,NTRIAL_2, & NPRINT,NSTEP,NTRIAL,NACPTR, & nacp_pt,ntri_pt,notemp real(8) :: RLAMSW_2,PE,PEFIN, & tempar(maxtemp),temp, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G real(8) :: XS(NSITE),YS(NSITE),ZS(NSITE), & BSMPE(NBL),BSMR2(NBL),BSMR2G(NBL) character(len=80) :: TITLE logical :: OVRLAP ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > maxtemp) & stop "PT_SWCHAIN_VACCUM: # of processes > max # of temperatures!" ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! INITIALISATION !
62
if(myrank == 0) then WRITE(*,’(’’ ENTER RUN TITLE ’’)’) READ(*,’(A)’) TITLE WRITE(*,’(’’ ENTER # OF EQL TRIALS ’’)’) READ(*,*) NTRIAL_1 WRITE(*,’(’’ ENTER # OF RUN TRIALS ’’)’) READ(*,*) NTRIAL_2 WRITE(*,’(’’ ENTER # OF STEPS BETWEEN OUTPUT LINES ’’)’) READ(*,*) NPRINT endif ! ! broadcast input ! call MPI_BCAST(TITLE,80,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(NTRIAL_1,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(NTRIAL_2,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(NPRINT,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ! ! Temperatures kB.T; T1 > T2 > ... > Tn ! notemp=5 ! Odd # if(mod(notemp,2) == 0) stop "PT_SWCHAIN_VACCUM: notemp Must Be Odd #!" tempar(1)=1.5d0 tempar(2)=1.0d0 tempar(3)=0.8d0 tempar(4)=0.6d0 tempar(5)=0.5d0 if(nprocs /= notemp) & stop "PT_SWCHAIN_VACCUM: # of processes /= # of temperatures!" ! ! Assign Actual Temperature ! temp=tempar(myrank+1) ! ! AUXILIARY VALUES ! RLAMSW_2=RLAMSW**2 IF(KMAX > KMAXQ) & STOP "PT_SWCHAIN_VACCUM: KMAX IS TOO BIG!" ! CHECK KMAX ! ! LABEL THE OUTPUT ! CALL OUTPUT(0,NBL,NBLOCK,KMAX,maxtemp,NSTEP, & NTRIAL,NACPTR,nacp_pt,ntri_pt, & TEMP,RLAMSW, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G, & TITLE, & myrank) ! ! SET INITIAL COORDINATES RANDOMLY ! CALL RANSTR(NSITE,idum,RLAMSW_2,XS,YS,ZS) ! ! CALCULATE INITIAL ENERGY ! CALL SUM_ENR(NSITE,RLAMSW_2,XS,YS,ZS,PE,OVRLAP) IF(OVRLAP) & STOP "PT_SWCHAIN_VACCUM: OVERLAP IN INITIAL CONFIGURATION!" if(myrank == 0) then WRITE(*,’(1X,A,E13.5)’) " INTIAL P.E. =",PE endif !
63
! INITIAL OUTPUT ! CALL OUTPUT(1,NBL,NBLOCK,KMAX,maxtemp,NSTEP, & NTRIAL,NACPTR,nacp_pt,ntri_pt, & TEMP,RLAMSW, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G, & TITLE, & myrank) ! ! ZERO ACCUMULATORS ! CALL ZERO(NBL,SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G) ! ! LOOP OVER MODE ! DO IMODE=1,2 ! ! ZERO COUNTERS ! NACPTR=0 nacp_pt=0 ntri_pt=0 ! ! BEGIN OF LOOP OVER TRIALS ! IF(IMODE.EQ.1) THEN NTRIAL=NTRIAL_1 ELSE NTRIAL=NTRIAL_2 ENDIF DO NSTEP=1,NTRIAL ! ! CONFORMATION USING CONFIGURATIONAL BIAS ! CALL CHCONF(NACPTR,NSITE,idum,KMAX,KMAXQ, & RLAMSW_2,TEMP,XS,YS,ZS,PE) if(IMODE == 2) then ! ! Parallel Tempering After ’nptemp’ Steps ! if(mod(nstep,nptemp) == 0) then call update_global(nprocs,myrank, & NSITE,idum,TEMP,tempar, & XS,YS,ZS,PE, & nacp_pt,ntri_pt) endif ! ! EVALUATE QUANTITIES OF INTEREST ! CALL QUANTS(NSTEP,NBLOCK,NBL,NSITE, & XS,YS,ZS,PE, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G) ! ! PERIODIC OUTPUT ! IF(MOD(NSTEP,NPRINT).EQ.0) THEN CALL OUTPUT(2,NBL,NBLOCK,KMAX,maxtemp,NSTEP, &
64
NTRIAL,NACPTR,nacp_pt,ntri_pt, TEMP,RLAMSW, STPPE,STPR2,STPR2G, SUMPE,SUMR2,SUMR2G, SSQPE,SSQR2,SSQR2G, BSMPE,BSMR2,BSMR2G, TITLE, myrank)
& & & & & & &
ENDIF endif ENDDO ENDDO ! ! END OF LOOP OVER TRIALS ! ! CALCULATE FINAL ENERGY ! CALL SUM_ENR(NSITE,RLAMSW_2,XS,YS,ZS,PEFIN,OVRLAP) IF(OVRLAP) & STOP "PT_SWCHAIN_VACCUM: OVERLAP IN FINAL CONFIGURATION!" if(myrank == 0) then WRITE(*,’(1X,A,E13.5)’) ’ FINAL P.E. =’,PEFIN ! ! FINAL CHECK OF ENERGY ! WRITE(*,’(1X,A,E13.5)’) ’ (MC P.E. - SUM P.E.) =’,DABS(PE-PEFIN) endif ! ! FINAL OUTPUT ! CALL OUTPUT(3,NBL,NBLOCK,KMAX,maxtemp,NSTEP, & NTRIAL,NACPTR,nacp_pt,ntri_pt, & TEMP,RLAMSW, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G, & TITLE, & myrank) ! ! shut down MPI ! call MPI_FINALIZE(ierr) ! STOP "PT_SWCHAIN_VACCUM: END OF CALCULATION!" END PROGRAM PT_SWCHAIN_VACCUM ! ! SET UP RANDOMLY SITES ! SUBROUTINE RANSTR(NSITE,idum,RLAMSW_2,XS,YS,ZS) implicit none integer,parameter :: ITMAX=1000 integer :: NSITE,idum,ISG,IT real(8) :: PE,RLAMSW_2,EX,EY,EZ real(8) :: XS(*),YS(*),ZS(*), & XSI(NSITE),YSI(NSITE),ZSI(NSITE) logical :: OVRLAP ! ! CHAIN RANDOMLY ! ! THE 1ST SEGMENT IN THE ORIGIN ! XS(1)=0.0d0 YS(1)=0.0d0 ZS(1)=0.0d0
65
! ! LOOP OVER REST OF SEGMENTS ! DO ISG=2,NSITE IT=0 do IT=IT+1 IF(IT > ITMAX) STOP "RANSTR: MAX INSERTION!" ! ! GENERATE RANDOM VECTOR ON UNIT SPHERE ! CALL SPH3E(EX,EY,EZ,idum) ! ! NEW POSITION OF THE SEGMENT ! XSI(ISG)=XSI(ISG-1)+EX YSI(ISG)=YSI(ISG-1)+EY ZSI(ISG)=ZSI(ISG-1)+EZ IF(ISG >= 3) THEN ! ! TEST OVRLAP ! CALL ENR(RLAMSW_2,ISG,XSI,YSI,ZSI,PE,OVRLAP) IF(OVRLAP) cycle ENDIF exit enddo ! ! SAVE POSITION ! XS(ISG)=XSI(ISG) YS(ISG)=YSI(ISG) ZS(ISG)=ZSI(ISG) ENDDO END SUBROUTINE RANSTR ! ! CONFORMATION USING CONFIGURATIONAL BIAS ! SUBROUTINE CHCONF(NACPTR,NSITE,idum,KMAX,KMAXQ, & RLAMSW_2,TEMP,XS,YS,ZS,PE) implicit none integer :: NACPTR,NSITE,idum,KMAX,KMAXQ,ISG,NGRWTH,K real(8) :: random real(8) :: RLAMSW_2,TEMP,PE,PENW,PEOL,WNW,WOL, & ACPROT real(8) :: XS(*),YS(*),ZS(*), & XSOL(NSITE),YSOL(NSITE),ZSOL(NSITE), & XSNW(NSITE),YSNW(NSITE),ZSNW(NSITE) logical :: NWCONF,OVRLAP ! ! SELECT RANDOMLY A SEGMENT FOR GROWTH ! NGRWTH=INT(DBLE(NSITE-1)*random(idum))+1 NGRWTH=NGRWTH+1 ! ! CALCULATE NEW ROSENBLUTH FACTOR ! NWCONF=.TRUE. ! ! SITE COORDINATES ! DO ISG=1,NGRWTH-1 XSNW(ISG)=XS(ISG) YSNW(ISG)=YS(ISG) ZSNW(ISG)=ZS(ISG)
66
ENDDO CALL GROWTH(NWCONF,NGRWTH,idum,NSITE,KMAX,KMAXQ, & TEMP,RLAMSW_2, & WNW,XSNW,YSNW,ZSNW,PENW,OVRLAP) IF(OVRLAP) RETURN ! ! CALCULATE OLD ROSENBLUTH FACTOR ! NWCONF=.FALSE. ! ! SITE COORDINATES ! DO ISG=1,NSITE XSOL(ISG)=XS(ISG) YSOL(ISG)=YS(ISG) ZSOL(ISG)=ZS(ISG) ENDDO CALL GROWTH(NWCONF,NGRWTH,idum,NSITE,KMAX,KMAXQ, & TEMP,RLAMSW_2, & WOL,XSOL,YSOL,ZSOL,PEOL,OVRLAP) ! ! CHECK FOR ACCEPTANCE: min(1,(W(n)/W(o))]} ! ACPROT=(WNW/WOL) IF(ACPROT > random(idum)) THEN ! ! BOOKKEEPING ! ! THERMO ! PE=PE+PENW-PEOL ! ! SITE COORDINATES ! DO ISG=NGRWTH,NSITE XS(ISG)=XSNW(ISG) YS(ISG)=YSNW(ISG) ZS(ISG)=ZSNW(ISG) ENDDO NACPTR=NACPTR+1 ENDIF END SUBROUTINE CHCONF ! ! ROSENBLUTH FACTOR FOR SW CHAIN ! SUBROUTINE GROWTH(NWCONF,NGRWTH,idum,NSITE,KMAX,KMAXQ, TEMP,RLAMSW_2, W,XSNW,YSNW,ZSNW,PE,OVRLAP) implicit none integer :: NGRWTH,idum,NSITE,KMAX,KMAXQ,K,ISG,N real(8),parameter :: EXPMAX=700.0D0 real(8) :: random real(8) :: TEMP,RLAMSW_2,PE,PEOLK,QPEK,EX,EY,EZ, & W,WiNW,WiOL,QARG,WiK real(8) :: XSNW(*),YSNW(*),ZSNW(*), & XSI(NSITE),YSI(NSITE),ZSI(NSITE), & XSK(KMAXQ),YSK(KMAXQ),ZSK(KMAXQ), & QPE(KMAXQ),Wi(KMAXQ) logical :: NWCONF,OVRLP1,OVRLAP ! ! AUXILIARY VARIABLES ! OVRLAP=.FALSE. ! ! SITE COORDINATES
67
& &
! DO ISG=1,NGRWTH-1 XSI(ISG)=XSNW(ISG) YSI(ISG)=YSNW(ISG) ZSI(ISG)=ZSNW(ISG) ENDDO ! ! OLD/NEW ROSENBLUTH FACTORS ! PE=0.0D0 W=1.0D0 IF(NWCONF) THEN ! ! CALCULATE NEW ROSENBLUTH FACTOR ! ! LOOP OVER SEGMENTS ! DO ISG=NGRWTH,NSITE ! ! LOOP OVER ORIENTATIONS ! WiNW=0.0D0 DO K=1,KMAX ! ! GENERATE RANDOM VECTOR ON UNIT SPHERE ! CALL SPH3E(EX,EY,EZ,idum) ! ! NEW POSITION OF THE SEGMENT ! XSI(ISG)=XSI(ISG-1)+EX YSI(ISG)=YSI(ISG-1)+EY ZSI(ISG)=ZSI(ISG-1)+EZ IF(ISG >= 3) THEN ! ! CALCULATE ENERGY ON THE SEGMENT ! CALL ENR(RLAMSW_2,ISG,XSI,YSI,ZSI,QPE(K),OVRLP1) ! ! CALCULATE w_i(n) FOR ISG > 2 ! IF(OVRLP1) THEN Wi(K)=0.0D0 ELSE QARG=QPE(K)/TEMP IF(QARG > EXPMAX) THEN Wi(K)=0.0D0 ELSE Wi(K)=DEXP((-1.0D0)*QARG) ENDIF ENDIF ELSE ! ! CALCULATE w_i(n) FOR ISG = 2 ! QPE(K)=0.0D0 Wi(K)=1.0D0 ENDIF WiNW=WiNW+Wi(K) ! ! SAVE SITE COORDINATES ! XSK(K)=XSI(ISG) YSK(K)=YSI(ISG) ZSK(K)=ZSI(ISG)
68
ENDDO ! ! UPDATE WNW ! IF(WiNW == 0.0D0) THEN OVRLAP=.TRUE. RETURN ENDIF W=W*WiNW ! ! SELECT A NEW ORIENTATION WITH PROBABILITY = w(k) / W(n) ! CALL SELECT(Wi,WiNW,idum,N) ! ! NEW SITE COORDINATES ! XSI(ISG)=XSK(N) YSI(ISG)=YSK(N) ZSI(ISG)=ZSK(N) ! XSNW(ISG)=XSI(ISG) YSNW(ISG)=YSI(ISG) ZSNW(ISG)=ZSI(ISG) ! ! UPDATE ENERGY ! PE=PE+QPE(N) ENDDO ELSE ! ! CALCULATE OLD ROSENBLUTH FACTOR ! ! LOOP OVER THE SEGMENTS ! DO ISG=NGRWTH,NSITE ! ! OLD ORIENTATION ! XSI(ISG)=XSNW(ISG) YSI(ISG)=YSNW(ISG) ZSI(ISG)=ZSNW(ISG) IF(ISG >= 3) THEN ! ! CALCULATE ENERGY ON THE SEGMENT ! CALL ENR(RLAMSW_2,ISG,XSI,YSI,ZSI,PEOLK,OVRLP1) IF(OVRLP1) THEN ! ! ERROR STOP ! STOP "GROWTH: OVRLAP IN OLD CONF.!" ENDIF ! ! CALCULATE w_i(o) FOR ISG > 2 ! QARG=PEOLK/TEMP IF(QARG > EXPMAX) THEN ! ! ERROR STOP ! STOP "GROWTH: OVRLAP IN OLD CONF.!" ENDIF WiOL=DEXP((-1.0D0)*QARG) ELSE !
69
! CALCULATE w_i(o) FOR ISG = 2 ! PEOLK=0.0D0 WiOL=1.0D0 ENDIF ! ! UPDATE ENERGY ! PE=PE+PEOLK ! ! LOOP OVER ORIENTATIONS ! DO K=2,KMAX ! ! GENERATE RANDOM VECTOR ON UNIT SPHERE ! CALL SPH3E(EX,EY,EZ,idum) ! ! NEW POSITION OF THE SEGMENT ! XSI(ISG)=XSI(ISG-1)+EX YSI(ISG)=YSI(ISG-1)+EY ZSI(ISG)=ZSI(ISG-1)+EZ IF(ISG >= 3) THEN ! ! ENERGY ! CALL ENR(RLAMSW_2,ISG,XSI,YSI,ZSI,QPEK,OVRLP1) ! ! CALCULATE w_i(o) FOR ISG > 2 ! IF(OVRLP1) THEN WiK=0.0D0 ELSE QARG=QPEK/TEMP IF(QARG> EXPMAX) THEN WiK=0.0D0 ELSE WiK=DEXP((-1.0D0)*QARG) ENDIF ENDIF ELSE ! ! CALCULATE w_i(o) FOR ISG = 2 ! WiK=1.0D0 ENDIF WiOL=WiOL+WiK ENDDO ! ! UPDATE WOL ! W=W*WiOL ! ! OLD ORIENTATION BACK ! XSI(ISG)=XSNW(ISG) YSI(ISG)=YSNW(ISG) ZSI(ISG)=ZSNW(ISG) ENDDO ENDIF END SUBROUTINE GROWTH ! ! Global Updates !
70
subroutine update_global(nprocs,myrank, & NSITE,idum,TEMP,tempar, & XS,YS,ZS,PE, & nacp_pt,ntri_pt) implicit none ! include ’mpif.h’ ! preprocessor directive integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr integer :: status(MPI_STATUS_SIZE) ! integer :: idum,NSITE,nacp_pt,ntri_pt,ii,ip1,itag,iacpt, & ninfo,j real(8) :: random real(8) :: temp,PE,PE_ii,Delta_beta,Delta_PE,acp_pt real(8) :: tempar(*), & XS(*),YS(*),ZS(*) real(8) :: vectrn0(3*NSITE+1),vectrn1(3*NSITE+1) ! ntri_pt=ntri_pt+1 call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! Synchronize All Processes ! ! Type of ’Switch’ ! if(mod(ntri_pt,2) == 0) then if(myrank == (nprocs-1)) return if(mod(myrank,2) == 0) then ii=myrank ! 0, 2, 4, ... ip1=myrank+1 ! 1, 3, 5, ... else ii=myrank-1 ! 0, 2, 4, ... ip1=myrank ! 1, 3, 5, ... endif else if(myrank == 0) return if(mod(myrank,2) /= 0) then ii=myrank ! 1, 3, 5, ... ip1=myrank+1 ! 2, 4, 6, ... else ii=myrank-1 ! 1, 3, 5, ... ip1=myrank ! 2, 4, 6, ... endif endif ! ! Send Potential from ’i+1’ on ’i’ ! itag=10 if(myrank == ip1) then call MPI_SEND(PE,1,MPI_REAL8,ii,itag,MPI_COMM_WORLD,ierr) else if(myrank == ii) then call MPI_RECV(PE_ii,1,MPI_REAL8,ip1,itag,MPI_COMM_WORLD, & status,ierr) endif ! ! Check for Acceptance on ’i’: ! min[1,exp(DeltaBeta.DeltaU)] (Delta = ’i+1’ - ’i’) ! if(myrank == ii) then Delta_beta=(1.0d0/tempar(ip1+1))-(1.0d0/temp) Delta_PE=PE_ii-PE acp_pt=dexp(Delta_beta*Delta_PE) if(acp_pt > random(idum)) then iacpt=1 nacp_pt=nacp_pt+1
71
else iacpt=0 endif endif ! ! Send ’iacpt’ from ’i’ on ’i+1’ ! itag=20 if(myrank == ii) then call MPI_SEND(iacpt,1,MPI_INTEGER,ip1,itag,MPI_COMM_WORLD,ierr) else if(myrank == ip1) THEN call MPI_RECV(iacpt,1,MPI_INTEGER,ii,itag,MPI_COMM_WORLD, & status,ierr) endif if(iacpt == 0) return ! ! Exchange Replicas between ’i’ and ’i+1’ ! ninfo=3*NSITE+1 ! ! Send Replica from ’i’ on ’i+1’ ! itag=30 if(myrank == ii) then do j=1,NSITE vectrn0(j)=XS(j) vectrn0(NSITE+j)=YS(j) vectrn0(2*NSITE+j)=ZS(j) enddo vectrn0(ninfo)=PE call MPI_SEND(vectrn0,ninfo,MPI_REAL8,ip1,itag,MPI_COMM_WORLD,ierr) else if(myrank == ip1) then call MPI_RECV(vectrn0,ninfo,MPI_REAL8,ii,itag,MPI_COMM_WORLD, & status,ierr) endif ! ! Send Replica from ’i+1’ on ’i’ ! itag=40 if(myrank == ip1) then do j=1,NSITE vectrn1(j)=XS(j) vectrn1(NSITE+j)=YS(j) vectrn1(2*NSITE+j)=ZS(j) enddo vectrn1(ninfo)=PE call MPI_SEND(vectrn1,ninfo,MPI_REAL8,ii,itag,MPI_COMM_WORLD,ierr) else if(myrank == ii) then call MPI_RECV(vectrn1,ninfo,MPI_REAL8,ip1,itag,MPI_COMM_WORLD, & status,ierr) endif ! ! Update Information ! if(myrank == ii) then ! ! On ’i’ ! do j=1,NSITE XS(j)=vectrn1(j) YS(j)=vectrn1(NSITE+j) ZS(j)=vectrn1(2*NSITE+j) enddo PE=vectrn1(ninfo) else if(myrank == ip1) then
72
! ! On ’i+1’ ! do j=1,NSITE XS(j)=vectrn0(j) YS(j)=vectrn0(NSITE+j) ZS(j)=vectrn0(2*NSITE+j) enddo PE=vectrn0(ninfo) endif end subroutine update_global ! ! CALCULATE ENERGY ON SW CHAIN SEGMENT ! SUBROUTINE ENR(RLAMSW_2,NSG,XSI,YSI,ZSI,PE,OVRLAP) implicit none integer :: NSG,ISG real(8) :: RLAMSW_2,PE, & XISG_NSG,YISG_NSG,ZISG_NSG,R2ISG_NSG real(8) :: XSI(*),YSI(*),ZSI(*) logical :: OVRLAP ! ! AUXILIARY VARIABLES ! OVRLAP=.FALSE. ! ! ZERO ACUMULATOR ! PE=0.0D0 ! ! GO OVER SEGMENTS ! do ISG=1,NSG-2 XISG_NSG=XSI(NSG)-XSI(ISG) YISG_NSG=YSI(NSG)-YSI(ISG) ZISG_NSG=ZSI(NSG)-ZSI(ISG) ! ! RANGE OF INTERACTION ! R2ISG_NSG=XISG_NSG**2+YISG_NSG**2+ZISG_NSG**2 IF(R2ISG_NSG < RLAMSW_2) THEN ! ! HS REPULSION ! IF(R2ISG_NSG < 1.0D0) THEN OVRLAP=.TRUE. RETURN ENDIF ! ! SW ATTRACTION ! PE=PE-1.0D0 ENDIF enddo END SUBROUTINE ENR ! ! CALCULATE INTRAMOLECULAR ENERGY FOR SW CHAIN ! SUBROUTINE SUM_ENR(NSITE,RLAMSW_2,XS,YS,ZS,PE,OVRLAP) implicit none integer :: NSITE,NSG,ISG real(8) :: RLAMSW_2,PE, & XISG_NSG,YISG_NSG,ZISG_NSG,R2ISG_NSG real(8) :: XS(*),YS(*),ZS(*) logical :: OVRLAP
73
! ! AUXILIARY VARIABLES ! OVRLAP=.FALSE. ! ! ZERO ACUMULATOR ! PE=0.0D0 ! ! GO OVER CHAIN ! DO NSG=3,NSITE ! ! GO OVER SEGMENTS ! do ISG=1,NSG-2 XISG_NSG=XS(NSG)-XS(ISG) YISG_NSG=YS(NSG)-YS(ISG) ZISG_NSG=ZS(NSG)-ZS(ISG) ! ! RANGE OF INTERACTION ! R2ISG_NSG=XISG_NSG**2+YISG_NSG**2+ZISG_NSG**2 IF(R2ISG_NSG < RLAMSW_2) THEN ! ! HS REPULSION ! IF(R2ISG_NSG < 1.0D0) THEN OVRLAP=.TRUE. RETURN ENDIF ! ! SW ATTRACTION ! PE=PE-1.0D0 ENDIF enddo ENDDO END SUBROUTINE SUM_ENR ! ! EVALUATE QUANTITIES OF INTEREST ! SUBROUTINE QUANTS(NSTEP,NBLOCK,NBL,NSITE, & XS,YS,ZS,PE, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G) implicit none integer :: NSTEP,NBLOCK,NBL,NSITE, & ISG,IB real(8) :: PE,STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & XCOM,YCOM,ZCOM real(8) :: XS(*),YS(*),ZS(*), & BSMPE(*),BSMR2(*),BSMR2G(*) ! ! CALCULATE INSTANTANEOUS VALUES ! STPPE=PE/DBLE(NSITE) ! ! END-TO-END DISTANCE ! STPR2=(XS(NSITE)-XS(1))**2 &
74
+(YS(NSITE)-YS(1))**2 +(ZS(NSITE)-ZS(1))**2
&
! ! COM (MASS OF EACH SEGMENT EQUAL TO 1) ! XCOM=0.0D0 YCOM=0.0D0 ZCOM=0.0D0 DO ISG=1,NSITE XCOM=XCOM+XS(ISG) YCOM=YCOM+YS(ISG) ZCOM=ZCOM+ZS(ISG) ENDDO XCOM=XCOM/DBLE(NSITE) YCOM=YCOM/DBLE(NSITE) ZCOM=ZCOM/DBLE(NSITE) ! ! RADIUS OF GYRATION ! STPR2G=0.0D0 DO ISG=1,NSITE STPR2G=STPR2G+((XS(ISG)-XCOM)**2) & +((YS(ISG)-YCOM)**2) & +((ZS(ISG)-ZCOM)**2) ENDDO STPR2G=STPR2G/DBLE(NSITE) ! ! UPDATE AVERAGES ! SUMPE=SUMPE+STPPE SUMR2=SUMR2+STPR2 SUMR2G=SUMR2G+STPR2G ! ! UPDATE FLUCTUATIONS ! SSQPE=SSQPE+STPPE**2 SSQR2=SSQR2+STPR2**2 SSQR2G=SSQR2G+STPR2G**2 ! ! UPDATE BLOCK AVERAGES ! IB=((NSTEP-1)/NBLOCK)+1 IF(IB > NBL) STOP "QUANTS: NBL TOO SMALL!" BSMPE(IB)=BSMPE(IB)+STPPE BSMR2(IB)=BSMR2(IB)+STPR2 BSMR2G(IB)=BSMR2G(IB)+STPR2G END SUBROUTINE QUANTS ! ! ZERO ACCUMULATORS ! SUBROUTINE ZERO(NBL,SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G) implicit none integer :: NBL real(8) :: SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G real(8) :: BSMPE(NBL),BSMR2(NBL),BSMR2G(NBL) ! ! AVERAGES ! SUMPE=0.0D0 SUMR2=0.0D0 SUMR2G=0.0D0 !
75
! FLUCTUATIONS ! SSQPE=0.0D0 SSQR2=0.0D0 SSQR2G=0.0D0 ! ! BLOCK AVERAGES ! BSMPE=0.0D0 BSMR2=0.0D0 BSMR2G=0.0D0 END SUBROUTINE ZERO ! ! OUTPUT ! SUBROUTINE OUTPUT(K,NBL,NBLOCK,KMAX,maxtemp,NSTEP, & NTRIAL,NACPTR,nacp_pt,ntri_pt, & TEMP,RLAMSW, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & BSMPE,BSMR2,BSMR2G, & TITLE, & myrank) implicit none integer :: K,NBL,NBLOCK,KMAX,maxtemp,NSTEP, & NTRIAL,NACPTR,nacp_pt,ntri_pt, & IBLOCK,IB,myrank,iunit real(8) :: TEMP,RLAMSW, & STPPE,STPR2,STPR2G, & SUMPE,SUMR2,SUMR2G, & SSQPE,SSQR2,SSQR2G, & NAVG,NAVGB,qpi,RAT_TRL, & rat_glo, & QRNR2,QRNR2G,QRNPI, & AVPE,AVR2,AVR2G, & FLCPE,FLCR2,FLCR2G, & BVPE,BVR2,BVR2G, & BFPE,BFR2,BFR2G real(8) :: BSMPE(NBL),BSMR2(NBL),BSMR2G(NBL) character(len=80) :: TITLE character(len=7) :: cnfile character(len=7),dimension(maxtemp) :: outfil ! ! File Names ! outfil(1) ="f01.out";outfil(2) ="f02.out";outfil(3) ="f03.out" outfil(4) ="f04.out";outfil(5) ="f05.out";outfil(6) ="f06.out" outfil(7) ="f07.out";outfil(8) ="f08.out";outfil(9) ="f09.out" outfil(10)="f10.out" ! ! LABEL THE OUTPUT ! IF(K == 0) THEN ! ! Define File Unit and Open Ooutput Files ! iunit=10+myrank*5 cnfile=outfil(myrank+1) open(unit=iunit,file=cnfile,status="unknown") ! WRITE(iunit,’(/,’’ **************************’’)’) WRITE(iunit,’(/,’’ SW CHAIN IN VACCUM ’’)’) WRITE(iunit,’(/,1X,A)’) TITLE WRITE(iunit,’(/,’’ **************************’’)’)
76
RETURN ENDIF ! ! INITIAL OUTPUT ! IF(K == 1) THEN write(iunit,fmt="(a)") "** SPECIFICATION **" write(iunit,fmt="(a,e13.5)") "kB.T: ",TEMP write(iunit,fmt="(a,e13.5)") "lambdaSW : ",RLAMSW write(iunit,fmt="(a,i4)") "k_max: ",KMAX RETURN ENDIF ! ! PERIODIC OUTPUT ! IF(K == 2) THEN ! ! print out every nprint steps ! ! STEP/RUNNING VALUES ! WRITE(iunit,’(/1X,’’ STEP/RUNNING VALUES’’)’) RAT_TRL=DBLE(NACPTR)/DBLE(NSTEP) if(ntri_pt == 0) then rat_glo=0.0d0 else rat_glo=DBLE(nacp_pt)/DBLE(ntri_pt) endif write(iunit,fmt="(a,a,a)") " # TRIAL "," ACP_TR "," ACP_PT " WRITE(iunit,’(1X,I10,6(3X,D8.2))’) NSTEP,RAT_TRL,rat_glo write(iunit,fmt="(a,a,a,a)") " PE ", & " R2 ", & " R2G ", & " Pi " QPI=STPR2/(6.0D0*STPR2G) WRITE(iunit,’(1X,5(D13.4))’) STPPE,STPR2,STPR2G,QPI QRNR2=SUMR2/DBLE(NSTEP) QRNR2G=SUMR2G/DBLE(NSTEP) QRNPI=QRNR2/(6.0D0*QRNR2G) WRITE(iunit,’(1X,5(D13.4))’) SUMPE/DBLE(NSTEP),QRNR2,QRNR2G,QRNPI RETURN ENDIF ! ! FINAL OUTPUT ! IF(K == 3) THEN ! ! CALCULATE THERMODYNAMIC AVERAGES ! NAVG=DBLE(NTRIAL) NAVGB=DBLE(NBLOCK) IBLOCK=NTRIAL/NBLOCK AVPE=SUMPE/NAVG AVR2=SUMR2/NAVG AVR2G=SUMR2G/NAVG ! ! CALCULATE ROOT MEAN SQUARE FLUCTUATIONS ! FLCPE=DSQRT(DABS(SSQPE/NAVG-AVPE**2)) FLCR2=DSQRT(DABS(SSQR2/NAVG-AVR2**2)) FLCR2G=DSQRT(DABS(SSQR2G/NAVG-AVR2G**2)) ! ! CALCULATE BLOCK AVERAGES ! BVPE=0.0D0
77
BVR2=0.0D0 BVR2G=0.0D0 DO IB=1,IBLOCK BSMPE(IB)=BSMPE(IB)/NAVGB BVPE=BVPE+BSMPE(IB) BSMR2(IB)=BSMR2(IB)/NAVGB BVR2=BVR2+BSMR2(IB) BSMR2G(IB)=BSMR2G(IB)/NAVGB BVR2G=BVR2G+BSMR2G(IB) ENDDO BVPE=BVPE/DBLE(IBLOCK) BVR2=BVR2/DBLE(IBLOCK) BVR2G=BVR2G/DBLE(IBLOCK) ! ! CALCULATE STANDARD DEVIATIONS ! BFPE=0.0D0 BFR2=0.0D0 BFR2G=0.0D0 IF(IBLOCK > 1) THEN DO IB=1,IBLOCK BFPE=BFPE+(BVPE-BSMPE(IB))**2 BFR2=BFR2+(BVR2-BSMR2(IB))**2 BFR2G=BFR2G+(BVR2G-BSMR2G(IB))**2 ENDDO BFPE=DSQRT(BFPE/DBLE(IBLOCK-1)) BFR2=DSQRT(BFR2/DBLE(IBLOCK-1)) BFR2G=DSQRT(BFR2G/DBLE(IBLOCK-1)) ENDIF WRITE(iunit,’(/1X,’’ RUN AVERAGES’’)’) write(iunit,fmt="(a,a,a,a)") " PE ", & " R2 ", & " R2G ", & " Pi " QPI=AVR2/(6.0D0*AVR2G) WRITE(iunit,’(1X,5(D13.4))’) AVPE,AVR2,AVR2G,QPI WRITE(iunit,’(/1X,’’ RMS FLUCTUATIONS’’)’) write(iunit,fmt="(a,a,a,a)") " PE ", & " R2 ", & " R2G ", & " Pi " WRITE(iunit,’(1X,5(D13.4))’) FLCPE,FLCR2,FLCR2G WRITE(iunit,’(/1X,’’ BLOCK AVERAGES’’)’) write(iunit,fmt="(a,a,a,a)") " PE ", & " R2 ", & " R2G ", & " Pi " QPI=BVR2/(6.0D0*BVR2G) WRITE(iunit,’(1X,5(D13.4))’) BVPE,BVR2,BVR2G,QPI WRITE(iunit,’(/1X,’’ STANDARD DEVIATIONS’’)’) write(iunit,fmt="(a,a,a,a)") " PE ", & " R2 ", & " R2G ", & " Pi " WRITE(iunit,’(1X,5(D13.4))’) BFPE,BFR2,BFR2G ! close(unit=iunit) ENDIF END SUBROUTINE OUTPUT ! ! RANDOM VECTOR ON THE SURFACE OF 3D SPHERE ! SUBROUTINE SPH3E(EX,EY,EZ,idum) implicit none integer :: idum
78
real(8) :: random real(8) :: EX,EY,EZ, & DZETA1,DZETA2,DZETSQ,QPAR ! do DZETA1=2.0D0*random(idum)-1.0D0 DZETA2=2.0D0*random(idum)-1.0D0 DZETSQ=DZETA1*DZETA1+DZETA2*DZETA2 IF(DZETSQ > 1.0D0) cycle QPAR=DSQRT(1.0D0-DZETSQ) EX=2.0D0*DZETA1*QPAR EY=2.0D0*DZETA2*QPAR EZ=1.0D0-2.0D0*DZETSQ exit enddo END SUBROUTINE SPH3E ! ! SELECT A NEW ORIENTATION WITH PROBABILITY = w(k) / Sum w(k) ! SUBROUTINE SELECT(W,SUMW,idum,N) implicit none integer :: idum,N real(8) :: random real(8) :: SUMW,WS,CUMW real(8) :: W(*) ! WS=random(idum)*SUMW CUMW=W(1) N=1 ! DO WHILE(CUMW < WS) N=N+1 CUMW=CUMW+W(N) ENDDO END SUBROUTINE SELECT ! ! RANDOM GENERATOR (Numerical Recipes in Fortran90) ! FUNCTION random(idum) ! ! "Minimal" random number generator of Park and Miller combined with a Marsaglia shift ! sequence. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint ! values). This fully portable, scalar generator has the "traditional" (not Fortran 90) calling ! sequence with a random deviate as the returned function value: call with idum a negative ! integer to initialize; thereafter, do not alter idum except to reinitialize. The period of this ! generator is about 3.1 x 10^18. ! IMPLICIT NONE ! ! Declaration of Variables ! INTEGER,PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B),INTENT(INOUT) :: idum REAL :: ran INTEGER(K4B),PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836 REAL,SAVE :: am INTEGER(K4B),SAVE :: ix=-1,iy=-1,k ! real(8) :: random ! ! Initialize ! if((idum <= 0).or.(iy < 0)) then am=nearest(1.0,-1.0)/IM iy=ior(ieor(888889999,abs(idum)),1)
79
ix=ieor(777755555,abs(idum)) ! ! Set idum Positive ! idum=abs(idum)+1 endif ! ! Marsaglia Shift Sequence with Period 2^32 - 1 ! ix=ieor(ix,ishft(ix,13)) ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) ! ! Park-Miller Sequence by Schrage’s Method, Period 2^31 - 2 ! k=iy/IQ iy=IA*(iy-k*IQ)-IR*k if(iy < 0) iy=iy+IM ! ! Combine the Two Generators with Masking to Ensure Nonzero Value ! ran=am*ior(iand(IM,ieor(ix,iy)),1) ! random=dble(ran) END FUNCTION random
Obr. 18 ukazuje příklad výpočtu R2 , Rg2 a < pro polymer tvořený 10 segmenty s λ = 1.5; TC−G = 2.85. Pro jednoduchost se uvažovalo σ = 1 a = 1. Rosenblutovo MC začalo z náhodné konformace, 100 000 MC kroků bylo použito pro „ekvilibraci” a následných 5 000 000 MC kroků bylo použito pro výpočet průměrů < . >; k = 10.14
14
Pro kontrolu správnosti programu uveďme konkrétní hodnoty pro kB T = 1: R2 = 6.33020 , Rg2 = 1.4602 a < = 0.7224.
80
Obrázek 18: R2 , Rg2 a < pro polymer tvořený 10 segmenty s λ = 1.5, σ = 1 a = 1. 81
3.3 3.3.1
Paralelní molekulární dynamika Úvod
Molekulární dynamika (MD) slouží k řešení pohybu N molekul, jejichž interakce je popsána potenciálem u. Uvažujme pro jednoduchost systém částic reagujích LennardovýmJonesovým (LJ) potenciálem
σ uLJ (rij ) = 4ε rij
!12
σ − rij
!6
(37)
kde rij = |ri − rj | je vzdálenost mezi molekulami i a j, ε a σ jsou resp. energetický a délkový LJ parametr. V MD programech se více než 90 % výpočtového času spotřebuje na výpočet celkové potenciální energie systému U , viriálu W a sil fij . Ty se spočtou pro LJ tekutinu následujícím způsobem: U = Uc + ULRC =
N −1 X
N X
"
u
LJ
i=1 j=i+1
8 1 σ (rij ) + πN ερσ 3 3 3 rc
9
σ − rc
3 #
(38)
kde rc je poloměr ořezání a ρ je hustota. W = Wc + WLRC " 3 # N −1 X 16 1 NX σ σ 9 LJ 3 2 w (rij ) + πN ερσ = − − 3 i=1 j=i+1 3 3 rc rc
(39)
kde wLJ = r duLJ /dr . "
fij
3.3.2
#
1 duLJ (rij ) = − rij rij drij = −fji
(40) (41)
Paralelizace
Paralelizací dvojnásobného cyklu do i=1,N-1 ! Vnejsi cyklus do j=i+1,N ! Vnitrni cyklus ... enddo enddo
ve výpočtu Uc , Wc a fij lze docílit značného urychlení běhu MD programů. Paralelizaci dvojnásobného cyklu lze provést následujícím způsobem: 82
1. Každý proces spočte jen část vnějšího cyklu t.j. spočte jen část z Uc , Wc a fij . To lze provést následujícím způsobem:
do i=1+myrank,N-1,nprocs ! Vnejsi cyklus do j=i+1,N ! Vnitrni cyklus ... enddo enddo
2. Použitím MPI příkazu MPI ALLREDUCE sečteme části Uc , Wc a fij z jednotlivých procesů a výsledné Uc , Wc a fij se „uloží” na všechny procesy. Program pro paralelní MD LJ tekutiny může vypadat následovně: ! ! Parallel NVT MD of Lennard-Jones Fluids ! ! uLJ=4*eps*[(sig/r)^12-(sig/r)^6] ! program ParNVTmdLJ implicit none include ’mpif.h’ ! preprocessor directive ! ! System Data ! integer,parameter :: nmax=5000 ! Max # of Molecules real(8),parameter :: dt=0.005d0 ! Time Step real(8) :: rcut=1.0d0 ! Cut-Off Radius integer,parameter :: nom=500 ! # of Molecules real(8),parameter :: temp=2.0d0 ! Reduced Temperature real(8),parameter :: dens=0.4d0 ! Reduced Number Density ! ! Declaration of Variables ! integer :: nrun,nprint,nstep real(8) :: vol,box,boxinv,rcutsq,time, & U,W,U_LRC,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P real(8),dimension(nmax) :: x,y,z, & vx,vy,vz, & fx,fy,fz character(len=3) :: yesno character(len=30) :: cnfile integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,ip ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! ! get my process rank !
83
call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! Input ! if(myrank == 0) then print*,"# of Time Steps:" read*,nrun print*,"# of Time Steps for Printing:" read*,nprint print*,"Read Configurational File (yes/no):" read*,yesno print*,"Name of Configurational File:" read*,cnfile endif ! ! broadcast input ! call MPI_BCAST(nrun,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nprint,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(yesno,3,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(cnfile,30,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) ! ! Simulation Set-Up Values ! call values(nom,rcut,dens,vol,box,boxinv,rcutsq,U_LRC,W_LRC) ! ! Initial Output ! if(myrank == 0) then call output(1,nstep,nrun,nom,rcut,temp,dens,vol,box, & time,u_ac,P_ac,sum_u,sum_P,ssq_u,ssq_P) endif ! ! Initialization ! if((yesno == ’yes’).or.(yesno == ’YES’)) then ! ! Read from File; One Processor at a Time ! do ip=0,nprocs-1 if(myrank == ip) then call readcn(nmax,x,y,z,vx,vy,vz,cnfile) endif call MPI_BARRIER(MPI_COMM_WORLD,ierr) enddo else ! ! Positions from fcc Lattice and ! Velocity from Maxwell-Boltzmann’s Distribution ! call init_conf(nom,temp,box,x,y,z,vx,vy,vz) endif ! ! Zeroth Thermo Accumulators ! call zero(sum_u,sum_P,ssq_u,ssq_P) ! ! MD cyclus ! time=0.0d0 do nstep=1,nrun ! ! Energy, Virial and Forces ! call force(nom,box,boxinv,rcutsq,U,W, &
84
x,y,z,fx,fy,fz, myrank,nprocs)
&
! ! Leap-Frog Algorithm ! call Leap_Frog(nom,dt,temp,box,boxinv, & x,y,z,vx,vy,vz,fx,fy,fz) ! ! Termodynamic Quantities ! call thermo(nom,temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) ! ! Periodic Output ! time=time+dt if(myrank == 0) then if(mod(nstep,nprint) == 0) then call output(2,nstep,nrun,nom,rcut,temp,dens,vol,box, & time,u_ac,P_ac,sum_u,sum_P,ssq_u,ssq_P) endif endif enddo ! ! Final Write to File ! if(myrank == 0) then call writcn(nmax,x,y,z,vx,vy,vz,cnfile) ! ! Final Output ! call output(3,nstep,nrun,nom,rcut,temp,dens,vol,box, & time,u_ac,P_ac,sum_u,sum_P,ssq_u,ssq_P) endif ! ! shut down MPI ! call MPI_FINALIZE(ierr) ! stop "ParNVTmdLJ: End of Calcs!" end program ParNVTmdLJ ! ! Simulation Set-Up Values ! subroutine values(nom,rcut,dens,vol,box,boxinv,rcutsq,U_LRC,W_LRC) implicit none ! ! Declaration of Variables ! integer :: nom real(8) :: rcut,rcutsq,dens,vol,box,boxinv, & pi,ircut3,ircut9,U_LRC,W_LRC ! ! Volume, Length of Box and Cut-Off Radius ! vol=dble(nom)/dens box=vol**(1.0/3.0) boxinv=1.0d0/box if(rcut > 1.0d0) rcut=1.0d0 rcut=rcut*(box/2.0d0) rcutsq=rcut**2 ! ! LJ LRC ! pi=dacos(-1.0d0)
85
ircut3=1.0d0/(rcutsq*rcut) ircut9=ircut3*ircut3*ircut3 ! U_LRC=(8.0d0/3.0d0)*pi*dble(nom)*dens*((1.0d0/3.0d0)*ircut9-ircut3) W_LRC=(16.0d0/3.0d0)*pi*dble(nom)*dens*((2.0d0/3.0d0)*ircut9-ircut3) end subroutine values ! ! Positions from fcc Lattice and ! Velocity from Maxwell-Boltzmann’s Distribution ! subroutine init_conf(nom,temp,box,x,y,z,vx,vy,vz) implicit none ! ! Declaration of Variables ! integer :: nom,ncc,m,ix,iy,iz,iref,i,iseed real(8) :: gauss real(8) :: temp,box, & x(*),y(*),z(*), & vx(*),vy(*),vz(*), & cell,cell2,rtemp, & sumvx,sumvy,sumvz,sumv2, & temp_conf,fs ! ! Check # of Molecules ! ncc=(nom/4)**(1.0/3.0)+0.1d0 if(4*ncc**3 /= nom) stop "init_conf: Wrong # of Molecules!" ! ! Positions of Molecules, fcc Lattice ! cell=1.0d0/dble(ncc) cell2=0.5d0*cell ! x(1)=0.0d0 y(1)=0.0d0 z(1)=0.0d0 ! x(2)=cell2 y(2)=cell2 z(2)=0.0d0 ! x(3)=0.0d0 y(3)=cell2 z(3)=cell2 ! x(4)=cell2 y(4)=0.0d0 z(4)=cell2 ! ! Construct Lattice from Unit Cell ! m=0 do iz=1,ncc do iy=1,ncc do ix=1,ncc do iref=1,4 x(iref+m)=x(iref)+cell*dble(ix-1) y(iref+m)=y(iref)+cell*dble(iy-1) z(iref+m)=z(iref)+cell*dble(iz-1) enddo m=m+4 enddo enddo enddo
86
! ! Shift Positions ! do i=1,nom x(i)=(x(i)-0.5d0)*box y(i)=(y(i)-0.5d0)*box z(i)=(z(i)-0.5d0)*box enddo ! ! Velocity from Maxwell-Boltzmann’s Distribution, ! Total Momentum and Kinetic Energy ! iseed=425001 rtemp=dsqrt(temp) sumvx=0.0d0 sumvy=0.0d0 sumvz=0.0d0 sumv2=0.0d0 do i=1,nom vx(i)=rtemp*gauss(iseed) vy(i)=rtemp*gauss(iseed) vz(i)=rtemp*gauss(iseed) sumvx=sumvx+vx(i) sumvy=sumvy+vy(i) sumvz=sumvz+vz(i) sumv2=sumv2+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i) enddo sumvx=sumvx/dble(nom) sumvy=sumvy/dble(nom) sumvz=sumvz/dble(nom) ! ! Actual Temperature ! temp_conf=sumv2/dble(3*nom-4) ! ! Scaling Factor ! fs=dsqrt(temp/temp_conf) ! ! Remove Total Momentum and Scale Velocities ! do i=1,nom vx(i)=fs*(vx(i)-sumvx) vy(i)=fs*(vy(i)-sumvy) vz(i)=fs*(vz(i)-sumvz) enddo end subroutine init_conf ! ! Calculate Energy, Virial and Forces ! subroutine force(nom,box,boxinv,rcutsq,U,W, & x,y,z,fx,fy,fz, & myrank,nprocs) implicit none include ’mpif.h’ ! preprocessor directive ! ! Declaration of Variables ! integer :: nom,i,j real(8) :: box,boxinv,rcutsq,U,W, & xij,yij,zij,fxij,fyij,fzij, & rijsq,sr2,sr6,uij,wij,fij real(8) :: x(*),y(*),z(*), & fx(*),fy(*),fz(*) integer :: nprocs, & ! # of processes
87
myrank, & ! my process rank ierr real(8) :: qU,qW real(8) :: qfx(nom),qfy(nom),qfz(nom) ! ! Zeroth Energy, Virial and Forces ! qU=0.0d0 qW=0.0d0 qfx=0.0d0 qfy=0.0d0 qfz=0.0d0 ! ! Outer Cycle ! do i=1+myrank,nom-1,nprocs ! ! Inner Cycle ! do j=i+1,nom ! ! Distance Between Molecules ! xij=x(i)-x(j) yij=y(i)-y(j) zij=z(i)-z(j) ! ! Periodic Boundary Conditions ! xij=xij-dnint(xij*boxinv)*box yij=yij-dnint(yij*boxinv)*box zij=zij-dnint(zij*boxinv)*box ! ! Cut-Off ! rijsq=xij*xij+yij*yij+zij*zij if(rijsq < rcutsq) then ! ! Auxiliary Variables ! sr2=1.0d0/rijsq sr6=sr2*sr2*sr2 ! ! Potential Energy ! uij=4.0d0*sr6*(sr6-1.0d0) qU=qU+uij ! ! Virial ! wij=(-1.0d0)*24.0d0*sr6*(2.0d0*sr6-1.0d0) qW=qW-(wij/3.0d0) ! ! Forces ! fij=(-1.0d0)*(wij/rijsq) ! - w(r_ij)/r_ij^2 fxij=fij*xij fyij=fij*yij fzij=fij*zij qfx(i)=qfx(i)+fxij qfy(i)=qfy(i)+fyij qfz(i)=qfz(i)+fzij qfx(j)=qfx(j)-fxij qfy(j)=qfy(j)-fyij qfz(j)=qfz(j)-fzij
88
endif enddo enddo ! ! Add Results Together ! call MPI_ALLREDUCE(qfx,fx,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfy,fy,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfz,fz,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qU,U,1,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qW,W,1,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) end subroutine force ! ! Leap-Frog Algorithm ! subroutine Leap_Frog(nom,dt,temp,box,boxinv, & x,y,z,vx,vy,vz,fx,fy,fz) implicit none ! ! Declaration of Variables ! integer :: nom,i real(8) :: dt,temp,box,boxinv,sum_vx,sum_vy,sum_vz, K,temp_ac,beta real(8) :: x(*),y(*),z(*), & vx(*),vy(*),vz(*), & fx(*),fy(*),fz(*) ! ! Zeroth Accumulators ! K=0.0d0 sum_vx=0.0d0 sum_vy=0.0d0 sum_vz=0.0d0 do i=1,nom ! ! v_u(t)=v(t-dt/2)+(dt/2)*[f(t)/m] ! vx(i)=vx(i)+(dt/2.0d0)*fx(i) vy(i)=vy(i)+(dt/2.0d0)*fy(i) vz(i)=vz(i)+(dt/2.0d0)*fz(i) ! ! Kinetic Energy ! K=K+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i) ! ! Total Momentum ! sum_vx=sum_vx+vx(i) sum_vy=sum_vy+vy(i) sum_vz=sum_vz+vz(i) enddo sum_vx=sum_vx/dble(nom) sum_vy=sum_vy/dble(nom) sum_vz=sum_vz/dble(nom) ! ! Actual Temperature ! temp_ac=K/dble(3*nom-4) !
&
89
! Scaling Factor ! beta=dsqrt(temp/temp_ac) ! ! v(t)=beta*v_u(t) ! v(t+dt/2)=(2-1/beta)*v(t)+(dt/2)*[f(t)/m] ! r(t+dt)=r(t)+dt*v(t+dt/2) ! K=0.0d0 do i=1,nom ! ! v(t)=beta*v_u(t) ! vx(i)=beta*(vx(i)-sum_vx) vy(i)=beta*(vy(i)-sum_vy) vz(i)=beta*(vz(i)-sum_vz) ! ! Kinetic Energy ! K=K+vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i) ! ! v(t+dt/2)=(2-1/beta)*v(t)+(dt/2)*[f(t)/m] ! vx(i)=(2.0d0-(1.0d0/beta))*vx(i)+(dt/2.0d0)*fx(i) vy(i)=(2.0d0-(1.0d0/beta))*vy(i)+(dt/2.0d0)*fy(i) vz(i)=(2.0d0-(1.0d0/beta))*vz(i)+(dt/2.0d0)*fz(i) ! ! r(t+dt)=r(t)+dt*v(t+dt/2) ! x(i)=x(i)+vx(i)*dt y(i)=y(i)+vy(i)*dt z(i)=z(i)+vz(i)*dt ! ! Periodic Boundary Conditions ! x(i)=x(i)-dnint(x(i)*boxinv)*box y(i)=y(i)-dnint(y(i)*boxinv)*box z(i)=z(i)-dnint(z(i)*boxinv)*box enddo end subroutine Leap_Frog ! ! Zeroth Thermo Accumulators ! subroutine zero(sum_u,sum_P,ssq_u,ssq_P) implicit none ! ! Declaration of Variables ! real(8) :: sum_u,sum_P,ssq_u,ssq_P ! ! Zeroth sum Accumulators ! sum_u=0.0d0 sum_P=0.0d0 ! ! Zeroth ssq Accumulators ! ssq_u=0.0d0 ssq_P=0.0d0 end subroutine zero ! ! Termodynamic Quantities ! subroutine thermo(nom,temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, sum_u,sum_P,ssq_u,ssq_P)
90
&
implicit none ! ! Declaration of Variables ! integer :: nom real(8) :: temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P ! ! Instantaneous Thermodynamic Quantities ! u_ac=(U+U_LRC)/dble(nom) ! energy per particle P_ac=dens*temp+(W+W_LRC)/vol ! pressure ! ! sum Accumulators ! sum_u=sum_u+u_ac sum_P=sum_P+P_ac ! ! ssq Accumulators ! ssq_u=ssq_u+u_ac**2 ssq_P=ssq_P+P_ac**2 end subroutine thermo ! ! Output ! subroutine output(kk,nstep,nrun,nom,rcut,temp,dens,vol,box, & time,u_ac,P_ac,sum_u,sum_P,ssq_u,ssq_P) implicit none ! ! Declaration of Variables ! integer :: kk,nstep,nrun,nom real(8) :: rcut,temp,dens,vol,box,time,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P, & avr_u,avr_P,flc_u,flc_P if(kk == 1) then ! ! Initial Output ! write(*,*) " *** Initial Output ***" write(*,fmt="(a,e13.5)") "Temperature: ",temp write(*,fmt="(a,e13.5)") "Density: ",dens write(*,*) write(*,fmt="(a,i4)") "# of Molecules: ",nom write(*,fmt="(a,e13.5)") "Cut-Off Distance: ",rcut write(*,fmt="(a,e13.5)") "Volume of Box: ",vol write(*,fmt="(a,e13.5)") "Length of Box: ",box return endif if(kk == 2) then ! ! Periodic Output ! write(*,*) write(*,fmt="(a,i8,a,e13.5)") " # of time steps:",nstep, & " time:",time write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") u_ac,P_ac write(*,fmt="(8(e13.5))") sum_u/dble(nstep),sum_P/dble(nstep) return endif if(kk == 3) then ! ! Final Output
91
! ! Ensemble Averages ! avr_u=sum_u/dble(nrun) avr_P=sum_P/dble(nrun) ! ! Ensemble Fluctuations ! flc_u=dsqrt(dabs((ssq_u/dble(nrun))-avr_u**2)) flc_P=dsqrt(dabs((ssq_P/dble(nrun))-avr_P**2)) write(*,*) write(*,*) " *** Final Output ***" write(*,*) " Ensemble Averages" write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") avr_u,avr_P write(*,*) " Ensemble Fluctuations" write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") flc_u,flc_P return endif end subroutine output ! ! Read from File ! subroutine readcn(nmax,x,y,z,vx,vy,vz,cnfile) implicit none ! ! Declaration of Variables ! integer :: nmax real(8) :: x(nmax),y(nmax),z(nmax), & vx(nmax),vy(nmax),vz(nmax) integer,parameter :: cnunit=10 character(len=*) :: cnfile ! ! Unformatted File ! open(unit=cnunit,file=cnfile,status="old",form="unformatted") read(cnunit) x,y,z read(cnunit) vx,vy,vz close(unit=cnunit) end subroutine readcn ! ! Write to File ! subroutine writcn(nmax,x,y,z,vx,vy,vz,cnfile) implicit none ! ! Declaration of Variables ! integer :: nmax real(8) :: x(nmax),y(nmax),z(nmax), & vx(nmax),vy(nmax),vz(nmax) integer,parameter :: cnunit=10 character(len=*) :: cnfile ! ! Unformatted File ! open(unit=cnunit,file=cnfile,status="unknown",form="unformatted") write(cnunit) x,y,z write(cnunit) vx,vy,vz close(unit=cnunit) end subroutine writcn ! ! RANDOM VARIATE FROM THE STANDARD NORMAL DISTRIBUTION.
92
! THE DISTRIBUTION IS GAUSSIAN WITH ZERO MEAN AND UNIT VARIANCE. ! KNUTH D, THE ART OF COMPUTER PROGRAMMING, ADDISON-WESLEY, 1978 ! function gauss(iseed) implicit none ! ! Declaration of Variables ! integer :: iseed,I real(8) :: gauss real(8),parameter :: A1=3.949846138d0,A3=0.252408784d0, & A5=0.076542912d0,A7=0.008355968d0, & A9=0.029899776d0 real(8) :: SUM,R,R2 ! SUM=0.0d0 do I=1,12 SUM=SUM+dble(ran(iseed)) enddo ! R=(SUM-6.0d0)/4.0d0 R2=R*R gauss=((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R end function gauss
3.4 3.4.1
Paralelní Monte Carlo Úvod
Monte Carlo metoda podobně jako MD slouží ke vzorkování fázového prostoru N molekul, které interagují pomocí potenciálu u. Uvažujme opět pro jednoduchost systém LJ částic. Ve standardní MC metodě založené na Metropolisově algoritmu pohybujeme postupně vždy jednou částicí. Pohyb částice i přijímáme s pravděpodobností, která závisí na teplotě a rozdílu potenciální energie na částici i v původní a nové poloze. Potenciální energii na LJ částici i spočteme jako Ui =
N X j=1,i6=j
"
LJ
u
1 σ 8 (rij ) + περσ 3 3 3 rc
9
σ − rc
3 #
(42)
V MC programech podobně jako v MD programech se více než 90% výpočtového času spotřebuje na výpočet Ui . Výpočet Ui , který se realizuje v jednoduchém cyklu do j=1,N if(i /= j) then ’vypocet Ui’ endif enddo
lze paralelizovat např. následujícím způsobem:
93
irank=-1 ! do j=1,N irank=irank+1 if(irank == nprocs) irank=0 if(myrank == irank) then if(i /= j) then ’vypocet Ui’ endif endif enddo ! call MPI_ALLREDUCE(...)
Vzhledem k tomu, že ve standardním MC paralelizujeme jednoduchý cyklus, je efekt takové paralelizace nevýznamný.
3.4.2
Hybridní Monte Carlo metoda
Vedle standardní MC metody existuje tzv. hybridní MC metoda, která kombinuje prvky MD a Metropolisova MC. Pohyb částic z původní konfigurace o do nové konfigurace n se v hybridním MC realizuje pomocí několika MD kroků a nová konfigurace se přijme obdobně jako v Metropolisově MC. Jeden cyklus hybridního MC pro systém N částic v objemu V a při teplotě T se skládá z následujících kroků. 1. Vycházíme z původní konfigurace, kde známe polohy částic {roi }N i=1 . 2. Vygenerujeme rychlosti částic {vio }N i=1 pomocí Maxwellova-Boltzmannova rozdělení odpovídající teplotě T a následně spočteme celkovou energii systému (součet kinetické a potenciální energie) jako Ho =
N 1X mi vio · vio + U (roi ) 2 i=1
(43)
3. Provedeme pohyb všech částic systému pomocí nMD MD kroků s krokem integrace ∆t. Pohyb částic je nutno vzhledem k podmínce mikroskopické reverzibility realizovat pomocí reverzibilního integračního algoritmu, např. pomocí „leap-frog” algoritmu: ∆t ) 2
=
ri (to + ∆t)
=
3∆t ) 2
=
ri (to + 2∆t)
=
vi (to +
vi (to +
vio +
∆t fio 2 mi
∆t ) 2 ∆t fi (to + ∆t) vi (to + ) + ∆t 2 mi 3∆t ri (to + ∆t) + ∆t vi (to + ) 2 roi + ∆t vi (to +
94
(44)
. . . ∆t vi (tn − ) 2
=
rni
=
vin
=
3∆t fi (tn − ∆t) ) + ∆t 2 mi ∆t ri (tn − ∆t) + ∆t vi (tn − ) 2 ∆t ∆t fin vi (tn − )+ 2 2 mi vi (tn −
4. Spočteme celkovou energii systému v nové konfiguraci jako Hn =
N 1X mi vin · vin + U (rni ) 2 i=1
(45)
a novou konfiguraci přijmeme, pokud platí exp [−β (Hn − Ho )] ≥ ζ
3.4.3
(46)
Paralelizace
Hybridní MC podobně jako MD vyžaduje výpočet sil a potenciální energie, které se realizují v dvojnásobném cyklu do i=1,N-1 ! Vnejsi cyklus do j=i+1,N ! Vnitrni cyklus ... enddo enddo
Paralelizací dvojnásobného cyklu lze docílit značného urychlení běhu hybridních MC programů. Paralelizaci dvojnásobného cyklu lze provést stejným způsobem jako v případě paralelní MD. Uveďme si pro zajímavost alternativní způsob paralelizace dvojnásobného cyklu: irank=-1 ! do i=1,N-1 ! Vnejsi cyklus do j=i+1,N ! Vnitrni cyklus irank=irank+1 if(irank == nprocs) irank=0 if(myrank == irank) then ... endif enddo enddo ! call MPI_ALLREDUCE(...)
95
Program pro paralelní hybridní MC LJ tekutiny může vypadat následovně:
! ! Parallel Hybrid Monte Carlo of Lennard-Jones Fluids ! ! uLJ=4*eps*[(sig/r)^12-(sig/r)^6] ! program ParHybridMC implicit none include ’mpif.h’ ! preprocessor directive ! ! System Data ! integer,parameter :: nmax=500 ! Max # of Molecules integer :: n_MD=10 ! # of MD Steps per MC Cycle real(8),parameter :: dt=0.005d0 ! Time Step real(8) :: rcut=1.0d0 ! Cut-Off Radius ! integer,parameter :: nom=500 ! # of Molecules real(8),parameter :: temp=2.0d0 ! Reduced Temperature real(8),parameter :: dens=0.4d0 ! Reduced Number Density ! integer :: idum=-471 ! Seed for Random Generator ! ! Declaration of Variables ! integer :: nrun,nprint,nadj,ncycl,nstep,nacpt real(8) :: random real(8) :: vol,box,boxinv,rcutsq,U,W,U_LRC,W_LRC,K,Kn,Ho,Un,Wn,Hn,u_ac,P_ac, delH,beta_delH,rat_acpt,sum_u,sum_P,ssq_u,ssq_P real(8),dimension(nmax) :: x,y,z, & xn,yn,zn, & vx,vy,vz, & fx,fy,fz, & fxn,fyn,fzn character(len=3) :: yesno character(len=30) :: cnfile integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,ip ! ! start up MPI ! call MPI_INIT(ierr) ! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) ! ! get my process rank ! call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) ! ! Input ! if(myrank == 0) then print*,"# of MC Cycles:" read*,nrun print*,"# of Time Steps for Printing:" read*,nprint print*,"Read Configurational File (yes/no):" read*,yesno print*,"Name of Configurational File:"
96
&
read*,cnfile print*,"# of Cycles for Adjustment:" read*,nadj endif ! ! broadcast input ! call MPI_BCAST(nrun,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nprint,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(yesno,3,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(cnfile,30,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(nadj,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ! ! Simulation Set-Up Values ! call values(nom,rcut,dens,vol,box,boxinv,rcutsq,U_LRC,W_LRC) ! ! Initial Output ! if(myrank == 0) then call output(1,ncycl,nrun,nacpt,nom,rcut,temp,dens,vol,box,u_ac,P_ac, sum_u,sum_P,ssq_u,ssq_P) endif ! ! Zeroth Thermo Accumulators ! call zero(sum_u,sum_P,ssq_u,ssq_P) ! ! Initialization ! if((yesno == ’yes’).or.(yesno == ’YES’)) then ! ! Read from File; One Processor at a Time ! do ip=0,nprocs-1 if(myrank == ip) then call readcn(nmax,n_MD,x,y,z,vx,vy,vz,cnfile) endif call MPI_BARRIER(MPI_COMM_WORLD,ierr) enddo else ! ! Initial Positions from fcc Lattice ! call fcc_lattice(nom,box,x,y,z) endif ! ! Initial Energy, Virial and Forces ! call force(0,n_MD,nom,box,boxinv,rcutsq,U,W,x,y,z,fx,fy,fz, & myrank,nprocs) ! ! MC Cyclus ! nacpt=0 ! do ncycl=1,nrun ! ! Velocity from Maxwell-Boltzmann’s Distribution ! call MB_velocity(nom,temp,vx,vy,vz,idum) ! ! Kinetic Energy ! call kinetic_energy(nom,vx,vy,vz,K)
97
&
! ! Total Energy in Old Configuration ! Ho=K+U+U_LRC ! ! MD Steps ! xn=x yn=y zn=z ! fxn=fx fyn=fy fzn=fz ! do nstep=1,n_MD+1 ! ! Energy, Virial and Forces ! if(nstep > 1) call force(nstep,n_MD,nom,box,boxinv,rcutsq,Un,Wn,xn,yn,zn,fxn,fyn,fzn, myrank,nprocs) ! ! Leap-Frog Algorithm ! call Leap_Frog(nom,nstep,n_MD,dt,box,boxinv,xn,yn,zn,vx,vy,vz,fxn,fyn,fzn) enddo ! ! Kinetic Energy ! call kinetic_energy(nom,vx,vy,vz,Kn) ! ! Total Energy in New Configuration ! Hn=Kn+Un+U_LRC ! ! Metropolis Test ! delH=Hn-Ho beta_delH=delH/temp ! if(beta_delH < 75.0d0) then if(beta_delH < 0.0d0) then nacpt=nacpt+1 ! ! Bookkeeping ! U=Un W=Wn ! x=xn y=yn z=zn ! fx=fxn fy=fyn fz=fzn else if(dexp((-1.0d0)*beta_delH) > random(idum)) then nacpt=nacpt+1 ! ! Bookkeeping ! U=Un W=Wn ! x=xn
98
&
y=yn z=zn ! fx=fxn fy=fyn fz=fzn endif endif ! ! Termodynamic Quantities ! call thermo(nom,temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) ! ! Periodic Output ! if(myrank == 0) then if(mod(ncycl,nprint) == 0) then call output(2,ncycl,nrun,nacpt,nom,rcut,temp,dens,vol,box,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) endif endif ! ! Adjust # of MD Steps ! if(mod(ncycl,nadj) == 0) then rat_acpt=dble(nacpt)/dble(ncycl) ! if(rat_acpt > 0.3d0) then n_MD=n_MD+1 else if(n_MD > 1) n_MD=n_MD-1 endif ! nacpt=0 endif enddo ! ! Final Write to File ! if(myrank == 0) then call writcn(nmax,n_MD,x,y,z,vx,vy,vz,cnfile) ! ! Final Output ! call output(3,ncycl,nrun,nacpt,nom,rcut,temp,dens,vol,box,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) endif ! stop " ParHybridMC: End of Calcs!" end program ParHybridMC ! ! Simulation Set-Up Values ! subroutine values(nom,rcut,dens,vol,box,boxinv,rcutsq,U_LRC,W_LRC) implicit none ! ! Declaration of Variables ! integer :: nom real(8) :: rcut,rcutsq,dens,vol,box,boxinv, & pi,ircut3,ircut9,U_LRC,W_LRC ! ! Volume, Length of Box and Cut-Off Radius !
99
vol=dble(nom)/dens box=vol**(1.0/3.0) boxinv=1.0d0/box if(rcut > 1.0d0) rcut=1.0d0 rcut=rcut*(box/2.0d0) rcutsq=rcut**2.0 ! ! LJ LRC ! pi=dacos(-1.0d0) ircut3=1.0d0/(rcutsq*rcut) ircut9=ircut3*ircut3*ircut3 ! U_LRC=(8.0d0/3.0d0)*pi*dble(nom)*dens*((1.0d0/3.0d0)*ircut9-ircut3) W_LRC=(16.0d0/3.0d0)*pi*dble(nom)*dens*((2.0d0/3.0d0)*ircut9-ircut3) end subroutine values ! ! Positions from fcc Lattice ! subroutine fcc_lattice(nom,box,x,y,z) implicit none ! ! Declaration of Variables ! integer :: nom,ncc,m,ix,iy,iz,iref,i real(8) :: box,cell,cell2 real(8) :: x(*),y(*),z(*) ! ! Check # of Molecules ! ncc=(nom/4)**(1.0/3.0)+0.1d0 if(4*ncc**3 /= nom) stop " fcc_lattice: Wrong # of Molecules!" ! ! Positions of Molecules, fcc Lattice ! cell=1.0d0/dble(ncc) cell2=0.5d0*cell ! x(1)=0.0d0 y(1)=0.0d0 z(1)=0.0d0 ! x(2)=cell2 y(2)=cell2 z(2)=0.0d0 ! x(3)=0.0d0 y(3)=cell2 z(3)=cell2 ! x(4)=cell2 y(4)=0.0d0 z(4)=cell2 ! ! Construct Lattice from Unit Cell ! m=0 do iz=1,ncc do iy=1,ncc do ix=1,ncc do iref=1,4 x(iref+m)=x(iref)+cell*dble(ix-1) y(iref+m)=y(iref)+cell*dble(iy-1) z(iref+m)=z(iref)+cell*dble(iz-1) enddo
100
m=m+4 enddo enddo enddo ! ! Shift Positions ! do i=1,nom x(i)=(x(i)-0.5d0)*box y(i)=(y(i)-0.5d0)*box z(i)=(z(i)-0.5d0)*box enddo end subroutine fcc_lattice ! ! Velocity from Maxwell-Boltzmann’s Distribution ! subroutine MB_velocity(nom,temp,vx,vy,vz,idum) implicit none ! ! Declaration of Variables ! integer :: nom,i,idum real(8) :: gauss real(8) :: temp,rtemp real(8) :: vx(*),vy(*),vz(*) ! ! Maxwell-Boltzmann’s Distribution ! rtemp=dsqrt(temp) ! do i=1,nom vx(i)=rtemp*gauss(idum) vy(i)=rtemp*gauss(idum) vz(i)=rtemp*gauss(idum) enddo end subroutine MB_velocity ! ! Kinetic Energy ! subroutine kinetic_energy(nom,vx,vy,vz,K) implicit none ! ! Declaration of Variables ! integer :: nom,i real(8) :: K real(8) :: vx(*),vy(*),vz(*) ! K=0.0d0 ! do i=1,nom K=K+(vx(i)**2.0+vy(i)**2.0+vz(i)**2.0) enddo ! K=0.5d0*K end subroutine kinetic_energy ! ! Calculate Energy, Virial and Forces ! subroutine force(nstep,n_MD,nom,box,boxinv,rcutsq,U,W,x,y,z,fx,fy,fz, myrank,nprocs) implicit none include ’mpif.h’ ! preprocessor directive ! ! Declaration of Variables
101
&
! integer :: nstep,n_MD,nom,i,j real(8) :: box,boxinv,rcutsq,U,W,xij,yij,zij,fxij,fyij,fzij, rijsq,sr2,sr6,uij,wij,fij real(8) :: x(*),y(*),z(*), & fx(*),fy(*),fz(*) integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,irank real(8) :: qU,qW real(8) :: qfx(nom),qfy(nom),qfz(nom) ! if(nstep == n_MD+1) then ! Energy, Virial and Forces ! ! Zeroth Energy, Virial and Forces ! qU=0.0d0 qW=0.0d0 do i=1,nom qfx(i)=0.0d0 qfy(i)=0.0d0 qfz(i)=0.0d0 enddo ! ! Outer Cycle ! irank=-1 ! do i=1,nom-1 ! ! Inner Cycle ! do j=i+1,nom irank=irank+1 if(irank == nprocs) irank=0 if(myrank == irank) then ! ! Distance between Molecules ! xij=x(i)-x(j) yij=y(i)-y(j) zij=z(i)-z(j) ! ! Periodic Boundary Conditions ! xij=xij-dnint(xij*boxinv)*box yij=yij-dnint(yij*boxinv)*box zij=zij-dnint(zij*boxinv)*box ! ! Cut-Off ! rijsq=xij*xij+yij*yij+zij*zij if(rijsq < rcutsq) then ! ! Auxiliary Variables ! sr2=1.0d0/rijsq sr6=sr2*sr2*sr2 ! ! Potential Energy ! uij=4.0d0*sr6*(sr6-1.0d0) qU=qU+uij ! ! Virial
102
&
! wij=(-1.0d0)*24.0d0*sr6*(2.0d0*sr6-1.0d0) qW=qW-(wij/3.0d0) ! ! Forces ! fij=(-1.0d0)*(wij/rijsq) ! - w(r_ij)/r_ij^2 fxij=fij*xij fyij=fij*yij fzij=fij*zij qfx(i)=qfx(i)+fxij qfy(i)=qfy(i)+fyij qfz(i)=qfz(i)+fzij qfx(j)=qfx(j)-fxij qfy(j)=qfy(j)-fyij qfz(j)=qfz(j)-fzij endif endif enddo enddo ! ! Add Results Together ! call MPI_ALLREDUCE(qfx,fx,nom,MPI_REAL8,MPI_SUM, MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfy,fy,nom,MPI_REAL8,MPI_SUM, MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfz,fz,nom,MPI_REAL8,MPI_SUM, MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qU,U,1,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qW,W,1,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) else ! Forces ! ! Zeroth Forces ! do i=1,nom qfx(i)=0.0d0 qfy(i)=0.0d0 qfz(i)=0.0d0 enddo ! ! Outer Cycle ! irank=-1 ! do i=1,nom-1 ! ! Inner Cycle ! do j=i+1,nom irank=irank+1 if(irank == nprocs) irank=0 if(myrank == irank) then ! ! Distance between Molecules ! xij=x(i)-x(j) yij=y(i)-y(j) zij=z(i)-z(j) ! ! Periodic Boundary Conditions ! xij=xij-dnint(xij*boxinv)*box
& & &
103
yij=yij-dnint(yij*boxinv)*box zij=zij-dnint(zij*boxinv)*box ! ! Cut-Off ! rijsq=xij*xij+yij*yij+zij*zij if(rijsq < rcutsq) then ! ! Auxiliary Variables ! sr2=1.0d0/rijsq sr6=sr2*sr2*sr2 ! ! Virial ! wij=(-1.0d0)*24.0d0*sr6*(2.0d0*sr6-1.0d0) ! ! Forces ! fij=(-1.0d0)*(wij/rijsq) ! - w(r_ij)/r_ij^2 fxij=fij*xij fyij=fij*yij fzij=fij*zij qfx(i)=qfx(i)+fxij qfy(i)=qfy(i)+fyij qfz(i)=qfz(i)+fzij qfx(j)=qfx(j)-fxij qfy(j)=qfy(j)-fyij qfz(j)=qfz(j)-fzij endif endif enddo enddo ! ! Add Results Together ! call MPI_ALLREDUCE(qfx,fx,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfy,fy,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) call MPI_ALLREDUCE(qfz,fz,nom,MPI_REAL8,MPI_SUM, & MPI_COMM_WORLD,ierr) endif end subroutine force ! ! Leap-Frog Algorithm ! subroutine Leap_Frog(nom,nstep,n_MD,dt,box,boxinv,x,y,z,vx,vy,vz,fx,fy,fz) implicit none ! ! Declaration of Variables ! integer :: nom,nstep,n_MD,i real(8) :: dt,box,boxinv real(8) :: x(*),y(*),z(*), & vx(*),vy(*),vz(*), & fx(*),fy(*),fz(*) ! ! Integration ! if(nstep == 1) then ! Case nstep=1 do i=1,nom ! ! v(to+dt/2)=v(to)+(dt/2)*[f(to)/m] !
104
vx(i)=vx(i)+(dt/2.0d0)*fx(i) vy(i)=vy(i)+(dt/2.0d0)*fy(i) vz(i)=vz(i)+(dt/2.0d0)*fz(i) ! ! r(to+dt)=r(to)+dt*v(to+dt/2) ! x(i)=x(i)+vx(i)*dt y(i)=y(i)+vy(i)*dt z(i)=z(i)+vz(i)*dt ! ! Periodic Boundary Conditions ! x(i)=x(i)-dnint(x(i)*boxinv)*box y(i)=y(i)-dnint(y(i)*boxinv)*box z(i)=z(i)-dnint(z(i)*boxinv)*box enddo else if(nstep == n_MD+1) then ! Case nstep=n_MD+1 do i=1,nom ! ! v(tn)=v(tn-dt/2)+(dt/2)*[f(tn)/m] ! vx(i)=vx(i)+(dt/2.0d0)*fx(i) vy(i)=vy(i)+(dt/2.0d0)*fy(i) vz(i)=vz(i)+(dt/2.0d0)*fz(i) enddo else ! Otherwise do i=1,nom ! ! v(t+dt/2)=v(t-dt/2)+dt*[f(t)/m] ! vx(i)=vx(i)+dt*fx(i) vy(i)=vy(i)+dt*fy(i) vz(i)=vz(i)+dt*fz(i) ! ! r(t+dt)=r(t)+dt*v(t+dt/2) ! x(i)=x(i)+vx(i)*dt y(i)=y(i)+vy(i)*dt z(i)=z(i)+vz(i)*dt ! ! Periodic Boundary Conditions ! x(i)=x(i)-dnint(x(i)*boxinv)*box y(i)=y(i)-dnint(y(i)*boxinv)*box z(i)=z(i)-dnint(z(i)*boxinv)*box enddo endif end subroutine Leap_Frog ! ! Zeroth Thermo Accumulators ! subroutine zero(sum_u,sum_P,ssq_u,ssq_P) implicit none ! ! Declaration of Variables ! real(8) :: sum_u,sum_P,ssq_u,ssq_P ! ! Zeroth sum Accumulators ! sum_u=0.0d0 sum_P=0.0d0 ! ! Zeroth ssq Accumulators !
105
ssq_u=0.0d0 ssq_P=0.0d0 end subroutine zero ! ! Termodynamic Quantities ! subroutine thermo(nom,temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) implicit none ! ! Declaration of Variables ! integer :: nom real(8) :: temp,dens,vol,U,U_LRC,W,W_LRC,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P ! ! Instantaneous Thermodynamic Quantities ! u_ac=(U+U_LRC)/dble(nom) ! Energy per Particle P_ac=dens*temp+(W+W_LRC)/vol ! Pressure ! ! sum Accumulators ! sum_u=sum_u+u_ac sum_P=sum_P+P_ac ! ! ssq Accumulators ! ssq_u=ssq_u+u_ac**2.0 ssq_P=ssq_P+P_ac**2.0 end subroutine thermo ! ! Output ! subroutine output(kk,ncycl,nrun,nacpt,nom,rcut,temp,dens,vol,box,u_ac,P_ac, & sum_u,sum_P,ssq_u,ssq_P) implicit none ! ! Declaration of Variables ! integer :: kk,ncycl,nrun,nacpt,nom real(8) :: rcut,temp,dens,vol,box,u_ac,P_ac,rat_acp,sum_u,sum_P,ssq_u,ssq_P, & avr_u,avr_P,flc_u,flc_P ! ! Initial Output ! if(kk == 1) then write(*,*) " *** Initial Output ***" write(*,fmt="(a,e13.5)") "Temperature: ",temp write(*,fmt="(a,e13.5)") "Density: ",dens write(*,*) write(*,fmt="(a,i4)") "# of Molecules: ",nom write(*,fmt="(a,e13.5)") "Cut-Off Distance: ",rcut write(*,fmt="(a,e13.5)") "Volume of Box: ",vol write(*,fmt="(a,e13.5)") "Length of Box: ",box return endif ! ! Periodic Output ! if(kk == 2) then write(*,*) write(*,fmt="(a,i8)") " # of Cycles:",ncycl rat_acp=dble(nacpt)/dble(ncycl) write(*,fmt="(a,e13.5)") " Acceptance Ratio:",rat_acp
106
write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") u_ac,P_ac write(*,fmt="(8(e13.5))") sum_u/dble(ncycl),sum_P/dble(ncycl) return endif ! ! Final Output ! if(kk == 3) then ! ! Ensemble Averages ! avr_u=sum_u/dble(nrun) avr_P=sum_P/dble(nrun) ! ! Ensemble Fluctuations ! flc_u=dsqrt(dabs((ssq_u/dble(nrun))-avr_u**2.0)) flc_P=dsqrt(dabs((ssq_P/dble(nrun))-avr_P**2.0)) write(*,*) write(*,*) " *** Final Output ***" write(*,*) " Ensemble Averages" write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") avr_u,avr_P write(*,*) " Ensemble Fluctuations" write(*,fmt="(5(a))") " u* "," P*" write(*,fmt="(8(e13.5))") flc_u,flc_P return endif end subroutine output ! ! Read from File ! subroutine readcn(nmax,n_MD,x,y,z,vx,vy,vz,cnfile) implicit none ! ! Declaration of Variables ! integer :: nmax,n_MD real(8) :: x(nmax),y(nmax),z(nmax), & vx(nmax),vy(nmax),vz(nmax) integer,parameter :: cnunit=10 character(len=*) :: cnfile ! ! Unformatted File ! open(unit=cnunit,file=cnfile,status="old",form="unformatted") read(cnunit) n_MD read(cnunit) x,y,z read(cnunit) vx,vy,vz close(unit=cnunit) end subroutine readcn ! ! Write to File ! subroutine writcn(nmax,n_MD,x,y,z,vx,vy,vz,cnfile) implicit none ! ! Declaration of Variables ! integer :: nmax,n_MD real(8) :: x(nmax),y(nmax),z(nmax), & vx(nmax),vy(nmax),vz(nmax) integer,parameter :: cnunit=10 character(len=*) :: cnfile
107
! ! Unformatted File ! open(unit=cnunit,file=cnfile,status="unknown",form="unformatted") write(cnunit) n_MD write(cnunit) x,y,z write(cnunit) vx,vy,vz close(unit=cnunit) end subroutine writcn ! ! RANDOM VARIATE FROM THE STANDARD NORMAL DISTRIBUTION. ! THE DISTRIBUTION IS GAUSSIAN WITH ZERO MEAN AND UNIT VARIANCE. ! REFERENCE: KNUTH D, THE ART OF COMPUTER PROGRAMMING, ! (2ND ED. ADDISON-WESLEY), 1978 ! function gauss(idum) implicit none ! ! Declaration of Variables ! integer :: idum,I real(8) :: gauss,random real(8) :: SUM,R,R2 real(8),parameter :: A1=3.949846138d0,A3=0.252408784d0, & A5=0.076542912d0,A7=0.008355968d0, & A9=0.029899776d0 ! SUM=0.0d0 do I=1,12 SUM=SUM+random(idum) enddo ! R=(SUM-6.0d0)/4.0d0 R2=R*R gauss=((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R end function gauss ! ! RANDOM GENERATOR (Numerical Recipes in Fortran90) ! FUNCTION random(idum) ! ! "Minimal" random number generator of Park and Miller combined with a Marsaglia shift ! sequence. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint ! values). This fully portable, scalar generator has the "traditional" (not Fortran 90) calling ! sequence with a random deviate as the returned function value: call with idum a negative ! integer to initialize; thereafter, do not alter idum except to reinitialize. The period of this ! generator is about 3.1 x 10^18. ! IMPLICIT NONE ! ! Declaration of Variables ! INTEGER,PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B),INTENT(INOUT) :: idum REAL :: ran INTEGER(K4B),PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836 REAL,SAVE :: am INTEGER(K4B),SAVE :: ix=-1,iy=-1,k ! real(8) :: random ! ! Initialize ! if((idum <= 0).or.(iy < 0)) then am=nearest(1.0,-1.0)/IM
108
iy=ior(ieor(888889999,abs(idum)),1) ix=ieor(777755555,abs(idum)) ! ! Set idum Positive ! idum=abs(idum)+1 endif ! ! Marsaglia Shift Sequence with Period 2^32 - 1 ! ix=ieor(ix,ishft(ix,13)) ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) ! ! Park-Miller Sequence by Schrage’s Method, Period 2^31 - 2 ! k=iy/IQ iy=IA*(iy-k*IQ)-IR*k if(iy < 0) iy=iy+IM ! ! Combine the Two Generators with Masking to Ensure Nonzero Value ! ran=am*ior(iand(IM,ieor(ix,iy)),1) ! random=dble(ran) END FUNCTION random
109
Reference [1] Y. Aoyama, J. Nakano, RS/6000: Practical MPI Programming, International Technical Support Organization, 1999. [2] W. Gropp, E. Lusk, A. Skjellum, Using MPI. Portable Parallel Programming with the Message-Passing Interface, 2. vydání, MIT Press, 1999. [3] P. S. Pacheco, W. C. Ming, MPI Users’ Guide in FORTRAN, 1997. [4] P. S. Pacheco, A Users’ Guide to MPI, 1998. [5] MPI: A Message-Passing Interface Standard, verze 1, 1995. [6] MPI-2: Extension to the Message-Passing Interface, verze 2, 1997. [7] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 1. vydání, Clarendon Press, Oxford, 1987. [8] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2. vydání, Academic Press, London, 2002. [9] U. H. H. Hansmann, Y. Okamoto, F. Eisenmenger, „Molecular dynamics, Langevin and hybrid Monte Carlo simulations in a multicanonical ensemble”, Chem. Phys. Letts. 259, 321, 1996.
110
