Action-CSA [algorithm, 컴퓨터 프로그램] by 바죠

Action-CSA

미시세계에서 원자들이 만드는 이동경로를 전산모사 하는 것은 아주 어려운 문제 중 하나이다. 특히, 화학반응같이 원자들이 집단적으로 움직이는 상황에 해당하는 전산모사는 더욱더 어려운 것이다. 왜냐하면, 화학반응은 순간적으로 일어나는 것이기 때문이다. 아주 오랜 시간 동안 원자들은  한 자리에 머물게 된다. 그러다가 어느 순간 갑자기 원자들은 집단적으로 이동해버린다. 실험적 측정에서도 마찬가지이다. 너무나도 빠른 시간 안에 화학반응이 순간적으로 일어나기 때문에 이를 관측하는 것은 매우 어렵다. 또한 자주 일어나는 것이 아니다. 아주 가끔씩 빠른 시간 안에 일들이 일어나버리는 것이다. 전산모사의 경우, 너무나도 오랜 시간을 기다려도 보고자 하는 것을 못 보는 경우가 대부분이다. 실험에서는 너무나도 희귀하게 일어나고 동시에 너무 빨라서 측정을 못하는 것이다.

분자동역학(molecular dynamics) 계산 방법은 주어진 온도, 주어진 압력에서 상호작용하는 원자들의 이동을 관찰할 수 있는 최고의 프로토콜이다. 아주 긴 시간(msec ~ sec)을 전산모사 해야만 단백질 접힘과 같은 화학반응을 전산모사 할 수 있다. 여전히 부족하지만, 그나마 화학반응을 관찰할 수 있는 기회를 긴시간을 통해서 확보하는 것이다. 이마저도 확률적인 것이어서 분자동역학 계산 중에 확실하게 기대하는 화학반응이 일어난다고 말할 수 없다. 한 가지 아이디어는 매우 높은 온도에서 일어나는 현상들을 관측하고 낮은 온도에서도 그러한 일들이 일어날 것이라는 추측은 할 수 있다. 하지만, 높은 온도에서 일들이 진행되는 방식과 낮은 온도에서 일들이 진행되는 방식 사이에는 많은 차이가 있다. 우선 낮은 온도에서는 매우 낮은 에너지 장벽이 수반되게 된다. 

이론적으로 특정 화학반응에 대한 전이상태를 탐색하거나 낮은 에너지 장벽을 가지는 집단 원자 이동경로를 계산하는 것은 매우 주요한 관심 사항이다. 하지만, 통계적 물리학적 관점, 고전역학적 원자 이동경로 계산은 매우 복잡한 것이다. 따라서 매우 많은 컴퓨터 자원을 필요로 하기 때문에 극도로 단순화된 방법들만이 실제로 사용되어 온 것이 사실이다. 전이상태 계산을 위해서 dimer method, nudged elastic band method, string method 등의 방법이 널리 이용되고 있다.  여러 가지 다양한 원자 이동경로를 고려하지 않기 때문에 어떠한 경로가 최적의 경로인지 논하지 않는다. 최적의 경로, 가장 확률적으로 가능성이 높은 경로를 논하기 위해서는 반드시 여러 경로들을 동시에 논해야만 한다. 당연히 계산량이 많아질 수밖에 없다.  

분자동역학 방법은 철저하게 시간에 따른 순차 사건 나열이기 때문에 시간 단계별로만 계산이 진행된다. 단계와 단계사이는 약 1 fs 수준의 시간 차이를 가지고 있다. 이 보다 더 긴 시간 차이를 사용하면 단계적 계산에 의존하는 분자동역학 계산의 정밀도는 없어져 버린다. 다시 말해서, 1 sec을 전산모사하려면 10^{15} 단계의 분자 동역학 계산을 수행해야만 한다. 2.936 msec 이 역사상 최장 전산모사 기록이다. 이마저도 특수 컴퓨터 Anton (D. E. Shaw 그룹에서 개발한 분자 동역학 계산에 최적화된 컴퓨터)을 활용한 것이다. 계산 시간 최소 단위(1 fs)와 관측하고자 하는 최소 시간 단위(1 sec) 사이의 시간 스케일 차이가 너무나 현저하다. 물론, 실제 계산에서는 상호작용 포텐셜의 정밀도 또한 문제가 될 수 있다. 이러한 포텐셜의 정밀도 문제를 언급하지 않는다고 하더라도 분자 동역학 계산에서의 근본적인 물리적인 시간 스케일의 차이(1 fs 와 1 sec)는 현재의 컴퓨팅 능력으로는 극복하기 힘든 것이 사실이다. 

