Harmonic Fourier Beads (HFB) 방법 (free energy) by 바죠


HFB (Harmonic Fourier Beads) 방법 [J. Chem. Phys. 125, 174108 (2006) ] [Ilja V. Khavrutskii, C. L. Brooks III, http://mccammon.ucsd.edu/~ikhavru/  저자 웹페이지에 manuscript가 올라와 있음. 2007년 5월2일 확인함.]은 유한 온도(T) 조건 하에서 최소 자유에너지 경로 (minimum free energy path)를 계산하는 방법이다. I. V. Khavrutskii 박사 등이 최근에 개발한 방법이다. 절대 온도 0 K로 세팅하여 계산된 경로는 최소-포텐셜 에너지 경로 (minimum energy path) 계산으로 볼 수 있다.  transition mechanism, transition rate를 이해할 수 있는 바탕으로 받아들여 진다.

반응전 상태에서 반응후 상태로  분자 모양이 변할 때, 각 원자들의 움직임을 계산하는 것을 넘어서 유한 온도 조건하에서 가능한 최소 자유에너지 경로를 계산할 수 있다는 것이다. 다양한 전이경로로 부터 얻어질 수 있는 자유에너지 계산을 가능하게 해준다는 점이 특징이다.

최소 포텐셜 에너지 경로를 계산하는 것과는 차이가 있다. 당연히, 더욱더 어려운 계산으로 보인다. 유한 온도 조건을 활용해야 하기 때문일것이다. 관련된 샘플링 때문에 계산량도 상당히 많을 것이다.HFB 방법은 병렬계산에 유리한 알고리듬이다.

유한 온도 조건하에서의 전이 경로 (transition pathway) 계산과 관련하여 많은 연구가 진행되어 왔다. HFB 방법은 Finite temperature string (FTS) 방법보다 쉽게 구현할 수 있다. hyperplane 사이의 '일'의 양를 계산하는, revirsible work 계산을 행하는, 방법보다 안정적으로 보인다. [J. Chem. Phys. 122, 114108 (2005). ]

최소 포텐셜 에너지 경로 계산 방법으로 NEB (Nudged Elastic Band), string, ADMD (Action-Derived Molecular Dynamics; 작용유도 분자동역학) 방법들이 개발되어 사용되고 있다. 이들 모두는 두 개의 분자 구조들(reactants, products)을 알 때, 전이경로 (transition pathway), 전이상태 (transition state)를 계산하는 방법으로 활용될 수 있다

당연히 기존의 분산된 분자 모양 변화들을 활용하는 방법과 밀접한 관련이 있다. 처음 인풋으로 {q_j}, 분자구조 내의 원자의 이동 경로를 줄 수 있다. 물론, 분산된 양으로 주어진다. j=0, j=1, j=2, ...., j=P.  ADMD, NEB, string 방법의 경우, j=0 과 j=P는 고정된 분자 구조를 나타낸다. 각각, 반응전 물질, 반응후 물질들을 reactants, products로 표시할 수도 있다. j 를 step index라고 부르고, 사용하였다. alpha= j Delta/ tau = j/P, Delta=tau/P 로 표시할 수도 있다. 전이경로를 연구하는 입장에서는 j=1,2,3,...,P-1에 해당하는 분자구조들을 찾아낸다는 것을 의미한다. 다시말해서, 계산에서 구해야하는 변수들의 수는 3N*(P-1) 가 됨을 의미한다. 여기서 N은 원자의 수를 나타낸다. 각 원자의 이동을 기술하는 데, 3*(P-1)개의 변수를 사용한다.

alpha는 0 에서 1까지 연속적으로 바뀔 수 있다. 여기서 alpha는 사실상 step index로 표시될 수 있는 양이기도 있지하지만, 그 이상의 의미를 나타낼 때, 사용되어질 수 있다. 바로 interpolation 에 활용될 수 있다는 점이다. sine함수들을 활용할 경우, 원자의 이동경로를 연속적인 함수로 표시할 때 사용되어 질 수 있다. 


 

http://en.wikipedia.org/wiki/Image:Activation2_updated.svg : 그림출처

HFB는 NEB에서 이용하는 spring 을 활용하지 않고(다른 방식으로 사용한다고 볼 수 있다.),  아래와 같이 원자가 특정 좌표에 머물게하는 추가적인 bias potential을 도입함으로써, 전이경로를 계산을 할 수 있게 설계된 것이다.  NEB 방법의 경우 스프링에 의하여 원자에 작용하는 힘과 분자 내부에서 자체적으로 생성된 원자에 작용하는 힘, 이들 둘 다 고려함으로써, 전이경로를 찾아가는 양식이다. 이 때, 특별히 전이 경로 방향에 수직인 성분만 활용한다는 것이 특징이다.  이 특징 때문에 NEB에 N이 있는것이다. 수직으로 찌러다의 뜻을 가진 nudged를 사용하였다. plain elastic band는 이 개념을 활용하지 않은 상태를 나타낸다. 전이상태 계산에 특별한 형식의 문제점이 생기게 된다.

아래에 표시된 것처럼 조화진동자 형식의 bias 된 포텐셜을 추가로 부여한다. 이러한  추가 포텐셜 하에서 얻어진 평균위치(<q>_b)를 계산함으로써 bias된 에너지 함수와 관련된 자유에너지 변화량을 계산할 수있다.  NEB, string 방법에서처럼 특수하게 정의된 "목적함수"는 정의되어 있지 않다. 여기서 지적한 자유에너지 변화량은 자유에너지를 step index에 대한 미분으로 연결 시킬 수 있음을 의미한다. ADMD 방법의 경우 목적함수를 정의하고 목적함수를 최소화시키므로써 전이경로를 계산하는 특징이 있다. 

연속된 분자구조들(이미지들) 사이의 특별한 관계 없이, 그야말로 독립적인 에너지 최적화 또는 MD(Molecular Dynamics; 분자동역학) 전산모사를 수행함으로써 최종적으로 전이경로를 계산할 수 있다. 이미지들 사이의 간격은 정확히 같게되도록 요구한다. 결국, 이 부분이 NEB의 스프링 역할을 한다고 볼 수 있다. 스트링 방법에서도 유사한 절차를 거친다. 이러한 작업이 없으면 모든 이미지들이 '낮은 에너지"쪽으로 쏠리게 될것이다. 이렇게 이미지들이 "쏠리는 현상"을 방지하기 위한 방법이 있다.

이러한 작업을 수행할 때, "직선+사인함수들"의 공식을 활용하여, 사실상 연속함수처럼 원자 이동경로를 활용한다. 거리계산에서 연속함수 적분하듯이 할 수 있다. 나중에 자유에너지 계산에서도 마찬가지이다. 적분할 때, 연속함수 적분하듯이 한다. 사실상 연속함수의 적분에 지나지 않는다. 그것도 1차원 적분, 심프슨 공식 을 그대로 적용할 수 있다. 


입력으로 주어진 q_j를 이용하여 분자동역학을 실행할 수 있다. 이렇게 하면, q_j에 해당하는 <q>_b를 구할 수 있다. 바이어스된 상황에서 얻어진 위치 평균값을 얻을 수 있다. q_j를 레퍼런스로 사용하고 이 레퍼런스에서 조화함수꼴로 bias 된 포텐셜 함수와 더불어 sampling을 취한다. 이러한 MD/MC 등의 방법을 활용하여  <q>_b를 얻을 수 있다. 각 이미지에 대해서 이러한 작업을 수행한다. 각 이미지에 대한 작업 수행은 병렬화되기 쉬운 형식의 작업이다.


HFB 방법과 관련하여 하나 꼭짚고 넘어가야할 점이 있다. 아래에 표시된 양으로 부터 자유에너지를 계산할 수 있다. 굉장히 간단한 공식으로 자유에너지 프로파일이 얻어질 수 있다.  처음으로 주어진 이미지들이 당연히 레퍼런스 위치이다. 그리고 추가적인 bias 포텐셜하에서, 유한온도 하에서 MD(분자동역학)로 부터 얻어진 평균 위치를 계산할 수 있다. 놀라운점은 이 양이 자유에너지를 한 번 미분한 양이라는 것이다. 자유에너지계산은 이것을 적분함으로써 얻을 수 있다는 뜻이다. [ J. Chem. Phys. 123, 144104 (2005) ] 물론, 이 자유에너지 미분은 근사적으로 얻어진 것이다. 해석적으로 얻어진 것이 아니다. 이 이론은 소위, Umbrella Integration이라고 하는 것이다. Thermodynamic Integration과는 구별되는 것이다.

이 양을 계산하여 레퍼런스의 위치를 바꿀 수 있다. HFB iteration (cycle)이라고 한다. 이 때, 다음 iteration을 위해서, 계산 중간에 동일간 이미지들 사이의 간격을 요구한다. 이렇게 하면 단순하게 에너지가 낮은쪽으로 이미지들이 이동하고, 인접한 이미지들 사이의 동일한 거리 조건에 의해서 전이경로가 자연스럽게 배열되어 나오게 된다. 여러 이미지들에서 이러한 계산들을 동시에 수행할 수 있다. 병렬계산이 용이하다.  이렇게 계산된 "힘" 양을 다시, 좌표에 적용했던것처럼,  "직선+사인함수들"로 변환하여 연속함수로 만들고 (이 번에는 <q>_b로 부터 얻어낸 힘을 연속함수로 만들었다.), \int F.dr을 수행하면 주어진 온도 조건하에서 자유에너지(일의 양) 프로파일을 얻을 수 있다.  \int F.dr은 결국, \int function(alpha) d alpha로 바뀔 수 있다. 1차원 적분이다. alapha의 함수로 자유에너지 프로파일을 얻어서 그릴 수 있다.

<q>_b는 bias potential 하에서 계산된 평균위치를 뜻한다. 이것 대신에, bias potential 하에서 에너지 최소화로 부터 얻어지는 위치, q^m_b 를 사용하면 최소에너지 경로를 따라서 에너지가 미분된 형식이라고 할 수 있다. 이 경우는 최소 에너지 경로를 계산하는 방법이된다. 사실은 T=0 조건으로 풀어도 될 것이다. 고속싸인변환을 이용할 수 있다. {q_j}---> {a_k}로 단박에 변환되는 것이 고속 싸인변환이다.

아래와 같이 사용하는것이 더 편리하다. 연속함수라서 그렇고, 원자 표시를 위해서 서브스크립트를 더 사용하지 말아야 할 듯, 아래의 수식을 보면 더욱 그렇다. "직선+사인함수들"로 표시된 원자이동 경로, 사실상 alpha의 연속함수임을 알 수 있다.



원자 이동 경로 표현, 연속함수로 표현할 수 있다. 동시에 위에서 논의한 자유에너지 변화량도 마찬가지로 alpha의 함수로 표현, 계산할 수 있다. 이때, {a_k} 대신에 {b_k}를 얻어 낼 수 있다. 이러한 식들을 이용하여 적분을 정밀하게 수행할 수 있다.
j=0, 1, 2, 3, ..., P. 

여기서, N개의 원자가 관련되어 있다. 3N 차원 벡터들의 내적을 일차원에서 적분하면 자유에너지 프로파일를 얻을 수 있다. 아래에 공식을 나타내었다.

HFB에서는{q_j}를 계속해서 바꿔어줘야한다. 물론 수렴할 때까지이다.  {q_j}를 가장 쉽게는 {<q>_b}로 대입해서 다음 단계의 HFB계산을 수행한다. 물론, 대입한 다음 '동일 이동거리' 조건을 부여한다. 이렇게 새롭게 생성된 {q_j}를 활용하여 다음 단계의 HFB를 수행한다. 다음 단계의 HFB 계산에서는 MD 계산을 수행하고, 다시 <q>_b를 계산한다. 각 HFB단계에서 A(alpha)를 계산할 수 있다. 이것이, alpha에 의존하는 자유에너지 프로파일이다. HFB단계에 따라서 자유에너지 프로파일의 수렴성을 조사함으로써 계산을 종료할 수 있다.

결국, HFB는 주어진 온도조건하에서 수렴된 {q_j}를 찾는것과 같으며, 다른 말로는 수렴된 <q>_b를 찾는것과 같다.

여기서, q는 3N 개의 원자위치 좌표를 표시하는 것이다. N 개의 원자가 있을 경우를 가정했다.

Harmonic Fourier Beads (HFB)방법을 사용할 경우 원자들의 이동경로는 "직선+사인함수들" 로 표현한다. 이 경우, 원자들의 이동경로를 사실상 연속함수로 사용함을 의미한다. 이렇게 할 경우, alpha는 0 에서 1 까지 연속적으로 바뀌는 실수이고, 이것이 반응경로를 나타내는 좌표로 생각할 수 있다. 이 때, 분자의 모양 변형에 따른 이동경로들의 길이, 거리를  매 단계, alpha마다 같게 둘 수있다. 즉, 새로운 alpha 값들을 순서대로 찾을 수 있다. 이동경로를 계산하면서, 특별한 alpha_j를 찾을 수 있다. 이 때, alpha_j는 이동경로들의 거리가 모든 단계 j에 대해서 똑같게 만들어 주는 alpha를 의미한다. 자유에너지 일차 미분에 해당하는 양도 "직선+사인함수들"로 표시할 수 있다. 이 양도 사실상 연속함수로 나타낼 수 있기 때문에, 1차적분을 수행하여 자유에너지 프로파일을 alpha의 함수로 계산할 수 있다.

이것은 결국, 일차원 해를 찾는 문제이다. zeroin 라이버러리 함수를 사용하여 해를 구할 수 있다. 이렇게 alpha_j를 구하면, 분자모형 변화의 단계를 '같은 이동거리'로 표현할 수 있다. 

이렇게 할 경우, HFB 방법에서는 이미지들 사이의 상관관계를 '같은 이동거리'라는 개념하에 인위적으로 도입할 수 있다. NEB의 경우, 스프링을 활용하여 이미지 사이의 관계를 도입했었다. 이미지들이 에너지가 낮은쪽으로 흘러 내려가는 것을 자연스럽게 방지할 수 있다. HFB 방법에서는 목적함수는 없다. 굳이 목적함수와 유사한 양을 선택하라고 하면은 bias potential과 기존의 포텐셜 함수의 함을 일시적인 목적함수라고 볼 수 있다. 일시적인 목적함수는 국소적으로 정의된다. 레퍼런스 위치에 의존하기 때문에 그렇게 볼 수 있다. 물론, 동등한 거리 조건은 이러한 일시적이고 국소적인 목적함수 최소화 뒤에 항상 따라오게 되어 있다.

주어진 온도 조건하에서 임의적으로 evolution하는 bead들이 결국에는 '같은 이동거리'조건에 의해서 전체적으로 재배치될 수밖에 없다. 이미지들을 반응경로 상에서 배열하는 방법이 HFB와 NEB에서 구조적으로 다름을 알 수 있다. 

분산된 이미지 사이들로 부터 추가로 유도되고 이어지는 계산에 사용되는 수학적 양이없다. 이런 양들은 정밀도에서 도움을 안 줄 수 있기도하다. 즉, hyperplane, normal vector같은 것들이 없다. forward, backward 의 방향성이 전혀 없다. 완전히 독립된 계산을 수행한 후, 동일 이동거리 조건으로 이미지들을 다시 배열한다.

HFB 방법의 경우, 이미지들 사이의 커플링을 요구하지 않는다. 하지만, 계산 단계에서 근접한 이미지 사이의 동일한 이동 거리를 추가로 요구한다. 이러한 단계의 계산을 통해서 이미지 사이의 커플링이 자연스럽게 들어오게 된다. 

NEB 방법이 유도된 이유중의 하나는, 동시에 NEB 방법과 유사한 방법들이 계속 발표되는 이유중 하나는, 바로 이미지들 사이의 간격 유지와 관련된 것이 많이 있다. 사실, 이 부분이 여러 이미지들을 활용하는 방법들에서의 핵심 사항이다. string 방법도 이와 유사한 처리가 있음을 알 수 있다.

레퍼런스 q_j, 온도 T 를 입력으로 하여 각 MD 계산을 통해서 수렴시켜야 할 양은 <q>_b 이다. 즉, bias potential 하에서, 주어진 온도 조건 하에서 허용된 원자의 위치. 이 양만 정확히 계산되면 현재 단계의 HFB 계산은 정밀하게 된것이라고 볼 수 있다.  물론, HFB 계산의 입력으로 {q_j}가 잘 주어져야 한다. 이때, 레퍼런스 이미지들을 준비함에 있어서, 인접 이미지간의 이동거리가 같게 주어져야 한다.


자유에너지 프로파일 계산 방식에는 여러 가지 방법들이 현재 사용되고 있다. 대표적인 계산 방법들은 reaction coordinates를 미리 정의하고, 즉 선택하고 이에 대한 활률적 접근법이다. 이러한 reaction coordinates 없이 미리 정의하지 않고 분산된 이미지들을 입력으로 해서 자유에너지 프로파일을 계산하는 방법은 앞서 논의한 방법들 보다는 일반적인 자유에너지 프로파일 계산 방법이라고 분류할 수 있다. 이러한 의미에서 HFB 방법은 그 자체의 높은 병렬효율성을 포함하여 평가할 때, 상당히 일반적인 자유에너지 계산 방법의 하나로 분류될 수 있다.

 

위의 그림에서는, HFB 방법을 이용하여 자유에너지 경로를 계산할 때, 순서도를 표시했다. 

물리적으로 알아본 HFB iteration(cycle)의 특성은 아래 그림에서 잘 설명되어 질 수 있다. 각 단계를 통하여 MD 계산이 수행된다. 이 계산에서의 입력은 레퍼런스 q_j와 온도 T이다. 이때, 레퍼런스들은 동등한 이동거리 조건으로 부터 유도된 레퍼런스 이미지들이다. 이 계산을 통하여 <q>_b가 계산된다. 이러한 과정을 magenta 화살표 선으로 표시했다. 각 이미지들은, 대략, 에너지가 낮아지는 쪽으로 몰리게 될 것이다. 인접한 이미지들 사이의 거리는 편차를 가지게 될것이다. 이것을 인위적으로 보정해 주는 작업을 수행한다. 이 작업은 HFB iteration이 들어가기 전에 해준다. 새로운 HFB iteration은 통상 <q>_b를 그대로 레퍼런스로 둔다고 생각하면 된다. 즉, {<q>_b}--->{q_j}.   {<q>_b}가 최소 자유에너지경로에 해당한다. 이에 대응하는 자유에너지 프로파일도 계산할 수 있다. 굉장히 간단한 알고리듬이다. 병렬계산에 적합한 알고리듬이다. HFB 방법에는 목적함수가 정의되어 있지 않다.




먼저, 아래에 나타낸 공식에서 처럼, 거리단위를 0 에서 1 까지 적분을 하면 우리는 총이동거리 "L" 을 얻어 낼수 있다. 우리가 사용하는 간격의 수는 P이기 때문에 동일한 거리로 L/P를 선택할 수 있다.
인접한 이미지들 사이의 거리가 똑같게 (L/P) 되도록 만들어주는 해, {alpha_j}를 결정할 수 있다. 아래의 순서처럼, {alpha_j}를 결정할 수 있다. 제일 먼저, 0 서 부터 적분을 해서 거리가 L/P가 되는 alpha_1을 결정한다. 이 때, zeroin함수를 사용한다. 그 다음, alpha_1 에서 적분을 시작해서 거리가 L/P가 되는 두 번째 구간을 설정해 주는 alpha_2를 찾아낸다.


아래의 그림은 하나의 예이다. 처음에는 구간별로 거리가 다르다. 위에서 언급한 순서대로 적절한 {alpha_j}를 순서대로 결정하고 나면, 동일한 인접 이미지 사이의 간격을 확보할 수 있다.

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

NEB 방법:


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

ADMD 방법:


FTS(finite temperature string) 방법: constrained dynamics, mean string--> centerline of the tube, mean force sampling on the final plane--> free energy, reaction coordinate is not known beforehand.

a collection of paths

transition tubes---> MEP at T=0, as T-->0.

transition pathway sampling: real  physical time,
J. Phys. Chem. B 109, 6688 (2005).

NEB, zero-temperature sting방법

conjugate peak refinement, dimer

cf.
Dr. W. Quapp, http://www.math.uni-leipzig.de/~quapp/

cf.
이전에 관심을 가지고 있었던 방법들:

http://incredible.egloos.com/2342691

작용유도 분자동역학(ADMD):
http://incredible.egloos.com/372174

cf.doltsinis2_free_energy_and_rare_events_in_molecular_dynamics.pdf


 



cf. HFB method
JChemPhys_125_174108_HFB.pdf
cf.
ggaHFB
JChemPhys_127_124901_gaHFB.pdf

cf. string method, collective variables
JChemPhys_125_024106_string_collective_variables.pdf

CPL_446_182_on_the_fly_STRING.pdf

핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax