parallel tempering (병렬 조질) [algorithm] by 바죠

Parallel Tempering
http://en.wikipedia.org/wiki/Parallel_tempering

Swendsen RH and Wang JS (1986) Replica Monte Carlosimulation of spin glasses Physical Review Letters 57 : 2607-2609


simulated annealing과 더불어서 매우 유용한 계산 방법으로 알려져 있다. 간단하게 말하면 simulated annealing을여러 가지 온도에서 동시에 수행하는 알고리듬이다. 뿐만 아니라, 서로 다른 온도들에서 만들어진 분자구조들은 확률적으로 바꿔 가면서 계속해서 시뮬레이션 하는 특징이 있다.  http://iopscience.iop.org/article/10.1088/0953-8984/27/23/233102


병렬 계산에 아주 적합한 계산방법이다. 여러 개의 독립된 계산을 병행해서 수행함으로써 굉장한 시너지를 유발해 낼 수 있는 특징이있다. 각각의 온도에서 보다 더 정확한 샘플링이 구축된다. 여러 온도 계산을 동시에 수행한다. 알고리듬에서 핵심적으로 활용되는 공식은 아래와 같다.

 p = \min \left( 1, \frac{ \exp \left( -\frac{E_j}{kT_i} - \frac{E_i}{kT_j} \right) }{ \exp \left( -\frac{E_i}{kT_i} - \frac{E_j}{kT_j} \right) } \right) = \min \left( 1, e^{(E_i - E_j) \left( \frac{1}{kT_i} - \frac{1}{kT_j} \right)} \right) ,


기본적으로 다중 시뮬레이션이다. 여러 온도에서의 시뮬레이션을 동시에 진행하는 것이 기본 시작점이다. 여러 온도에서의 앙상블을 각자 계산하는 것이 기본 세팅이다. 하지만, 동시에 수행하면서 상승효과 시너지 효과를 볼 수 있다는 점에 유의 해야만 한다.  시스템 크기가 커질 경우, 보다 많은 온도들에서 샘플링을 해야만한다.  그렇지 않을 경우, 서로 다른 온도들 사이의 맞교환 확률이 낮아 진다. 맞교환 확률은 해당 상태밀도 함수들의 상당한 중복이 있을 때 허용된다. 즉, 서로 다른 온도들에서 얻어낼 수 있는 상태밀도 함수들에서 중첩이 있어야만한다.

주어진 모든 온도들 각자의 앙상블을 더욱 더 정확하게 계산할 수 있다. 시너지 효과가 있다는 점에 주의해야 한다.

또한, 이 시뮬레이션은 global optimization의방법으로 해석할 수도 있다.
http://incredible.egloos.com/4784592

매트로폴리스 알고리듬에서로 다른 온도 분포가 있다는 점에서 simulated annealing의 변형으로 볼 수 있다. 사실은 super simulated annealing이라고 부를 수 있다. 왜냐하면, 온도 분포를 잘 세팅 하고 고정적/지속적으로높은 온도에서 만들어진 conformation을 낮은 온도에 공급할 수 있다. 따라서, 온도를 낮추는 방식이 고정되어 있는 simulated annealing이라고 볼 수 있다. global optimization을 수행하는 하나의 방법이다. 오래된 알고리듬이다. 하나의 구조(입자)가 특정 구조(위치)에서 에너지를 주고, 또다른 구조에서 또 다른 에너지를 줄 경우, 온도에 따라서 이 새로운 구조를 수용하거나 거부한다. 메트로폴리스 알고리듬을 그대로 수용한다. 고정된 온도 분포들을 가정한다. 온도가 천천히 바뀌는 것이 아니다. 고정된 온도들이 있다. 이 부분에서 simulated annealing 방법과 차이가 난다. 또한, 여러 가지 conformation을 동시에다루는 것을 의미한다. 최소한 각 온도에서 하나의conformation이 변화함을 의미한다.

주어진 온도에서 충실히 매트로폴리스 알고리듬을 수행한다. 서로 다른 온도에 적응된 conformation들을 바꾸어주는 작업이 추가된다. T1의 conformation이T2에서의 conformation과 서로 바꾸어지는 동작이 가능할 수 있다. 매우 유사한 에너지를 가진다면 그것은 가능하다. 또한 T1 T2차이가 그렇게 크지 않을 경우에도 가능하다. 실제로는 매트로폴리스 알고리듬으로 돌아 갈 필요가 있다. 메트로폴리스 알고리듬에서, 특히, 온도 T1 에서 conformation이 변환될 때, 즉, X ---> X' 으로 변환 가능한지를 따지게된다. 이 때, conformation X' 이 다른 온도(T2 /= T1)에서 만들어진 것으로 대체할 수 있다는 것이다.
http://en.wikipedia.org/wiki/Metropolis_algorithm





       module pt_mod
