복제품-맞바꿈 분자동역학 (replica-exchange molecular dynamics; REMD) [algorithm] by 바죠

분자동역학 (molecular dynamics) 입문


복제품-맞바꿈 분자동역학(replica-exchange molecular dynamics; REMD) parallel tempering이라고도 불리우는 계산 방법이다.  주로 주어진 온도에서 관심있는 시스템의 열역학적 물리량을 계산할 때 사용된다.

(1)
몇 개의 서로 다른 온도{ T_1, T_2, T_3, ..., T_M }분포들에서 독립적인 분자동역학(또는 Monte Carlo) 계산을 동시에 수행한다.
{ T_1 <T_2 < T_3, ..., < T_M }분포가 여러 개 있어도 좋다. 1개만 있는 경우가  통상의 REMD/REMC 방식이다. 동일한 온도 분포 세트가 여러 개 있는 경우를 특별히 Multiplexed REMD/REMC라고 한다.
(2)
특정 시간(또는 Monte Carlo의 경우 스텝 수) 이후에 인접한 온도들 사이의 분자 모양들을, Metropolis 형식으로, 확률적으로 직접 맞바꾸어 준다. 이러한 계산들을 반복한다. 이웃한 두 개의 온도 분포에 대해서 분자 모양의 맞바꿈 시도를 준비한다.  이웃한 두 개의 온도는 랜덤으로 선정한다. 현실적으로는 인접한 온도 사이에서만 높은 복제품 맞교환 확률이 높다.
 

다양한 온도 분포에서 각각 sampling을 수행한다. 어차피 각 온도에서의 시스템의 특성을 조사하는 것이 목표이다.  이러한 계산에서의 sampling과 관련된 계산 수렴성은 일반적인 단일 온도, 분자동역학 (또는 Monte Carlo) 계산에서 얻어지는 sampling의 계산 수렴성 보다 뛰어나다.  그 이유는 분자들이 특히, 상대적으로 낮은 온도 분포에서 에너지 골짜기에서 쉽게 빠져 나오지 못하는 볼쯔만 분포의 특성 때문이다. 결국, 높은 온도에서 상대적으로 자유롭게 포텐셜 에너지 골짜기를 빠져 나온 분자 구조가 낮은 온도 분포에, REMD/REMC에서는, 자연스럽게 스며들게 하기 때문에 전체 계산의 수렴성은 높아지게 된다.  보다 넓은 conformation space를 효율적으로 탐색하게 해준다.

높은 온도 분자동역학 계산에서 자동적으로 새롭고 다양한 분자 구조를 지속적으로, 확률적으로, 낮은 온도 분자동역학 계산에 제공하기 때문에 사실상 simulated annealing (SA), 모의 담금질의 특성을 가지고 있다. 통상 SA에서는 온도를 줄여나가면서 분자 구조를 결정하지만, 이러한 온도 조절 자체가 까다로운 계산 파라미터가 된다. 담금질 온도 설정에 따라서 계산이 달라질 수 있다.  또한, SA에서는 단 하나의 복제품만을 사용한다.  결국, REMD 계산은 초 모의 담금질, Super SA로 불러도 무방하게 된다. 즉, global optimization 방법으로 활용될 수 있는 특징이 있다.  

특별히, 여러 개의 온도 분포들이 존재할 수 잇다.  이렇게 할 경우 병렬 효율성은 증대될 것이다.  이러한 상황에서도 마찬가지로 REMD/REMC 계산 방식은 그대로 적용될 수 있다. 결국, 독립적인 REMD/REMC를 수행하는 것에 지나지 않기 때문이다. 병렬 효율성이 극대화 되는 장점도 있다. 결국, 서로 다른 온도 분포들 사이에서 이웃한 온도들 사이의 맞바꿈은 정당하게, 동등하게 여전히 적용될 수 있는 것이다. 이것이  Multiplexed REMD라는 방법으로 알려진 것이다.  {T_1, T_2, T_3,..., T_M}
 Refernce: Understanding Molecular Simulation, D. Frenkel and B. Smit, p. 391   
  A typo is found in the book.  min(1.d0, exp(+(b_i-b_j)*(e_i-e_j))) 분자 구조들을 서로 바꾸는것과 온도를 서로 바꾸는것 사이에 차이가 없다. 결국, 위의 식은 이러한 상황을 만족하게 되어 있다.