action은 상호작용하는 원자들이 만드는 이동경로들의 함수이다. 각각의 원자들은 자신의 초기와 말기 위치가 있다. 각 원자들은 특정한 이동경로를 취할 수 있다. 물론, 원자들은 상호작용을 한다. 이렇게 상호작용을 하면서 모든 원자들이 동시에 이동하는 경우를 생각한다. 특정한 시간이 소요될 것이다. 이러한 이동경로들의 수는 무수히 많을 수 있다. 그러한 이동경로들 중에서 action 값이 가장 작은 경로를 찾고자하는 것이다. action 값이 작을수록 그러한 경로가 실현될 확률이 높다. 그 확률은 exp{-S/(kBT)}에 비례한다.  여기서 S는 이동경로의 함수이고 Onsager-Machlup action이라고 불린다. 특별히, S는 에너지 차원을 가지고 있다. S=S(x0, x1, x2, x3, x4, x5,...., xP)로 나타낼 수 있다. 초기 원자들의 위치는,  즉, t=0 일 때, x0 이다. 말기 위치는 xP이다. 이 때 흘러간 시간은 delta * P 처럼 주어질 수 있다. 유한한 시간 간격을 delta라고 정의하면 그렇다. 이것은 마치 exp{-V/(kBT)} 처럼 볼쯔만 분포를 가지는 원자들의 집단을 연상시킨다. 물론, 상호작용하는 에너지는 V=V(x)라는 함수로 표현된다. 원자들의 위치를 x 로 표시했다. 우리의 관심사인 자유도는 X 이다.  원자들의 집단 위치를 표시하는 것이다. 브라운 운동에서처럼 유체속에 있는 분자들과의 충돌로 부터 random 운동을 하는 자유도를 가정한 것이다. 시간 순서에 따라서 자유도 X 가 바뀌어 가게 된다.
https://en.wikipedia.org/wiki/Action_(physics)

action이 최적화(극대, 극소, saddle point) 될 때의 상호작용하는 원자들의 이동경로는 정확히 뉴튼의 운동 방정식을 만족한다.  특정 위치에서 또 다른 특정 위치로, 특정 시간 동안, 원자들이 상호작용하면서 다 함께 움직일 때 action은 정의된다. action을 최적화하는 문제는 매우 도전적인 문제이다. 왜냐하면 원자들의 이동경로는 너무나도 많은 경우의 수를 내포하고 있기 때문이다. 그밖에도 몇 가지 기술적 어려움들이 있다. 우선 모든 컴퓨터 프로그램에서는 운동 방정식이 반드시 이산화 되어야만 한다. 따라서 해석적인 조건들이 만족이 안 될 수 있다. 예를 들어 시간 간격은 반드시 특정한 값이 되어야 하고 0으로 수렴되는 그러한 값이 될 수 없다. 따라서 action 은 사용하는 시간 간격 값에 따라서 그 표현과 의미가 달라질 수 있다. 이러한 수치적 어려움을 극복하기 위한 방법이 필요하다. 또한, 이산화된 action을 사용할 경우, 계산된 이동경로를 따라서 총에너지(운동에너지+포텐셜 에너지)를 계산할 경우, 이동경로를 따라서 총에너지가 보존되지 않는 약점이 있다. 즉, action 값이 적절히 국한되지 않는 문제가 있다. 또한, 얻어진 이동경로가 뉴튼의 운동 방정식 풀이의 성질을 만족하지 못하는 것이다. 이들 몇 가지 수치적 문제점들은 해결이 가능하다. 이러한 문제점들을 해결한 방법이 ADMD  (action-derived molecular dynamics)  방법이다. 주어진 원자구조에 대해서 힘과 에너지를 계산할 수 있으면 ADMD 를 수행할 수 있다. ADMD 방법은 병렬적으로 계산이 된다.

ADMD 계산은 단백질 접힘 과정과 같은 매우 복잡한 원자들의 집단 이동경로를 탐색할 때 유용하다. 다른 계산들 보다 더 일반적인 원자들의 이동경로들이 자연스럽게 탐색될 수 있다. 하지만, 여전히 몇 개의  이동경로들을 탐색할 수밖에 없다. 계산량이 많고 이동경로들 사이의 통계적 의미가 확실하지 않다. 또한, 유한 온도 전산모사가 될 수 없다. 다만, 닫혀진 시스템의 총에너지 관점에서 이동경로를 탐색하는 특징이 있다. 

주어진 온도(T)에서 시스템은 소위 볼쯔만 분포를 가지게 된다. 시스템의 에너지 V가 주어질 때, exp(-V/(kBT)) 이 그 확률이 된다. 당연히 V=V(x)로 표시된다. x 는 가능한 원자들의 배열 방식을 표시한다. 물론, 다른 물리량이 될 수 있지만 그 표현 방식에는 변화가 없다. 사실, 분자 동역학 계산은 위의 볼쯔만 분포를 가정한 것이고 시스템의 안정한 상태를 찾아내는 방법으로 볼 수 있다. 매우 다양한 배열들 ({x}) 효율적으로 찾아낼 수 있으면 좋다.   
   

다양한 이동경로에 해당하는 action 값들에 대한 ranking을 매길 수 있다면, 보다 더 물리적으로 중요한 원자들의 이동경로를 얻을 수 있다. 다시 말해서, 매우 다양한 원자들의 이동경로를 action 값 기준으로 광역 최적화(global optimization)하는 것이 가능하다. 광역 최적화된 action을 기준으로 가장 확률이 높은 원자들의 집단 이동경로를 얻어낼 수 있다.

상호작용 포텐셜 함수의 명시적인 2차 미분 없이, 명시적 1차 미분만으로 고전적 action을 계산하였고 Onsager-Machlup action을 직접 계산하여 원자들의 집단 이동경로에 대한 계산된 action의 ranking 을 매겼다. 얻어진 경로들의 변이와 교차를 포함하여 원자들의 경로를 광역 최적화 하였다. 구체적인 이동경로 계산은 Action-Derived Molecular Dynamics (ADMD) 방법을 사용하였다. 광역 최적화 방법으로서 Conformational Space Annealing (CSA) 방법을 사용하였다.   

ADMD 계산과 CSA 계산은 각각 병렬적으로 수행이 가능하다. 왜냐하면, ADMD 는  통상의 분자 동역학 계산과 달리 초기조건 문제 풀이가 아니고 경계조건 문제 풀이이기 때문이다. 경계조건 안쪽에 있는 다양한 다수의 분자모양에 대한 에너지와 힘 계산이 동시에 진행될 수 있다. 이 점을 활용하면 병렬계산이 가능하다. 아울러, CSA 계산은 독립적인 ADMD 계산들 다수를 동시에 취급할 수 있기 때문에 계산이 병렬적이다. 시도 경로 형태로 주어진 분자모양들에 대해서 동시에 에너지와 힘을 계산할 수 있다. 이렇게 계산된 에너지와 힘을 이용하면 action 을 국소 최적화시킬 수 있다.   

예를 들어, 매우 유연한 분자인 단백질 분자의 접힘 과정을 고려해 볼 수 있다. 상호작용하는 원자들이 어떠한 순서로 분자형태를 만들고 또 바꾸면서 에너지가 낮은 상태인 고유의 모양을 찾아가게 될까? 고려하고 있는 단백질에 따라서 상황은 달라진다.
https://en.wikipedia.org/wiki/Protein_folding#/media/File:Protein_folding.png



1953년 Onager-Machlup action에 관한 논문이 발표된다. 특정 온도에서, 특정 시간 동안, 한 위치에서 다른 위치로 상호작용하는 입자들이  이동하는 경로가 실현될 확률을 이론적으로 제시한다.  사실은 이 이론에 제일 충실한 방법이 있을 수 있다. Onsager-Machlup action을 직접적으로 광역 최적화 시킬 수 있다면 좋다. 하지만, 실제 계산에서 이는 상당한 컴퓨터 자원을 요구한다. 왜냐하면 이러한 계산들은 결국 상호작용 포텐셜 함수에 대한 고차 미분을 요구하기 때문이다. 예를 들어, 어떠한 함수의  국소 최적화를 위해서는 그 함수의 미분이 필요하다. 물론, 여기서 미분이 필요한 이유는 매우 효율적인 함수 최적화를 의미하기 때문이다. 미분없이 함수를 최적화 할 수 있지만, 미분을 활용하는 것과 비교하면 그 효율성은 매우 떨어진다.  따라서, 가능하다면, 상호작용 포텐셜 함수의 1차 미분 수준에서 계산을 마무리할 수 있으면 좋다. 이것은 정확하게 분자 동역학 계산에서 요구하는 수준의 것이다. 에너지와 힘을 요구하는 분자 동역학 계산에서처럼 간단하게 함수, 1차 도함수 수준에서 모든 계산을 마무리하고자 하는 것이다.

