Metoda paralelního „temperingu” - Často používané MPI příkazy

2.8 Často používané MPI příkazy

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. Jednot-livé MC procesy C

₁

, C

₂

, . . ., C

seřadíme tak, že

T

₁

> T

₂

> ... > T

(β

₁

< β

₂

< ... < β

)

Základní myšlenkou paralelního „temperingu” je provádět vedle standardních MC „po-hybů” („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 {C

}

ⁿ_i=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 C

z {C

}

ⁿ⁻¹_i=1

;

• když

exp (∆) ≥ ζ (16)

provede se výměna konfigurací C

↔ C

_i+1

. V rovnici (16) je

∆ = ∆β∆U (17)

∆β = β

_i+1

− β

(18)

∆U = U

_i+1

− U

(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 C

z {C

}

ⁿ⁻¹_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 (C

_i+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í C

↔ C

_i+1

t.j. x(C

) ↔ x(C

_i+1

) a U (C

) ↔ U (C

_i+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 vypa-dat 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

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

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)

if(nprocs /= notemp) stop "ParPT: # of processes /= # of temperatures!"

! Auxiliary Variables

dxmax=0.1d0 dx=0.025d0

! Assign Actual Temperature

! Iput Values (Hard Wired)

x=-1.2d0 ! Initial x (1st valley) beta=1.0d0/temp ! 1/(kB.T)

! Zero Accumulators

! 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

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))

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

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

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

! 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 · 10

⁸

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 ∆x

_max

= 0.1 a

∆x = 0.025 a MC procesy začaly s částicí v x

ⁱⁿⁱ

= −1.2. MC procesy probíhaly při čtyřech teplotách: k

T = 2, k

T = 0.3, k

T = 0.1 a k

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.

Obrázek 16: Pravděpodobnost nalezení částice v místě x, P (x), pro teploty: k

T = 2,

k

T = 0.3, k

T = 0.1 a k

T = 0.05; k

je Boltzmannova konstanta. P (x) byla získána

z 1 · 10

⁸

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 ∆x

_max

= 0.1 a ∆x = 0.025 a Monte Carlo

procesy začaly s částicí v x

ⁱⁿⁱ

= −1.2.

In document Univerzita Jana Evangelisty Purkyně v Ústí nad Labem Přírodovědecká fakulta (Stránka 49-57)