위에서 보인 그림들은 1차원 포텐셜 하에서 움직이는 하나의 입자를 생각한 경우이다. 서로 다른 온도 분포에서 생각하고 있는 하나의 입자가 발견될 분포를 계산하고자 한다.

계산된 확률은 결국 exp(-E/(kBT)) 에 비례하는 것이 되어야 할 것이다. 하지만, 실제, Monte Carlo, Molecular dynamics 계산에서 이러한 분포를 얻어내기는 매우 어렵다.  그것은 위에 표시된 포텐셜 에너지가 온도에 비해서 상당히 높은 포텐셜 장벽에너지를 가지고 있기 때문이다.

하나의 포텐셜 에너지 골짜기에서 입자가 초기에 위치했다면 어지간히 높은 온도가 아닌 이상 옆쪽의 포텐셜 에너지 골짜기에서 발견되기는 쉽지 않다. 왜냐하면, 입자는 에너지 장벽을 반드시 거쳐서 직접 넘어가야만 하기 때문이다. 양자역학적 터널링 효과는 없다. 이것은 단순히 온도가 매우 높아야만 쉽게 허용되는 것이다. 이렇게 높은 온도를 REMD/REMC에서는 사용할 필요가 있다. 아주 다양한 분자 구조를 얻어 내기 위해서 그렇다.

이러한 현실적인 MC/MD의 문제점/비효율성을 상당 부분 해결해 주는 계산 방법이 바로 REMC/REMD 계산 방법이다.
위의 그림은 입자가 발견될 확률을 REMC 방법으로 온도에 따라서 계산한 결과이다. 볼쯔만 분포에서 예측한 것처럼 동일한 포텐셜 에너지 골짜기에 대해서 거의 동등한 입자 발견 확률 분포를 보여주고 있다. 굉장히 높은 온도에서 움직이는 입자는 상대적으로 쉽게 포텐셜 에너지 장벽을 넘어 갈 수 있다. 그런데, 이런 입자는 바로 인접한 낮은 온도에서 시뮬레이션 되는 분자동역학 계산의 입력으로 들어 갈 수 있다. 마찬가지로, 이러한 상황은 계속해서 낮은 온도 시뮬레이션에서도 일어나게 된다. 결국, 높은 온도에서 만들어진 분자 구조가 낮은 온도에서 시뮬레이션 되어질 수 있다. 통상의 낮은 온도 시뮬레이션에서 만들어지기 쉽지 않은 분자 구조를 입력으로 시뮬레이션이 진행되어 질 수 있다는 결론이 나온다. 최종적으로는 관심이 있는 모든 온도에서의 시뮬레이션이 자연스럽게 진행되는 것을 도와준다. 모든 온도에서 각각 보다 좋은 샘플링을 가능하게해 준다. 열역학적 샘플링에 효율성이 증대된다.  

랜덤으로 선정된 이웃한 두개의 온도 T_i  그리고 T_j에 대해서 분자구조 맞바꿈 시도를 준비한다. 확률적으로 선택되어 지면 실제 맞바꿈 동작을 수행하게 된다.
 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) ,


random number  < p일 때, 분자 구조들의 맞바꿈 작업이 일어 난다.  여기서 random number는 0에서 1.0 사이의 막수이다.

분자구조에 의해서 결정되어지는 E_i, E_j를 서로 바꾸는 것이 T_i, T_j를 서로 바꾸는 것과 다른 차이를 주지 않는다. 이들의 변화에 대칭으로 수식이 적혀져 있는것에 주목할 필요가 있다.  매우 단순한 알고리듬이다. 또한 상대적으로 높은 병렬 효율성을 쉽게 확보할 수 있다. 다양한 응용 연구가 가능하도록 알고리듬이 충분히 단순하다.   

한 마디로 REMC/REMD를 활용하면 "Accelerationg Monte Carlo Sampling"을 이루어낼 수 있다. 다른 말로 표현하면, super simulated annealing으로 볼 수 있다.  왜냐하면, 온도가 천천히 낮아 지는 단계를 허용하는 simulated annealing을 자동으로 포함하고 있다. 그리하여, super simulated annealing이라고 불리울 수 있다. 즉, 낮은 온도에서 얻어진 구조만을 생각하면, 특히, 온도를 굉장히 낮은 온도로 가정하면, 이것은 global optimization의 방법으로 생각할 수 있다.


참고:
http://en.wikipedia.org/wiki/Parallel_tempering
http://en.wikipedia.org/wiki/Simulated_annealing
http://en.wikipedia.org/wiki/Metropolis-Hastings
http://en.wikipedia.org/wiki/Global_optimization
http://en.wikipedia.org/wiki/Monte_Carlo_method


참고 아티클:
1566_article_physics_and_high_technology.pdf




       module remc
!      Written by In-Ho Lee, KRISS/KIAS, December 28 (2008).
!      Refernce: Understanding Molecular Simulation, D. Frenkel and B. Smit, p. 391
!      A typo is found in the book.  min(1.d0, exp(+(b_i-b_j)*(e_i-e_j)))
       implicit none
       private
       save
       integer ntemper,nbin
       real*8 delx,xmin,xmax
       real*8, allocatable :: xtemper(:),prob(:,:)
       real*8, allocatable :: xyz(:),energy(:)
       logical l_first

       public :: remc_initial,remc_final, remc_accrej

       contains

       subroutine remc_initial(n,t1,t2)
       implicit none
       integer n
       real*8 t1,t2
       integer itemper
       real*8 s1,s2,ds
       real ranmar

       nbin=1000 ; xmin=-2.d0 ; xmax= 2.d0 ; delx=(xmax-xmin)/float(nbin-1)

       ntemper=n
       write(6,*) ntemper,'ntemper'
       allocate(xyz(ntemper)) ; allocate(energy(ntemper)) ; allocate(xtemper(ntemper)) ; allocate(prob(nbin,ntemper))

       xyz=-1.d0+ranmar() ; energy=0.d0 ; prob=0.d0

       s1=min(t1,t2) ; s2=max(t1,t2) ; ds=(s2-s1)/float(ntemper-1)
       do itemper=1,ntemper
       xtemper(itemper)=s1+ds*float(itemper-1)
       write(6,'(i8,2f18.8)') itemper,xtemper(itemper),1.d0/xtemper(itemper)
       enddo
       l_first=.true.
       end subroutine remc_initial

       subroutine remc_final()
       implicit none
       integer ibin,itemper
       real*8 tmp

       do itemper=1,ntemper
       tmp=0.d0
       do ibin=1,nbin
       tmp=tmp+prob(ibin,itemper)
       enddo
       if(abs(tmp) < 1.d-9) tmp=1.d0 ; tmp=1.d0/tmp
       prob(:,itemper)=prob(:,itemper)*tmp
       enddo

       open(11,file='fort.11',form='formatted')
       do itemper=1,ntemper
       write(11,'(a,2x,2f18.8)') '#', xtemper(itemper),1.d0/xtemper(itemper)
       do ibin=1,nbin
       tmp=xmin+float(ibin-1)*(xmax-xmin)/float(nbin-1)
       write(11,'(2f18.8)') tmp,prob(ibin,itemper)
       enddo
       if(itemper /= ntemper) write(11,*) '&'
       enddo
       close(11)
       deallocate(xyz) ; deallocate(energy) ; deallocate(xtemper) ; deallocate(prob)
       end subroutine remc_final
       subroutine remc_accrej
       implicit none
       real*8 delu,enew,xnew,beta,dell,delb
       integer itemper
       real ranmar

       do itemper=1,ntemper

       xnew=xyz(itemper)+(ranmar()-0.5)*0.2d0
       call potential(xyz(itemper),energy(itemper))
       call potential(xnew,enew)
       delu=enew-energy(itemper)
       beta=1.d0/xtemper(itemper)
       dell=beta*delu
       if(dell >  50.d0) dell=50.d0
       if(dell < -50.d0) dell=-50.d0
       if( ranmar() < min(1.d0, exp(-dell)) )then
       xyz(itemper)=xnew
       energy(itemper)=enew
                                             endif
       call sampling(xyz(itemper),itemper)

       enddo