!      Written by In-Ho Lee, KRISS, December, 2012.
       implicit none
       private
       save
       integer ndata,ndeg, nseed,nsim,ntem, nsize
       integer, allocatable :: iseed(:)
       real*8 obj_lowest
       real*8 ptact,ptrjt
       real*8, allocatable :: xpt(:),ypt(:)
       real*8, allocatable ::qq(:,:),pp(:),qq_lowest(:),xtem(:)
       real*8, allocatable :: obj(:)
       real*8 poly


       public ::initial,final,do_pt


       contains


       subroutineinitial(nsim0,ntem0,nseed0,temp0)
       implicit none
       integer nsim0,ntem0,nseed0
       real*8 temp0
       real*8 tmp,t1,t2,dt
       integer k,item
       logical lexist


       nsim=nsim0
       ntem=ntem0
       nseed=nseed0

       call random_seed(size=nsize)
       allocate(iseed(nsize))
       do k=1,nsize
       iseed(k)=nseed+k
       enddo
       call random_number(put=iseed)
       deallocate(iseed)


       allocate(xtem(ntem))
       allocate(obj(ntem))


       t1=temp0
       t2=t1*(0.5d0)**8
       if(t2 < 1.d-6) t2=1.d-6
       dt=(t2-t1)/float(ntem-1)
       do k=1,ntem
       xtem(k)=t1+dt*float(k-1)
       print*, xtem(k)
       enddo


       ndata=100
       allocate(xpt(ndata))
       allocate(ypt(ndata))
       call preparation()


       ndeg=2
       allocate(pp(0:ndeg))
       allocate(qq_lowest(0:ndeg))
       allocate(qq(0:ndeg,ntem))
       qq=0.d0
       do item=1,ntem
       do k=0,ndeg
       call random_number(tmp)
       tmp=tmp-0.5d0
       qq(k,item)=tmp
       enddo
       enddo



       inquire(file='fort.11',exist=lexist)
       if(lexist)then
       print*, 'file fort.11 is present'
       open(11,file='fort.11',form='formatted')
       read(11,*) item
       read(11,*) qq_lowest
       read(11,*) tmp
       close(11)
       do item=1,ntem
       qq(:,item)=qq_lowest(:)
       enddo
                endif


       end subroutine initial


       subroutine final()
       implicit none
       integer j


      open(7,file='fort.7',form='formatted')
       do j=1,ndata
       write(7,*) xpt(j),ypt(j)
       enddo
       write(7,*) '&'
       do j=1,ndata
       write(7,*)xpt(j),poly(xpt(j),qq_lowest,ndeg)
       enddo
       close(7)


      open(11,file='fort.11',form='formatted')
       write(11,*) ndeg
       write(11,*) qq_lowest
       write(11,*) obj_lowest
       close(11)


       deallocate(xpt,ypt)
       deallocate(xtem,obj)
       deallocate(qq,pp,qq_lowest)
       end subroutine final


       subroutinepreparation()
       implicit none
       real*8 tmp,xi,xf,dx
       integer j


       xi=0.d0
       xf=10.d0
       dx=(xf-xi)/float(ndata-1)
       do j=1,ndata
       xpt(j)=xi+dx*float(j-1)
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*2.0d0
       ypt(j)= 0.1d0 +0.2d0*xpt(j)+0.3d0*xpt(j)**2    +tmp
       enddo
       end subroutine preparation


       subroutineperturb1(item)
       implicit none
       integer item
       integer k
       real*8 tmp


       pp(:)=qq(:,item)
       call random_number(tmp)
       k=tmp*(ndeg+1)
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*sqrt(xtem(item))
       pp(k)=pp(k)+tmp
       end subroutine perturb1


       subroutineperturb2(item)
       implicit none
       integer item
       integer k
       real*8 tmp


       pp(:)=qq(:,item)
       do k=0,ndeg
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*sqrt(xtem(item))
       pp(k)=pp(k)+tmp
       enddo


       end subroutineperturb2


       subroutine do_pt()
       implicit none
       integer, parameter :: ntry=2
       real*8 act0(0:ntry-1),rjt0(0:ntry-1)
       real*8 act,rjt
       real*8 objqq,objpp,dell, tmp
       integer item,isim,itry,k,ii
       integer nnpt


       callcal_obj(qq(1,1),objqq)
       obj_lowest=objqq ; qq_lowest(:)=qq(:,1)


       ptact=0.d0 ;ptrjt=0.d0
       nnpt=100
       do ii=1,nnpt


       do item=1,ntem