랑쥐방 방정식은 볼츠만 분포를 만들어야만 한다. 마찰항에 비해서 질량항이 무시될 수 있는 경우를 생각할 수 있다. 여전히 랑쥐방 방정식에 해당된다. 이동경로는 볼츠만 분포를 따를 것이다. https://en.wikipedia.org/wiki/Langevin_equation
질량 항이 마찰 항 보다 작은 경우, 랑쥐방 방정식은 아래와 같이 적어진다. 이러한 경우, 시스템은 에너지가 낮을수록 더 많이 발견되는데, 그 확률은 정확하게 볼츠만 분포를 따르게 된다. 이러한 상황에서 입자들은 움직이게 되고 경로를 만들게 된다. 경로의 시작 xi, 경로의 마지막 점은 xf 로 표시될 수 있다. 이 때, 소요된 시간은 t 로 나타낼 수 있다. 이렇게 한정된 경로가 발견될 확률은 P(xf|xi; t) 로 표시하고 아래와 같다.

모든 이동경로에 대한 적분을 완료하면, 특정 온도에서, 특정시간 동안에, 상호작용하는 입자들이 초기위치에서 출발하여 말기위치에 도달할 확률을 계산할 수 있다. 이 때, 확률적으로 더 중요한 입자들의 경로를 선별할 수 있다. 이것은 단순히 action(Onsager-Machlup) 값으로 정해진다. action 값이 낮을수록 더욱더 중요한 입자들의 이동경로가 된다. 왜냐하면, 그 경로가 채택될 확률이 높기 때문이다. 가장 작은 action 값을 가지는 이동경로는 확률적으로 매우 중요한 주요 이동경로가 될 수 있다.

이론적 관점에서 새로운 점:  다양한 반응경로를 동시에 취급한다. 다양한 여러 반응 경로들로부터 가장 action이 작은 값을 직접 찾아낸다. 광역 최적화 작업을 통해서 최소 작용을 직접 찾아낸다. Action-CSA 방법은 사실상 포텐셜 함수의 2차 미분 효과를 수치적으로 고려하고 있다. 주어진 온도에서 이상적인 반응경로(reaction pathway)를 찾아낸다. 이론적으로 잘 정의된 반응경로를 수치적으로 찾아낸다. 반응경로에 대한 교차, 변이를 활용한다. 통계역학적으로 가장 이상적인 원자 이동 경로 계산방법이다.

장점: 이론적으로 가장 완전한 계산 방법이다. 주어진 온도에서 가장 실현가능성이 높은 이동경로를 얻을 수 있다. 가능한 모든 종류의 포텐셜 함수에 적용이 가능하다. 에너지, 힘 계산만을 필요로 한다.
약점: 막대한 계산 자원을 요구한다. 병렬 컴퓨터가 반드시 동원되어야 한다.

참고 문헌:
https://arxiv.org/abs/1610.02652
https://www.nature.com/articles/ncomms15443

ADMD (Action-Derived Molecular Dynamics):
http://incredible.egloos.com/372174
http://incredible.egloos.com/3127566
http://incredible.egloos.com/3947906
findingtransitionpathways.pdf
RareEvent_simulations1-49.pdf
RareEvent_simulations50-132.pdf
1566_article_physics_and_high_technology.pdf
CSA (Conformational Space Annealing):
http://incredible.egloos.com/6208916
https://arxiv.org/pdf/1003.0971.pdf
http://incredible.egloos.com/7304654
http://incredible.egloos.com/4520344
라르스 온사게르:
http://incredible.egloos.com/7304654



핑백

덧글

  • csa 2017/06/12 11:26 # 삭제

    지난 주에 이주용 박사가 KIST에 세미나 다녀갔습니다. 이 주제를 가지고 발표했고, 흥미롭게 들었습니다. 볼츠만 분포와의 상관관계는 아직 미지수라고 하더군요. KMC 연구하는 사람으로 재미난 방법이었습니다.
  • 바죠 2017/06/12 11:32 #

    볼츠만 샘플링을 추구하는 것는 통상의 MC 로 취급하면 됩니다.
    그 것보다 더 중요한 가장 빈도가 높은 경로를 찾는 것이 급선무입니다. 그래서, dominant pahthway를 추구하는 것이 Action-CSA라고 할 수 있습니다.
    Action-CSA 이전에는 그냥 몇 가지 경로들을 찾는데만 만족했었습니다.
    그것이 ADMD 라는 방법입니다.
※ 로그인 사용자만 덧글을 남길 수 있습니다.

최근 포토로그