!
       if(.not. l_first)   then
       if( ranmar() < 0.10)then

       itemper=int(ranmar()*dble(ntemper-1))+1
       call potential(xyz(itemper),energy(itemper))
       xnew=xyz(itemper+1)
       call potential(xnew,enew)
       delu=enew-energy(itemper)
       delb=1.d0/xtemper(itemper)-1.d0/xtemper(itemper+1)
       dell=delb*delu
       if(dell >  50.d0) dell=50.d0
       if(dell < -50.d0) dell=-50.d0
       if( ranmar() < min(1.d0, exp(-dell) ))then
       delb=xnew
       delu=xyz(itemper)
       xyz(itemper+1)=delu
       xyz(itemper)=delb
       delb=enew
       delu=energy(itemper)
       energy(itemper+1)=delu
       energy(itemper)=delb
                                             endif
                           endif
                           endif
       l_first=.false.
       end subroutine remc_accrej


       subroutine sampling(q,itemper)
       implicit none
       integer itemper
       real*8 q
       integer ibin,i1,i2
       real*8 r1,r2

       ibin=int(q/delx)+1
       i1=nint(q/delx)-3
       i2=nint(q/delx)+3

       i1=1
       i2=nbin

       if(i1 < 1) i1=1 ; if(i2 > nbin) i2=nbin
       do ibin=i1,i2
       r1=xmin+float(ibin-1)*delx
       r2=r1+delx
       if( q >= r1 .and. q < r2)then
       prob(ibin,itemper)=prob(ibin,itemper)+1.d0
                                exit
                                endif
       enddo

       end subroutine sampling

       end module remc

       subroutine potential(q,e)
       implicit none
       real*8 q,e
       real*8 pi,arg

       pi=4.d0*atan(1.d0)
       arg=2.d0*pi*q

       if(q < -2.d0)                            then
       e=1.d10
       elseif( q >= -2.d0   .and. q <= -1.25d0 )then
       e=1.d0*(1.d0+sin(arg))
       elseif( q >  -1.25d0 .and. q <= -0.25d0 )then
       e=2.d0*(1.d0+sin(arg))
       elseif( q >  -0.25d0 .and. q <=  0.75d0 )then
       e=3.d0*(1.d0+sin(arg))
       elseif( q >   0.75d0 .and. q <=  1.75d0 )then
       e=4.d0*(1.d0+sin(arg))
       elseif( q >   1.75d0 .and. q <=  2.00d0 )then
       e=5.d0*(1.d0+sin(arg))
                                                else
       e=1.d10
                                                endif
       return
       end
       program parallel_tempering_test
       USE remc, ONLY :  remc_initial,remc_final, remc_accrej
       implicit none
       integer n
       real*8 t1,t2,q,e
       integer iseed1,iseed2,isim,jsim,nsim
       integer itemp,irate,itemq

       call timestamp ( ) ; call system_clock(itemp,irate)

       iseed1=933 ; iseed2=339

       n=10
       t1=0.05d0
       t2=2.0d0

       call rmarin(iseed1,iseed2)

       open(19,file='fort.19',form='formatted')
       do isim=1,1000
       q=(-2.d0)+float(isim-1)*(2.d0-(-2.d0))/float(1000-1)
       call potential(q,e)
       write(19,'(2f18.8)') q,e
       enddo
       close(19)

       call remc_initial(n,t1,t2)

       nsim=1000000

       do isim=1,20
       do jsim=1,nsim
       call remc_accrej()
       enddo
       enddo

       call remc_final()

       call timestamp ( ) ; call system_clock(itemq)
       write(6,'(2e15.8,2x,a6)') float(itemq-itemp)/float(irate)/60., float(itemq-itemp)/float(irate)/3600.,' min,h'
       end program parallel_tempering_test


------------------------------------------------------------------------------------------------
유전 알고리듬: http://incredible.egloos.com/4856511


Conformational Space Annealing (CSA) : http://incredible.egloos.com/6208916

------------------------------------------------------------------------------------------------




핑백

덧글

  • dyanos 2011/06/14 14:18 # 삭제

    좋은글감사합니다 ㅎㅎ
  • 바죠 2011/06/14 15:31 #

    도움이 되셨나요?
    좋은글은 아닌것 같은데.....
※ 로그인 사용자만 덧글을 남길 수 있습니다.

최근 포토로그