!      print*, xtem(item),'xtem'
       act=0.d0 ; rjt=0.d0  ; act0=0.d0 ;rjt0=0.d0
       call cal_obj(qq(1,item),objqq)


       do isim=1,nsim/100
       itry=mod(isim,ntry)


       if(itry ==0)   then
       call perturb1(item)
       elseif(itry ==1)then
       call perturb2(item)
                      endif


       call cal_obj(pp,objpp)
       dell=(objpp-objqq)/xtem(item)
       if(dell > 50.d0) dell=50.d0
       if(dell < -50.d0) dell=-50.d0
       call random_number(tmp)
       if(tmp < min(1.d0,exp(-dell)))then
       qq(:,item)=pp(:) ; objqq=objpp
       act=act+1.d0
       act0(itry)=act0(itry)+1.d0
                                    else
       rjt=rjt+1.d0
       rjt0(itry)=rjt0(itry)+1.d0
                                    endif
       obj(item)=objqq


!      if(item == ntem)then
       if(objqq < obj_lowest)then
       obj_lowest=objqq ; qq_lowest(:)=qq(:,item)
                            endif
!                     endif


       enddo
!      do itry=0,ntry-1
!      write(6,'(f18.8,2x,a3,2x,2f20.3)')act0(itry)/(act0(itry)+rjt0(itry)),'(%)',act0(itry),rjt0(itry)
!      enddo
!      write(6,'(f18.8,2x,a3,2x,2f20.3)') act/(act+rjt),'(%)',act,rjt
!      print*, xtem(item),objqq,obj_lowest
       enddo

       call do_exchange()


       enddo
       print*, ptact/(ptact+ptrjt),ptact,ptrjt
       print*, obj_lowest
       print*, qq_lowest


       end subroutine do_pt


       subroutinedo_exchange()
       implicit none
       integer i,j
       real*8 dell,tmp


       do i=1,ntem-1
       j=i+1
      dell=-(1.d0/xtem(i)-1.d0/xtem(j))*(obj(i)-obj(j))
       if(dell > 50.d0) dell=50.d0
       if(dell < -50.d0) dell=-50.d0
       call random_number(tmp)
       if(tmp < min(1.d0,exp(-dell)))then
       pp(:)=qq(:,i)
       qq(:,i)=qq(:,j)
       qq(:,j)=pp(:)
       ptact=ptact+1.d0
                                    else
       ptrjt=ptrjt+1.d0
                                    endif
       enddo


       end subroutinedo_exchange


       subroutine cal_obj(q,objq)
       implicit none
       real*8 q(0:ndeg),objq
       integer i


       objq=0.d0
       do i=1,ndata
       objq=objq+(ypt(i)-poly(xpt(i),q,ndeg))**2
       enddo
       objq=objq/float(ndata)
       end subroutine cal_obj


       real*8 function poly(x,q,ndeg)
       implicit none
       integer ndeg
       real*8 x,q(0:ndeg)
       real*8 tmp
       integer i


       tmp=0.d0
       do i=0,ndeg
       if(i == 0)then
       tmp=tmp+q(i)
                else
       tmp=tmp+q(i)*x**float(i)
                endif
       enddo
       poly=tmp
       end function poly


       end module pt_mod
!234567890
       program test
       USE pt_mod, ONLY : initial,final,do_pt
       implicit none
       integer nsim0,ntem0,nseed0
       real*8 temp0
       integer itemp,irate,itemq


       callsystem_clock(itemp,irate)


       nsim0=100000
       ntem0=10
       nseed0=4321
       temp0=1.d0


       callinitial(nsim0,ntem0,nseed0,temp0)


       call do_pt()


       call final()


       callsystem_clock(itemq)
       write(6,'(2e15.8,2x,a6)')float(itemq-itemp)/float(irate)/60., float(itemq-itemp)/float(irate)/3600.,'min,h'
       end program test


 



복제품 맞바꿈 시뮬레이션:
http://incredible.egloos.com/4587643


그림 출처:
http://pubs.rsc.org/en/content/articlepdf/2005/cp/b515592b/2005-11-16?page=search

 




 







핑백


최근 포토로그