전이경로를 찾아서 by 바죠

 
 

전이경로를 찾아서


 

서론

  물질들의 근본 특성들을 이해하기 위해서는 미세구조 분자세계에서 표현되는 원자들의 배열, 전자와 핵 간의 상호작용, 전자들 간의 상호작용을 이해해야 한다. 원자들의 움직임과 전자들의 간의 상호작용, 원자 내부의 핵과 전자들이 만들어 내는 상호작용으로 많은 기계적, 열역학적, 전기적, 광학적 물성들을 설명할 수 있다. 양자역학을 적용한 분자세계의 이론적, 실험적 기술은 매우 성공적으로 수행되어 왔다. 현재는 물질의 물성을 이해하는 차원을 넘어 새로운 물질을 만들고 디자인하는 단계에 이르고 있다고 평가할 수 있다. 이러한 의미에서 분자세계에서 일어나는 원자들의 운동과 분자들의 운동을 이해한다는 것은 물리, 화학, 생물, 재료과학 분야를 포함하는 포괄적인 과학기술 분야에서 매우 중요한 연구주제가 될 수 있다.

  최근 나노물질들을 연구하는 연구자들이 늘고 있다. 대량생산과 정보화 사회의 수요에 부응한다는 논리를 따르면 당분간 이러한 추세는 계속될 것 이다. 나노물질들을 만들어 내고 그 성질들을 예측하거나 확인하고 응용하는 것이 관련 과학자들의 연구주제이다. 나노물질들은 이미 자연계에 존재하고 있는 물질들을 포함한다. 물론, 이전에 존재하지 않았을 가능성이 있는 물질들을 포함한다. 경험적으로 보았을 때, 나노물질들은 새로운 가능성과 위험성을 동시에 사람들에게 전해주고 있다고 보는 것이 일반적인 관점일 것이다. 나노물질들의 합성, 물성 등을 완전히 이해하기 위해서는 다양한 실험 방법과 이론적인 이해가 뒷받침되어야 할 것이다. 많은 과학자들은 분자수준에서 일어나는 원자배열들의 변환과 관련된 물질의 특성, 원자 수준의 기작들에 대한 구체적인 정보에 목말라하고 있다.

  포괄적 의미로 원자들의 재배열 또는 화학반응에 상존하는 원자들의 움직임들을 신뢰할 수 있는 이론적 체계 내에서 설명하는 것은 현실적으로 많은 과학기술자들이 얻고 싶어 하는 중요한 정보 중의 하나이다. 기능성 나노물질의 합성, 화학반응의 이해, 조절 등과 맞물려 원자들의 이동경로 그 자체는 매우 중요한 연구주제가 되어왔다. 관심 있는 반응이 무엇인가에 따라서 반응에서 얻고 싶은 정보의 수준은 달라질 수 있을 것이다. 원자들은 어떻게 움직임으로써 하나의 화학반응을 이루어낼까? 또는 어떻게 생명활동의 단위과정을 원자들의 재배열로 설명할 수 있을까? 그 자체로서도 충분히 연구해 볼 가치가 있는 연구주제가 아닐 수 없다. 원시 화학기술의 형식을 취한 연금술 시대보다는 분명히 현 시대의 과학적 수준이 높아졌다고 말할 수 있다. 하지만, 여전히 과학자들은 원자의 운동을 이해하고 제어하기 위한 궁극의 목표에 도달하지 못하고 있는 상황이다. 현대 과학자들도 이전 시대와 마찬가지로 과학자의 기본적 상상력을 실제적, 기술적 요소들에 혼합할 수밖에 없다.

  본 논설에서는 최근 화학반응 또는 원자 재배열과 같은 미시세계 수준에서 원자들의 움직임을 결정하는 계산방법론 개발의 의미와 전산모사 방법의 아이디어들, 그리고 관련된 응용연구들을 개괄적으로 살펴보기로 한다. 최근에 Passerone, Parrinello 박사들에 의해서 제안된 작용유도 분자동역학(action-derived molecular dynamics; ADMD) 방법의 이론, 해당 이론의 확장, 다른 방법들과의 비교, 그리고 응용들에 대해서 자세히 기술하기로 한다.


전이경로란?

  전이경로(transition pathway)란 일반적인 화학반응과 같은 분자구조(원자배열) 변화에 있어서 원자들의 집단적인 이동경로들을 포괄적으로 표현하는 말하는 것이다. 이로부터 자연스럽게 분자구조 변화의 단계들을 살펴 볼 수 있을 것이다. 일반적으로 화학반응과 관련된 전이경로는 하나의 경로를 가지지 않고 그 자체로 다양성을 내포하고 있다. 단백질이 하나의 생명현상과 관련된 고유의 단위기능을 수행하는 경우를 생각할 수 있다. 주변의 원자, 분자들과 상호작용하면서 특정한 분자가 한 상태에서 또 다른 분자상태로 변화하는 과정을 연속적으로 상상할 수 있다. 이러한 연속적인 분자모양의 변화들도 전이경로의 한 형태로 볼 수 있다. 같은 화학반응이라고 하더라도 서로 다른 전이경로들도 생각해 볼 수 있다. 이러한 전이경로들의 모음을 전이경로 앙상블이라고 한다.

  잘 알려진 것처럼 화학반응의 형식 또는 원자 재배열의 형식은 매우 다양하다. 일반적 상식으로 분자의 모양을 바꾸는 현상, 원자들의 집단적인 움직임, 원자들 간의 결합과 분리 등을 상상하는 데는 우리가 어느 정도 동의할 수 있으나 구체적인 원자이동의 정보를 우리가 알아내는 것은 무척 어려운 일이다. 원자들의 재배열, 상전이, 좀 더 구체적으로는 이성질화, 화합, 분해, 치환, 복분해 등으로 분류할 수 있는 화학반응들을 일반적인 하나의 계산방법으로 기술하는 것도 힘든 것이 사실이다. 촉매반응의 경우, 촉매가 원자 수준에서 어떠한 역할을 수행함으로써 전체 반응의 활성화 에너지(activation energy)를 낮출 수 있는지는, 원자 수준에서 전이경로를 밝혀내지 않는 한 구체적인 촉매반응의 정보를 얻기는 힘들다.

그림 1)_ 일반적인 화학반응에서의 시간 스케일과 반응좌표와의 관계를 도식적으로 나타내었다. 일반으로 화학반응(원자 재배열)과 연관된 반응좌표의 선택에 관한 일반적인 기준은 없다. 반응좌표는 화학반응 전후뿐만 아니라 반응의 단계를 잘 기술할 수 있어야 한다. 하지만, 시간 스케일 상으로 주의를 요구한다. 실제 화학반응은 매우 짧은 순간 돌발적으로 일어나는 경우가 많다. 일반으로 아주 긴 시간 동안에 점진적으로 반응이 진행되지 않는다.

반응좌표, 전이상태

  반응좌표(reaction coordinate)라는 것은 화학반응과 연관하여 개념적으로 선택한 1차원 또는 다차원의 좌표계를 말하는 것으로 해당 화학반응의 단계를 잘 기술할 수 있어야 한다. 많은 경우 반응좌표의 선택은 연구자의 명상적, 사변적 선택에 의존할 수밖에 없다. [그림1]에서와 같이 반응전 상태(reactants)와 반응후의 상태(products)를 잘 기술함은 물론이고 반응의 중간 단계들도 잘 기술할 수 있어야 좋은 반응좌표라고 할 수 있다. 전이경로를 모르고 있는 상황에서 반응좌표를 선택한다는 것은 해당 화학반응이 매우 간단한 경우가 아니면 상당히 어려운 문제이다. 반응좌표를 설정하고 그 반응좌표를 따라서 분자의 모양을 결정하고 해당 에너지를 계산하여 화학반응을 기술할 수 있다. 하지만, 반응좌표가 적절하게 도입되지 못할 경우 해당 화학반응을 제대로 기술할 수 없음은 명백하다. 단백질의 접힘 과정에서와 같이 많은 원자들이 집단적으로 움직이면서 화학반응(원자 재배열)이 진행된다면 해당 반응을 잘 기술하는 반응좌표를 설정하는 것은 매우 어려운 작업이 될 수 있다. 간단한 화학반응의 경우, 원자간 결합 길이, 원자들 사이의 결합각도 등이 반응좌표로 사용된다.

  집단적인 원자들의 움직임과 관련된 시간 스케일은 매우 다양하다. 피코 초에서 몇 초대에 걸쳐있는 다양한 화학반응 시간 스케일들의 예를 찾는 것은 전혀 어려운 일이 아니다. 물질이 빛을 흡수하고 원자구조를 바꾸는 과정은 대략 피코초 정도에 걸쳐서 일어나는 사건이다. DNA 헬릭스가 풀리는 시간은 대략 마이크로 초에 해당한다. 효소-촉매반응의 경우 밀리 초, 단백질 합성의 경우 수 초의 시간이 소요될 수도 있다. 일반으로 화학반응에서 반응이 일어나는 과정은 시간상에서 비정규적인 분포를 가지며 순간적으로 일어난다. 대부분의 시간에는 반응이 일어나기를 기다리는 상황이 전개된다. 이것이 화학반응을 희귀사건(rare event)으로 볼 수 있는 이유이다.

   분자내의 전자들도 원자들의 움직임과 일반으로 관계를 가지고 있으나 많은 경우 전자들의 움직임을 원자들의 움직임과 분리해서 해석하는 경우가 많다. 많은 원자들의 위치를 독립변수들로 하여 정의되는 포텐셜 에너지 함수는 전자구조 연구 분야에서는 본-오펜하이머(Born-Oppenheimer) 표면이라고 부른다. 즉, 전자들의 자유도를 양자역학적으로 최적화 한 후 생성되는 다차원 공간에서의 에너지 함수 표면을 수학적으로 살펴 볼 수 있다. 여기서 양자역학적으로 전자들의 자유도를 최적화 했다는 것은 전자들에 대한 슈뢰딩거(Schrödinger) 방정식을 특정 근사 내에서 풀어서 결정했다는 것을 의미한다.

  전이상태(transition state)는 화학반응 중간의 특정한 분자구조를 말하는 것이다. 이러한 전이상태에서는 모든 원자에 작용하는 힘이 없으면서 모든 원자들의 위치에 대한 에너지 표면의 이차미분들로부터 만들어지는 동적행렬이 하나의 음수 고유값을 가진다. 원자들의 위치를 독립변수들로 하는 포텐셜 에너지 함수가 안장점(saddle point)을 만드는 상황에 해당하는 특별한 원자배열을 전이상태라고 한다. 전이상태에서 얻어진 음수 고유값에 해당하는 고유벡터는 분자구조가 변형해가는 방향을 의미한다고 볼 수 있다. 전이상태에서 적절한 속도 분포를 원자들에 할당함으로써 반응전 또는 반응후 상태로 도달하게 하는 분자동역학 반응을 수행할 수 있다. 마찬가지로, 에너지 함수의 기울기만을 따라서 천천히 이동하여 반응전 그리고 반응후의 상태에 도달할 수도 있다(intrinsic reaction coordinate method). 하지만, 일반적인 화학반응에서는 다단계의 안장점을 포함할 수 있어 단일 안장점 모델로 하나의 화학반응을 설명하기가 부족한 경우들도 많이 있다.

  간단한 분자동역학 연구에 따르면 실제 화학반응에서는 다양한 경로를 통하여 반응이 진행될 수 있고, 수학적으로 정의된 전이상태를 꼭 지나지 않을 수 있다는 것이 밝혀져 있다. 에너지적으로는 전이상태가 가장 유리하지만, 실제 반응에서는 전이상태 근처를 지나가는 경우가 많다는 것을 의미한다. [그림2]에서는 반응전, 반응후 상태들을 연결하는 가능한 전이경로들을 포텐셜 에너지 함수의 단면과 비교해서 도식적으로 나타내었다. 실질적인 동역학적 전이경로 모델에서는 다양한 전이경로들이 존재할 수 있다. 또한, 모든 전이경로들이 수학적으로 정의된 전이상태를 정확히 지나가는 것은 아니다. 실제 실험적으로 관측되는 화학 반응율은 자유에너지 프로파일을 얻어야 계산할 수 있는 물리량이다. 이 경우 다양한 전이경로들의 존재에 의한 엔트로피 효과를 고려해야 함을 의미한다. 단일 전이경로 모델로는 화학 반응율를 계산할 수 없다.


그림 2)_ 반응전 물질들(reactants), 반응후 물질들(products) 사이의 전이경로들을 도식적으로 표현했다. 반응전후의 물질들은 에너지적으로 낮은 값들을 유지하고 상대적으로 안정된 상태를 유지하고 있다. 확률적으로 가능성이 높은 전이경로는 낮은 에너지 골짜기를 따르는 반응경로들이다. 뒤에서 논의될 Parallel replica, transition pathway sampling, hyperdynamics 방법 등에서는 반응전, 반응후 상태를 다소 광범위하게 정의한다. NEB, string, ADMD 방법에서는 각각 하나의 반응전, 반응후의 원자배열 상태들을 활용한다. Dimer, ART방법에서는 반응전의 원자배열 상태만 알려져 있을 때, 근처에 있는 하나의 전이상태를 찾는 방법들이다.

실험적 한계 

  실험적 전이경로 관측의 어려움에 앞서 대표적인 희귀반응들을 살펴본다. 실험적으로 전이경로 관측에 대한 노력들을 살펴보고자 한다. 세포 내에서 단백질은 DNA 속에 보존되어 있는 유전정보들로부터 얻어진다. 단백질 생합성은 리보솜 위에서 행해진다. 실제 단백질의 합성과 기능은 아미노산 서열에 의해서 결정되어 지고 있다고 알려져 있다. 합성된 단백질은 그 고유의 기능을 펼치기 위해서 단백질 고유의 3차원 구조를 가져야 한다. 이 과정을 단백질 접힘 과정이라고 한다. 생명현상을 관할하는 단백질의 구조 변형(접힘, 잘못 접힘), 고유기능 수행 상태를 연구한다는 것은 분자세계에서 일어나는 전이경로 연구의 전형적인 한 예가 될 수 있을 것이다. 단백질이 잘못 접히는 현상은 여러 가지 질병과 연관되어 있다고 알려져 있다. 단백질 엔지니어링에서는 의도적으로 단백질의 서열상의 변이를 유도할 수 있다. 단백질 응용연구의 폭은 굉장히 넓다. 이는 분자생물학의 궁극적 목표 중의 하나에 맞닿아 있는 문제로 분류될 수 있다.

  또 다른 하나의 희귀사건의 예는 나노물질의 아이콘들인 탄소60/탄소나노튜브들이 형성되는 과정 등을 들 수 있겠다. 이와 관련된 전이경로를 밝히는 것은 실험이나 이론계산 두 분야에서 공히 결코 만만한 작업이 아니다. 실험의 경우, 원자를 볼 수 없을 뿐만 아니라 빨리 움직이는 원자들의 이동을 관찰할 수 있는 시간 분해능을 가지고 있지 못하다. 최고 분해능을 가지고 있는 고배율 투과전자현미경(TEM)을 활용해도 분자구조는 어렴풋이 확인 할 수 있을지 몰라도 원자를 구별해 낼 수는 없다. 시간 분해능 역시 문제이다. 매우 특별한 시스템이 아닐 경우 msec 이하의 시간 분해능을 이용하여 분자들의 구조변화를 관측하기란 대단히 힘들다. 놀랄 일도 아닌 것이 실제 분자 내의 원자들은 대략 초당 1012 번 정도 진동한다.

  1999년 노벨 화학상은 즈웨일 교수에게 주어졌다. 그의 연구업적은 초고속 레이저를 활용하여 화학반응의 결정적 순간을 조사하는 방법론을 개발한 것 이였다. 사실 많은 과학자들은 화학반응의 과정을 실험적으로 확인하고 싶어 했었다. 예로 들면, 화학 반응의 중간에 각 원자들은 어떠한 상태에 놓이게 될까? 상당히 복잡할 수도 있다. 어떠한 순서로 원자간 결합이 끊어지고 또한 생성될까? 나아가 촉매물질은 무슨 작용을 하고 있을까? 아무튼, 화학반응의 순간은 매우 재미있는 아주 특별한 순간임에 틀림없다.

  실험실에서 첫 번째 펨토초 펄스로 분자를 여기(excite)시켜 반응을 시작시키고, 정해진 시간이 흐른 후에 두 번째 펄스로 변화를 관찰하는(probe), 소위 펌프-프로브 (pump-probe)법을 이용하여 실험을 진행한다. 두 펄스 모두 펨토초 대의 시간 폭을 가지므로 두 펄스 사이의 경로차를 이용해 펨토초, 피코초 대의 짧은 시간 간격을 주는 것이 가능하다. 시간 간격을 다양하게 변화시키면서 얻은 프로브 신호를 종합하면, 펌프 펄스에 의해 시작된 분자 운동을 실시간으로 관찰하는 것과 같다. 즈웨일 박사팀은 한 분자의 분해반응을 연구했다. 이들은 처음으로 원자 수준에서 원자 간의 결합이 끊어지는 순간을 봤다고 봐도 무방하겠다. 하지만, 현재 이 방법을 활용하여 일반적인 화학반응을 관찰할 수 있을 정도의 수준은 아니다.

  핵자기공명(NMR) 분광학은 주로 작은 분자의 구조 분석에 이용되었다. 현재는 작은 단백질 또는 단백질 일부 부위의 구조 연구에 사용되고 있다. 소량의 농축된 순수 단백질 용액만으로도 구조 결정이 가능하다는 특징이 있다. 용액을 강한 자기장 속에 노출시키면 분자 속에 있는 수소 핵으로부터 나온 신호를 분석 하여 상호작용하는 수소 원자 간의 거리를 알아내는 것이 NRM 분광학 아이디어의 핵심이다. NMR 분광학을 활용하여 단백질 분자 속에서 형성되는 H-H 사이의 거리를 간접적으로 측정할 수 있다. 뷔트리히 박사는 NMR 분광학의 개발 공로로 2002년 노벨 화학상을 수상했다. 하지만, NMR 분광학이 거대 단백질의 분석에는 적용되기에는 어려움이 있는 것이 사실이다. 하지만, 관련된 기술들은 계속해서 발전하고 있다.


이론적 한계

  시간에 따른 원자들의 움직임을 기술하는 계산방법으로는 분자동역학(molecular dynamics) 방법이 널리 사용된다. 통상의 분자동역학 방법은 주어진 힘장(force field) 하에서 각 원자들의 운동방정식을 풀어서 전체 시스템을 시간의 함수로 기술하는 방법을 지칭한다. 다체계(many-body system)인 원자들의 운동을 일반적으로 기술할 수 있는 방법으로서 전산장비의 발달과 유용한 전산 알고리듬의 개발들에 힘입어 성공적으로 응용되어왔으며 보편화된 연구방법으로 알려져 있다. 원자들이 뭉쳐있는 분자들의 움직임을 운동 방정식으로 풀기 위해서는 분자가 가지고 있는 에너지 그리고 각 원자에 작용하는 힘을 계산해야 한다. 수소원자의 경우처럼 터넬링 운동을 하지 않는 한 원자들의 운동을 고전적으로 기술하는 것은 매우 좋은 근사임이 알려져 있다. 

  컴퓨터를 이용할 경우 원자 이동경로는 반드시 분산된 표현을 사용해야 한다. 분자동역학 방법은 분자생물, 재료물리 등의 분야를 비롯하여 일반적인 응집물질들의 운동학적, 열역학적 성질을 모델링 하는데 널리 이용되었다. 응용의 폭을 제한하는 근본 요소 중 하나는 전산모사 가능한 시간 스케일과 관련된 것이다. 현재 최장 시간동안 전산모사 된 예는 수천 ns 정도이다. 이는 막을 통과하는 이온 수송시간(마이크로 초), 단백질 접힘 시간(밀리 초)에 비하면 여전히 짧은 시간이다.

  원자수가 많은 경우에 화학반응과 관련된 전이경로를 계산하는 문제가 많은 연구자들에게 관심을 끄는 연구주제인 이유는 다음의 예를 통해서 쉽게 짐작할 수 있다. 11 fs 정도가 C-H 사이의 진동 주기인 것을 고려하면 수소 원자의 운동을 제대로 기술하기 위해서는 1 fs 수준의 이론적 시간 분해능을 확보해야한다. 그렇지 않으면 계산된 원자 이동경로의 정밀도를 유지할 수 없을 것이다. 1초를 전산모사 하려면 1015 단계가 필요하다는 뜻이다. 분자동역학 계산 한 단계가 컴퓨터 상에서 1초에 계산이 완료된다고 가정하면 전산모사 시간은 317.0979x105 년이 된다. 결국, 화학반응과 같이 정규적이지 않으며 희귀하게 일어나는 사건을 기술하기에는 분자동역학 방법이 효율적이지 못하다. 분자동역학 방법을 막연하게 생각하고 ‘전산모사 하면 되질 않느냐?’라고 쉽게 이야기 할 수 있는 입장이 결코 못 된다. 컴퓨터의 성능이 105배 좋아진다고 해도 여전히 어려움은 계속된다.

  원자들이 집단적으로 움직이면서 하나의 반응을 일으키기 위해서는 넘어야 할 에너지 장벽이 일반적으로 도사리고 있는데, 이 에너지 장벽이 높을수록 주어진 온도 하에서 전체 반응이 일어날 확률은 지수적으로 감소한다. 실온에서 관심이 있는 화학반응의 활성화 에너지 장벽이 대략 10 kcal/mol(=0.4336 eV) 정도라고 가정할 경우, 분자동역학 방법으로는 10-6 sec 이상의 전산모사를 통해서도 전이경로를 확인할 가능성은 여전히 낮다. 온도를 의도적으로 높이는 전산모사가 가능할 수 있다. 이렇게 하면 보다 높은 확률로 전이경로를 얻을 수 있다. 하지만, 높은 온도에서 얻은 전이경로는 낮은 온도에서 원자들이 취하는 전이경로와 다를 수 있다.

   분자동역학의 한 단계에서는 분자의 에너지 그리고 각 원자에 작용하는 힘을 계산하는 과정에 가장 많은 컴퓨터 자원을 사용한다. 분자동역학 방법의 병렬화는 대개 이 부분에 할애되기 마련이다. 분자동역학 방법은 시간적으로 완전히 순차적이다. 현실적으로 분자동역학 방법이 성공적으로 화학반응을 기술하기 위해서는 매우 많은 단계의 전산모사를 짧은 시간 안에 얻을 수 있어야한다고 말 할 수 있다. 반응좌표 없이, 가설을 활용한 반응상의 특별한 순서의 설정 없이 화학반응을 연구할 수 있으면 더욱 일반적인 연구방법이 될 수 있다. 아주 낮은 활성화 에너지를 가지는 몇 개의 경로들이라도 찾을 수 있다면, 화학반응을 정성적으로 기술할 수 있는 상당히 유용한 정보를 제공할 수 있을 것이다. 시간과 공간에 대해서 확장된 다중스케일(multiscale) 모델링 방법이 활발히 연구되고 있는 이유이다.

  통계 물리학적으로 각 원자들은 주어진 온도 하에서 전체적으로 특별한 분포를 가지는 운동 에너지들을 가지고 있다. 이러한 운동 에너지 분포는 분포 자체로부터 얻을 수 있는 평균 속력보다 매우 빠른 원자들과 느린 원자들의 운동을 포함할 수 있다. 또한 그 운동 방향도 상당히 특별할 수 있다. ‘이 특별한’ 운동들이 소위 반응을 불러일으킬 수 있다. 따지고 보면 분자동역학 방법은 초기조건 문제이다. 초기의 원자들의 위치와 속도를 알 때, 원자들의 위치를 시간의 함수로 계산하는 것이다. 냉정하게 이야기 하면, 초기조건으로서 원자들의 속도를 우리는 알지 못한다. 설사 안다고 한들 그것은 정말 무수히 많은 속도들 무리 중에서 하나의 속도들에 지나지 않는다. 우리는 주어진 온도 하에서 고전적으로 원자들에게 허용된 분포를 알 뿐이다. 무수히 많은 속도 분포들에 대한 분자동역학 전산모사는 사실상 불가능하다.

  완전한 화학반응의 이해를 위해서는 많은 수의 전이경로들이 필요하다. 특히, 전이경로들 중에서 상대적으로 포텐셜 에너지 함수값 상으로 높은 곳 주변에서 평평한 곳들이 매우 넓은 영역에 걸쳐 있을 경우 엔트로피 효과가 크게 나타날 수 있다. 이러한 경우, 최소-에너지-경로(minimum-energy path)로부터 구한 활성화 에너지보다 낮은 자유에너지 장벽이 나타날 수 있다. 궁극적으로는 최소 자유에너지 경로 계산이 필요한 이유가 여기에 있다.

  분자동역학 방법의 응용에서 또 하나의 어려운 점은 해당 반응이 일어난 후를 전산모사 중에 판단하여 전산모사를 기록해야 하는 점이다. 이 때, 화학 반응이 일어났다는 것을 판단하는데 사용되는 인공적인 하나의 조건이 필요하다. 사실, 원하는 희귀사건이 꼭 일어난다는 보장도 없다. 반응완료 확인 조건은 반응좌표를 선택하는 것보다는 용이한 기준이나 여전히 일반적 선택의 기준이 없는 것이 현실이다.

  지금까지 지적한 분자동역학 방법의 한계는 지금의 컴퓨터 성능보다 1000 배 뛰어난 컴퓨터가 있다고 하더라도 마찬가지 문제로 비추어질 수 있다. 이 문제가 컴퓨터의 계산능력과 상관이 있는 것은 분명하나, 단순히 그것만의 문제가 아님을 쉽게 알 수 있다. 실험적으로는 너무나도 많은 전이경로들이 가능하다는 것을 염두에 두고 있어야한다. 이론적 모델링에 있어서 확실히 새로운 아이디어가 필요한 분야이다. 결국, 새로운 이론적인 접근 방식이 필요하며, 몇 가지 보편타당한 가정들도 도입될 수 있겠다.


다양한 이론적 접근들

  실제로 많은 연구자들이 희귀사건 연구에 적용할 수 있은 다양한 계산 방법들을 개발하고 있다. 또한 관련된 응용 연구를 통해서 과거에는 다룰 수 없었던 희귀사건과 관련된 물리적, 화학적 해석이 유효함을 많은 연구자들이 보여 주고 있다. 분자동역학 방법이 기본단위가 되기도 하고, 다양한 수학적 방법들을 동원하여 희귀사건을 기술하기도 한다. 전이경로/전이상태를 결정하는 계산방법들은 단일 고정단 그리고 이중 고정단 형식의 두 부류로 크게 나눌 수 있다. 단일 고정단 형식의 계산 방법으로는 Dimer, ART, hyperdynamics, parallel replica 방법, 온도-가속화된 분자동역학 등의 방법을 들 수 있다. 이중 고정단 형식의 계산 방법으로는 nudged elastic band(NEB) 방법, string 방법, ADMD 방법, step-slide 방법, transition-path sampling 방법 등을 들 수 있다. 여기서 언급된 방법들은 포텐셜 에너지 함수의 이차미분을 활용하지 않는 방법들이다. 포텐셜 에너지 함수의 이차미분을 활용하는 방법들도 많이 개발되어 사용되고 있다. 포텐셜 에너지 함수의 이차미분을 활용할 경우 거대 분자계에 적용하기가 어렵다는 단점이 있으나 소규모 분자의 전이상태 결정에는 매우 유용한 방법임에 틀림없다.

  연구자들이 이론적으로 접근하는 방식은 매우 다양하고 그것들 자체가 흥미로운 아이디어들로 장식되어 있다. 가장 단순한 양식으로서 고려해볼 가치가 있는 것 중 하나는 동등한, 그렇지만, 초기조건이 다른 분자동역학 계산을 독립적으로 굉장히 많이 수행하는 방법이다(Pande 교수 그룹; Folding@home). 예를 들어, 약 1000 개의 독립적인 분자동역학 전산모사를 동시에 수행한다. 이를 parallel replica(병렬 복제) 방법이라고 한다. 물론, 이러한 응용연구에서는 활성화 에너지가 매우 크지 않은 경우에만 적용할 수 있을 것이다. 화학반응이 확률적으로 발생한다는 가정을 사용할 수 있을 것이다. 같은 종류의 많은 단백질 분자들 중에서 접힘 과정을 분석하고 있을 때를 가정한다. 독립적으로 한 단백질이 접히는 확률은 시간에 대하여 지수 함수적으로 변하는 분포를 가진다. 지수분포의 표준 편차는 단위시간 당 사건의 발생 횟수(사건 발생 평균 빈도)와 같다. 각각의 전산모사를 10 ns 수행했을 때, 1000개의 전산모사 중에서 4개의 전산모사에서 단백질이 접혔다는 것을 확인 했다면 간단하게 단백질이 접히는데 소요되는 시간을 우리는 근사할 수 있을 것이다. 물론, 4개의 전산모사 과정을 우리가 알고 있기 때문에 어떻게 단백질이 접혀지는지도 조사할 수 있을 것이다. 단백질이 접히는데 소요되는 시간에 대한 근사값은 (1000x10 ns)/4 로 할 수 있다. 이는 단백질 구조가 접힐 확률이 지수함수 꼴을 갖는다는 가정을 하고 접힐 확률이 매우 낮다는 가정을 사용하면 쉽게 얻을 수 있는 공식이다. 이러한 연구방법은 실제 다양한 소규모 단백질에 적용되었다. 연구결과 실제 단백질의 접힘 과정을 잘 설명하고 있다.

  이러한 연구 방법은 굉장히 일반적인 연구 방법일 수 있으나, 해당 화학반응의 활성화 에너지가 상대적으로 높은 경우에는 여전히 적용하기가 어렵다. 왜냐하면 보다 많은 컴퓨터들을 동원한다고 하더라도 관측하고자 하는 분자세계의 희귀사건이 일어날 확률은 매우 낮기 때문이다. 이러한 계산은 주로 스크린 세이버 형식으로 실행되며 대량의 PC 컴퓨터들을 동원할 수 있을 때 유용하다. 이러한 계산 방식은 분산 컴퓨팅(distributed computing) 기술의 일종으로 PC 간의 정보 교환은 무시할 정도 이거나 아예 없는 경우의 응용계산에 적합하다.

   또 달리 희귀사건을 연구하는 방법으로는 포텐셜 에너지 표면을 인위적으로 재구성해서 통상적인 분자동역학 전산모사를 행하는 방법이 있다. Voter 박사 그룹에서 개발된 hyperdynamics 방법은 위에서 언급한 포텐셜 에너지 표면을 바꾸어서 분자동역학 계산을 수행하는 방법의 전형적인 예이다. 포텐셜 에너지가 낮은 구간을 확인하고 그 구간에서는 의도적으로 포텐셜 에너지 표면을 변형시켜서 또 다른 에너지 구간으로 원자들이 쉽게 전이할 수 있도록 만들어 주는 것이다. 이 때, 상대적으로 높은 에너지 구간에 해당하는 곳들에서는 전이경로가 당초에 정의된 포텐셜 에너지 함수로부터 유도되도록 한다. 즉, 변형된 포텐셜 에너지 함수와 변형되지 않은 포텐셜 에너지 함수를 이용하여 전이경로를 만들어 내는 계산방법이다. 전이경로에 해당하는 부분에서는 변형되지 않은, 조작되지 않은 포텐셜 에너지 함수를 이용하여 분자동역학 계산을 수행한다. 이렇게 분자동역학 방법을 수행할 경우, 물리적으로 일어날 가능성이 있는 희귀사건의 출현이 빈번해 질 수 있다. 분자동역학이 가속화되었다고 생각할 수 있다. 실제 분자동역학 계산을 수행한 시간과 물리적으로 해석 가능한 희귀사건들이 발생하면서 지나간 시간 사이에는 엄청난 차이를 보일 수 있다. 예를 들어 10 ps를 가속화된 분자동역학 방법으로 전산모사 했지만, 실제로 물리적으로는 10 ns 이상의 물리적 시간이 흘러간 것으로 해석할 수 있다. 전이경로를 연구하는 입장에서 전이경로가 만들어지지 않을 상황들, 즉, 에너지 웅덩이 근처에서 물리적 시간을 소모하는 상황들을 의도적으로 피해보자는 것이다. 이것이 에너지 웅덩이 근처에서 의도적으로 추가된 포텐셜 에너지 함수가 하는 일이다.

  포텐셜 에너지 함수값이 매우 낮은 곳들에서는 인위적으로 분자구조가 오래 머물러 있지 못하게 하는 포텐셜 에너지 함수 항을 도입한다. 따라서 전산모사 시간에 비해서 실질적으로는 오랜 시간을 전산모사 한 것과 같은 효과를 볼 수 있다. 의도적으로 온도를 높여서 희귀사건이 잘 일어나도록 해서 분자동역학 방법을 수행할 수 있다. 이러한 아이디어 자체는 좋은 아이디어이다. 다만, 높은 온도에서 일어나는 희귀사건과 낮은 온도에서 일어나는 희귀사건의 전이경로들이 다를 경우 실제 전산모사에서는 높은 온도와 관련된 전이경로만 관측될 가능성이 있다. 이러한 인위적인 온도의 도입과 온도 차이를 보정하여 가속화된 분자동역학을 수행하는 계산방법도 발표 되어 있다. Parallel replica 방법에서는 주어진 온도에서 가능한 서로 다른 속도 분포들을 사용하여 완전히 독립적인 분자동역학 계산을 여러 대의 컴퓨터에서 수행한다. Hyperdynamics 방법과 parallel replica 방법을 혼용하여 전산모사를 수행할 수도 있을 것이다. 전이경로가 만들어지는 초기부분은 비정규적인 다소 불분명한 상태에서 시작했지만 전이경로의 중요한 부분, 핵심 부분은 미리 가정한 포텐셜 에너지 함수에 의해서 생성된다. 결국, 가속화된 분자동역학 방법은 전이경로 계산의 한 형태로 볼 수 있다. 다차원의 공간에서 인위적으로 희귀사건을 유도하게 만들어주는 포텐셜 에너지 함수를 설계하는 것이 단순한 형태의 일이 아니기 때문에 hyperdynamics 방법의 적용 범위가 광범위하게 넓지는 못하다. 향후 좋은 아이디어가 기대되는 연구 방법이라고 볼 수 있다. 

  희귀사건을 연구할 수 있는 일반적인 방법으로 Chandler 교수 그룹에서 개발된 transition-pathway sampling 방법이 있다. 이 방법에서는 반응전 그리고 반응후의 분자상태(원자배열)들을 미리 정의하고 이들 사이의 동역학적인 원자 이동경로들을 분자동역학 방법으로 직접 찾는 방법이다. 이 계산방법은 계산의 편의상 선택된 하나의 편향된 전이경로로부터 시작하여 분자동역학 방정식 조건을 만족하는 전이경로를 찾는 방법을 포함한다. 반복적인 계산을 통해서 하나의 전이경로로부터 에너지적으로 매우 유리한 전이경로를 찾는 방법이다. 매우 일반적인 계산방법이나, 매우 많은 분자동역학 계산수행이 필수적이다. 여전히 분자동역학 방법을 활용하지만 화학반응이 일어나기전의 ‘기다리는 시간’을 빼고 계산할 수 있는 장점이 이 방법에는 있다. 각종 화학반응들, 결정 응결, 단백질 접힘, 그리고 상전이 현상 등에 성공적으로 응용되었다. 다양한 전이경로들을 조사할 수 있기 때문에 반응기구, 활성화 에너지, 반응율 등을 계산할 수 있다.

  NEB(nudged elastic band) 방법(Jo'nsson 교수 그룹), string 방법(E 교수 그룹)은 두 개의 분자 구조들이 주어져 있을 때 그 들 중간에 위치한 전이상태를 찾는 방법이다. Dimer 방법(Jo'nsson 교수 그룹)은 주어진 하나의 분자 구조로부터 근처에 위치한 전이상태를 찾는 방법이다. 실제 계산에서 가상의 분자 구조 두 개를 사용하기 때문에 붙여진 계산방법 이름이다. Dimer 방법과 마찬가지로 ART 방법(Mousseau 교수 그룹)도 주어진 하나의 분자구조 근처에 위치한 전이상태들 중 하나를 찾는 방법이다. 포텐셜 에너지의 이차 미분을 필요로 하지 않은 계산 방법들이 원자수를 많이 포함한 분자계에 적용하기 쉽다는 면에서 유리하다. 작용유도 분자동역학(action-derived molecular dynamics; ADMD) 방법은 최소작용 원리를 이용한 계산 방법이다. ADMD(Parrinello 교수 그룹, Elber 교수 그룹), NEB, string 방법들은 동시에 여러 가지 모양의 분자 구조들(반응이 일어나는 순서에 따른 분자구조들)을 한꺼번에 다루는 계산방식을 취한다. 분자구조들이 하나의 분자 이동경로를 나타내는 양식으로 계산이 수행되어 자연스럽게 병렬화가 가능하다. 이것은 작용유도 분자동역학 방법이 시간적으로 순차적인 분자동역학 방법과 구별되는 특징이다.

  Dimer나 ART방법으로 유도된 전이상태를 이용하여 분자구조들을 변형시키고 Monte Carlo 방법을 활용하여 시간의 흐름을 도입하면 일반적인 kinetic Monte Carlo 방법이 될 수 있다. 다양하고도 많은 전이상태들을 넘어서 최종적으로 안정된 하나의 상태에 이르는 전산모사를 수행할 수 있다. 이 때 시간의 흐름은 관련된 전이상태를 통과하는데 걸리는 시간을 확률적으로 환산하여 책정한다. 하나의 전이상태를 지나는데 소요되는 시간은 조화 전이상태 (harmonic transition state) 이론을 활용하여 계산한다. 어떠한 전이경로를 지나면서 원자배열 상태가 바뀌느냐에 따라서 해당 원자배열의 변화에 걸리는 물리적인 시간은 다르게 책정된다. 이러한 kinetic Monte Carlo 방법은 미리 정해둔 사건들만을 활용하는 kinetic Monte Carlo 방법과는 구별된다는 점을 지적해야한다. 미리 정해둔 가능한 분자모양 변화들 또는 사건들의 리스트를 사용하지 않는 양식의 전산모사라는 점이 Dimer 방법을 활용한 kinetic Monte Carlo 방법의 특징이다. 일어날 사건은 프로그램 실행 중에 계산한다는 뜻이다. 통과할 전이상태들, 각각의 전이상태들을 통과하는 데 소요되는 시간 등을 전산모사를 진행 하면서 계산한다는 것이다. Dimer 방법을 활용한 kinetic Monte Carlo 방법은 다양한 형식의 분자모양 변화를 허용하는 보다 일반적인 전산모사 방식으로 평가받을 수 있다.

  금속 표면에 추가로 원자들이 장시간에 걸쳐서 증착될 때, 금속 표면에 형성되는 안정된 원자들의 모임을 실험적으로 관찰할 수 있다. 이와 같은 장시간에 걸쳐서 진행되는 실험적 상황은 기존의 이론적 접근법으로는 전산모사가 거의 불가능하다. Jo'nsson 교수 그룹에서는 Dimer 방법과 kinetic Monte Carlo 방법을 활용하여 장시간(1 msec)에 걸쳐서 원자들이 재배열되는 상황을 전산모사 했다. 원자 배열이 어떠한 양식들로, 어떠한 상태들을 거쳐서 최종적으로 안정화된 모양에 이르는지를 전산모사를 통해서 알아내었다. 또한 소요되는 시간도 계산할 수 있음을 보였다. 격자모형 없이 그리고 미리 정해둔 사건목록 없이 kinetic Monte Carlo 계산이 가능함을 보여주었다.


수직으로 찔려진 고무 밴드(NEB)

  수직으로 찔려진 고무 밴드 이것은 전이상태를 계산하는 하나의 방법 이름이다. 동쪽과 서쪽을 가르는 산맥이 있다고 가정할 때, 동쪽 평원에 있는 마을에서 서쪽 평원에 있는 또 다른 마을 까지를 하나의 고무 밴드로 연결한다고 생각하면 이 계산방법을 이해하는데 편하다. 가장 낮은 높이의 골짜기를 찾아서 고무 밴드를 연결할 수 있을 것이다. 처음에는 고무 밴드를 동서 방향으로 늘여놓고, 고무 밴드를 고무 밴드가 놓인 동서 방향에 수직한 남북방향으로 움직여서(찔러서) 가장 높이가 낮은 골짜기를 고무 밴드가 지나가게 만들 수 있을 것이다. 고무 밴드위의 여러 점들을 각각 고무 밴드의 수직 방향으로 움직일 수 설정할 수 있다. 이것이 nudged elastic band(NEB) 방법의 핵심 아이디어이다.

  양쪽 끝 단 두 지점에 고무 밴드의 고정시키고 고무 밴드 위의 여러 점을 생각한다. 고무 밴드 위에서 선택된 점들은 복제된 원자배열(또는 분자구조)을 나타낸다. 고무 밴드 위에서 첫 번째와 마지막에 선택된 점들은 각각 반응전, 반응후의 원자배열을 나타낸다. 그 밖의 많은 점들은 각각 하나의 원자배열을 표시하고 NEB 계산을 통해서 다시 결정되어야 할 원자배열들을 나타낸다고 볼 수 있다. 고무 밴드 중간에 위치한 원자배열을 복제품이라고 하여 replica 또는 image라는 용어를 사용하는 연구자들이 있다.

   고무 밴드 위의 각각의 점들은 고무 밴드의 수직 방향으로 찔려져서(nudged) 변형될 수 있다. 수직 방향으로만 변형시키지 않을 경우 고무 밴드 위의 점들 사이의 간격이 심하게 불규칙적으로 바뀔 수 있다. 고무 밴드의 수직한 방향은 고무 밴드 위에 있는 점의 위치에 따라 다르며 인접한 고무 밴드 위의 점들을 활용해서 정의할 수 있다. 이 때, 고무 밴드 위의 한 점이 찔러서 그 위치를 바꾸어지게 해주는 힘은 두 가지 형태가 있다. 하나는 고무 밴드 자체의 탄성으로부터 오는 것이다. 또 다른 하나는 고무 밴드 위의 한 점, 즉, 원자배열 자체에서 원자들이 받는 힘들로부터 나온다. 이 두 가지 힘을 이용하면 고무 밴드 위의 각 원자들은 힘을 받지 않은 다른 곳으로 이완 시킬 수 있다. 일반적으로 여러 번의 이완 과정을 통해서 고무 밴드의 펼쳐진 모양이 바뀌게 된다. 처음 고무 밴드가 일직선처럼 놓여있었다고 하더라도 고무 밴드 위의 각점들이 놓인 상황(포텐셜 에너지 함수값)이 같지 않기 때문에 고무 밴드는 굽어지고 전체 모양이 굽이쳐 펼쳐질 수 있을 것이다. 이러한 고무 밴드의 휘어짐은 주어진 포텐셜 에너지가 낮은 쪽으로 휘어질 수밖에 없다. 물론, 고무 밴드의 탄성 때문에 포텐셜 에너지가 낮은 쪽으로 무작정 휘어질 수는 없다. 전체적으로 고무 밴드가 모양을 바뀌어 가다가 더 이상 고무 밴드의 전체의 모양이 바뀌지 않을 때, 최고로 높은 골짜기를 근처의 한 점(원자배열)을 전이상태로 생각할 수 있다. 물론, 얼마나 많은 image들을 사용하는가에 따라서, 전이상태 결정의 정밀도가 정해질 것이다. 전이상태에 해당하는 에너지는 참값 보다 낮게 또는 높게 계산될 수 있다. 전이상태가 대체로 간단한 경우에는 이 방법을 이용하면 쉽게 전이상태를 결정할 수 있다. 상당히 큰 분자구조에서 일부분에서만 원자들이 활발히 움직이는 것이 확실 할 때에는 반응에서 활발하지 않은 원자들을 분리해서 방관하는 원자들로 만들 수 있다. 이렇게 할 경우 실제 고무 밴드 방법에서 실질적으로 결정해야하는 원자들의 움직임과 관련된 자유도가 줄어들기 때문에 보다 빨리 전이상태를 결정할 수 있다.  

  Jo'nsson 교수 그룹에서 개발된 NEB 방법은 그 알고리듬의 편리한 적용성 때문에 많은 연구자들에 의해서 성공적으로 응용되고 있다. NEB 방법의 기본 아이디어를 활용하여 각종 응용문제의 특성에 맞추어 개발된 유사한 방법들도 많이 존재한다. NEB 방법은 결국 "chain-of-states"를 직접 취급하는 방식이다. 고무 밴드 위의 점들은 실제로는 ‘스프링‘으로 연결되어 있다. 연속한 분자구조들 사이의 실질적인 거리는 ‘스프링’에 의해서 결정된다. 적절한 ‘스프링’의 강도 선택을 위해서는 몇 가지 테스트 계산이 필요하다. NEB 방법의 특징은 목적함수 없이 원자배열들에 대한 이완작업을 수행한다는 점이다. NEB 방법은 두 분자구조들 사이의 가장 에너지적으로 쉬운 연결 방식을 찾고 궁극적으로는 전이상태를 찾는 것을 목표로 하는 알고리듬이다. NEB 방법은 각 원자배열들에 대해서 계산적으로 독립된 에너지, 일차미분(force) 계산들을 필요로 한다. 따라서 NEB 방법은 쉽게 병렬화가 가능하다.  NEB 방법은 병렬 효율성이 매우 뛰어나다. "스프링"은 분자구조들을 conformational space 안에서 일정하게 배열시키려고 마련해 둔 장치이다. 특별히 "스프링" 에 의해서 연결된 분자구조들 사이에서 "스프링" 방향에 수직한 방향으로의 구조적 이완(relaxation)을 허용한다. 분자구조들 사이의 일정한 거리를 위해서 수직한 힘 성분들에 의한 이완을 이용한다. 사실, 이 이완들을 통해서 보다 에너지적으로 유리한 원자배열들을 확보할 수 있다.

  string 방법도 NEB방법과 유사한 방법이다. string 방법에서는 원자배열들에 대한 직접적인 이완을 활용하는 대신에 string(유연하게 움직이는 수학적 곡선을 표현하는 것으로 NEB 방법의 elastic band라고 볼 수 있다.)의 이완을 중요시 한다는 점이 NEB 방법과의 차이점이라고 할 수 있다. 추가적으로 string 길이에 대한 내부적 조건을 알고리듬 스스로 결정, 활용한다는 점이 특징이다. string 방법도 NEB 방법과 마찬가지로 목적함수 없이 원자배열들에 대한 이완작업을 수행한다.

  유한온도 string 방법은 string 방법을 유한온도 영역으로 확장한 계산방법이다. string의 방향벡터 수직한 hyperplane 상에서 분자동역학이나 Monte Carlo 샘플링을 통해서 주어진 온도 조건에서 가능한 원자배열들을 총체적으로 계산하여 하나의 string이 아닌 하나의 에너지적으로 허용된 ‘전이경로튜브’를 계산한다. 이러한 ‘전이경로튜브’는 전이경로들의 모임이 유한온도에서 고유의 폭을 가지고 분포하고 있는 것을 의미한다. 온도가 0 K로 수렴할 때, ‘전이경로튜브’의 직경은 0으로 수렴한다. 이 때의 전이경로가 바로 string 방법에서 찾은 전이경로(최소-에너지-경로)가 됨을 의미한다. 이러한 계산으로부터 자유에너지 프로파일을 계산할 수 있다.  

  Elber 박사가 제안한 전이경로 계산방법 중에는 고무 밴드 방법과 유사한 개념으로 전이경로를 계산하는 방법이 있다. 각각의 replica에서 고무 밴드를 따라가면서 분산된 양식의 고무 밴드의 방향벡터를 정의한다. 이 방향벡터에 수직한 다차원 공간상에서의 평면(hyperplane)을 생각할 수 있고 이 평면 위에서 원자배열을 에너지적으로 최적화 할 수 있다. 이 방법은 직관적이고 쉽게 응용할 수 있은 장점을 가지고 있어서 현재 폭넓은 분야에서 응용되고 있다. 조화 전이상태 이론(harmonic transition state theory)에 의하면 안정한 구조, 전이상태 구조에서 가능한 분자 진동 스펙트럼을 조사하면 주어진 온도 하에서 활성화 에너지와의 조합으로 해당 전이상태를 통과할 확률을 계산할 수 있다.


작용유도 분자동역학(ADMD)

   화학반응은 원자간 결합의 파괴와 재생성 등이 연관된 현상으로 매우 복잡하고 다양하다. 이러한 화학반응 경로에 대해서 지속적으로 응용된 하나의 계산방법은 존재하지 않았다. 체계적이고 근원적 이론에 근간을 둔 계산방법이 필요한 시점이다. 이러한 관점에서 최근에 활발히 응용되기 시작한 작용유도 분자동역학(action-derived molecular dynamics) 방법은 희귀사건 기술에 많은 가능성을 보여주는 계산방법으로 기대를 모으고 있다. 작용유도 분자동역학 방법은 최소 활성화 에너지를 가지는 전이경로 모델을 구축하는데 매우 유용한 방법이다. 최소작용 원리(least action principle)를 컴퓨터 상에서 구현한 것이다. 고전역학에서 논의되는 물리량 작용(action)은 원자들의 이동경로의 함수이다. 작용 함수가 최적화되는 조건에서 원자들의 이동경로를 결정할 수 있다. 이것이 최소 작용원리 (해밀튼 원리)이다. 이렇게 구한 원자 이동경로는 뉴튼의 운동 방정식을 만족한다. 컴퓨터 상에서는 원자들의 위치를 시간에 대해서 연속적으로 표현할 수 없다. 원자들의 속도, 운동 에너지는 유한차이(finite-difference) 법을 활용하는 경우가 많다. 이러한 상황에서는 작용함수를 해석적인 함수처럼 취급하기 곤란 점들이 있다. 시간 간격의 크기 선택에 따라서 작용함수의 전체적인 모양이 심하게 바뀔 수 있기 때문이다. 예를 들면, 시간 간격의 크기의 선택에 따라서, 특정한 원자 이동경로에서 최소점을 가지던 작용함수가 최대점을 가지는 작용함수로 바뀔 수 있다. 일반적인 응용연구에서 작용함수의 형태를 알지 못한 상태에서 문제를 풀기 때문에 작용함수의 최소화 또는 최대화를 결정하기가 힘들다.

  작용유도 분자동역학에 사용되는 운동에너지를 가장 낮은 차수의 유한차이 방법으로 표현했을 때에는 해석적으로는 Verlet 운동 방정식을 얻을 수 있으나, 실제 수치적 응용에서는 총에너지 보존이 이루어지지 않음을 확인할 수 있다. 이러한 문제점을 해결하기 위해서 Passerone, Parrinello 박사들은 총에너지 보존이 되도록 벌칙항(penalty term)을 추가로 도입하였다. 결국, 목적함수는 작용함수와 벌칙항의 합으로 구성되어 있다. 이렇게 할 경우, 작용함수의 최적화문제는 벌칙항을 포함하는 확장된 작용함수(목적함수)의 최소화 문제로 귀결되고 Verlet 적분경로와 유사한 근사된 원자들의 운동경로를 계산할 수 있다. 물론, 이렇게 구한 원자이동 경로들로부터는 총에너지가 보존됨을 확인 할 수 있다.

   Passerone, Parrinello 박사들이 제안한 분자동역학 방법에 의한 최종적인 전이경로는 계산을 시작할 때 사용하는 초기경로에 상당히 의존하는 특징이 있다. 이러한 약점을 보완하기 위한 하나의 확장된 목적함수가 발표되었다. 이 목적함수는 Passerone, Parrinello 박사들이 제안한 기존의 작용함수에 추가로 운동에너지를 조절하는 항이 있다. 최종적으로 계산된 전이경로의 초기경로 의존성을 상당히 완화 시킬 수 있음을 확인했다. 더욱이 각 원자에 대한 운동에너지 조절항을 목적함수에 추가적으로 도입하면 최종적으로 계산된 원자들의 운동경로가 Verlet 운동경로 모형에 더욱 근접할 수 있음을 보였다. 이러한 벌칙항들은 원자 이동경로 계산에서 나타날 수 있는 몇 가지 수치적 어려움들을 해결한 것으로 해석할 수 있다.

  통상의 분자동역학 방법에서처럼 포텐셜 에너지 함수와 각 원자에 작용하는 힘 계산만으로 원자의 이동경로를 계산할 수 있는 작용유도 분자동역학 방법은 폭넓은 분야에서 응용이 가능하다. 거대 시스템에 적용 가능한 계산방법이 되기 위해서는 그 계산방법이 단순한 복잡도를 가져야한다. 이점에서 앞서 Elber 박사 등이 제안한 작용유도 분자동역학 방법은 다소 불리한 면이 있다. 이 방법에서는 포텐셜 에너지 함수의 원자들 위치에 대한 이차미분을 필요로 하여 수치적 계산 시 비교적 컴퓨터 자원을 많이 사용하게 되는 불리한 측면이 있다. 하지만 이론적으로 목적함수 자체가 잘 정의되고 목적함수의 최소값이 존재한다는 특징들이 있다.

  작용유도 분자동역학은 화학반응을 기술함에 있어서 반응좌표를 미리 설정할 필요가 없는 방법으로 상당히 일반적인 계산방법이다. 화학반응이 복잡할수록, 해당 반응을 단순하게 쉽게 기술하는 반응좌표를 선택하기가 매우 어렵다. 일반적인 선택기준이 알려져 있지 않은 반응좌표의 설정 없이 화학반응을 연구할 수 있다는 것은 작용유도 분자동역학 방법이 매우 선진화된 계산방법임을 의미한다. 이것은 작용유도 분자동역학 방법의 중요한 특징으로 향후 많은 응용연구가 수행될 수 있음을 의미한다.

  통상 분자의 크기가 큰 경우에 전산모사 하기가 힘들어진다. 이는 공간적인 크기에 의해서 분자동역학 연구가 제한된다는 뜻이다. 시간 순차적으로만 계산이 진행되어야만 하는 특성 때문에 분자동역학의 병렬화 또한 단순하지 않고 그 효율성에 있었어도 한계가 있다. 실제 분자동역학 응용에서는 "크기"에 관한 제한뿐만 아니라, "시간"에 관한 제한 또한 만만하지 않다. 다룰 수 있는, 즉, 전산모사 할 수 있는 "시간"에 제한이 있다. [그림3]에서는 작용유도 분자동역학 방법이 응용될 수 있는 적용범위에 대해 도식적으로 나타내었다. 작용유도 분자동역학에서는 통상의 분자동역학 방법에서와는 달리 상당히 큰 시간 간격을 사용할 수도 있다. 이렇게 할 경우, 분자모양의 변화를 기술함에 있어서 사용된 시간 간격보다도 짧은 시간 안에 나타나는 원자들의 움직임은 무시된다.

그림 3)_ 작용유도 분자동역학 방법의 적용범위를 다른 계산방법들과 비교했다. 화학반응(원자 재배열)과 같이 원자의 진동 주기에 비해서 매우 긴 시간 스케일에서 일어나는 원자들의 집단 운동을 기술하기에 편리한 방법이 작용유도 분자동역학(ADMD) 방법이다. 제일원리 분자동역학(ab initio Molecular Dynamics)방법과 고전 분자동역학(classical Molecular Dynamics) 방법에서는 원자들의 이동경로 적분을 위해서 매우 짧은 시간 간격(1 fs)을 사용해야한다. ADMD 방법에서는 상대적으로 긴 시간 간격을 활용할 수 있다. 운동 몬테칼로 (kinetic Monte Carlo) 방법의 경우, 새로운 원자배열이 이루어질 때 지나는 전이상태의 특성과 전이상태 전의 원자배열 상태에 따라서 원자배열에 소요된 시간이 책정된다. 운동 몬테칼로에서는 시간의 흐름이 불규칙적으로 일어남을 의미한다. 원자배열(N 개의 원자들)이 주어졌을 때, 포텐셜 에너지 함수를 계산하는 과정의 복잡도(complexity)는 제일원리 분자동역학의 경우 O(N3)에 해당한다. 고전 분자동역학의 경우 복잡도는 O(N log N)인 경우가 많다. 운동 몬테칼로와 작용유도 분자동역학 방법들의 경우 사용하는 포텐셜 에너지 함수의 종류에 따라서 달라 질수 있다.


  작용유도 분자동역학 방법은 다양한 분야에 적용될 수 있는 일반적인 계산방법이다. 다양한 포텐셜 에너지 함수들에 적용될 수 있다. 표면 위에서의 원자들의 집단 운동, 단백질 접힘, 분자들과 원자들의 결합 분리, 촉매역할을 포함한 일반적인 화학반응 기작 연구에 사용될 수 있다. 또한 에너지적으로 유리한 반응경로(최소-에너지-경로)를 적극적으로 찾아낼 수 있다. 총에너지 보존과 관련된 벌칙항에 운동학적으로 허용된 총에너지 값을 사용자가 직접 넣어 줄 수 있다. 즉, 관심이 있는 에너지 영역에서 하나의 낮은 활성화 에너지를 가지는 전이경로를 적극적으로 찾아낼 수도 있다. 아레니우스 공식에서 알 수 있듯이 실현 가능성이 높은 전이경로는 해당 활성화 에너지로 쉽게 평가할 수 있다.

  
작용유도 분자동역학 방법은 태생적으로 병렬화에 매우 유리한 형식의 계산방법이다. 작용유도 분자동역학 방법에서는 하나의 원자 이동경로가 다른 형식의 원자 이동경로로 변환하는 단계가 하나의 계산단위로 요구되는 특성을 가지고 있다. 작용항, 벌칙항들의 계산 그리고 관련된 목적함수 기울기 벡터의 계산 부분이 가장 많이 컴퓨터 자원을 소모되는 부분이다. 원자의 이동경로가 순차적으로 표현되는 분자동역학 계산과 좋은 대조를 이룬다. [그림4]와 [그림5]에서 확인할 수 있는 것처럼 작용유도 분자동역학 방법은 쉽게 병렬화 되고 병렬 효율성도 매우 높다.

  작용유도 분자동역학에서의 문제풀이는 수학적으로 경계조건 문제 (boundary value problem) 형식이다. 반면, 분자동역학 방법은 초기조건 문제(initial value problem) 형식이다. 질량 중심의 이동과 회전운동에 대한 자유도는 쉽게 제거할 수 있다[S. J. Kearsley, J. Comput. Chem. 11, 1187 (1990)]. 연속해서 만들어지는 분자구조들 사이의 병진, 회전운동에 따른 잉여분의 자유도를 제거하는 조건인 Eckart 조건을 활용한 계산방법도 충분히 고려해 볼 수 있다. Eckart 조건을 만족하는 image들을 만드는 것은 NEB, string, ADMD 방법들에서 공히 활용 가능하다.

  작용유도 분자동역학 방법으로 계산한 원자 이동경로들을 Verlet 운동 방정식을 어느 정도 충실히 만족하는지를 살펴볼 필요가 있다. 사실, 유한차이 공식을 이용하여 분산된 작용을 정의했다. 이러한 분산된 작용을 사용하여 해석적으로 운동방정식을 유도하면 Verlet 경로를 얻을 수 있다. 하지만, 작용유도 분자동역학 방법을 활용하여 얻어진 원자 이동경로는 Verlet 운동경로와 다소 차이를 보인다. 대신에 계산된 전이경로는 총에너지보존 법칙을 만족한다. 총에너지보존이 이루어지지 않고 Verlet 경로를 최대한 따르도록 최적화된 전이경로를 찾는 방법이 Elber 박사가 제안한 계산방법이다.

  한 예제 문제에 대해서 작용유도 분자동역학을 적용한 후 원자 이동경로가 Verlet 경로로부터 벗어난 정도(error variable)의 분포를 [그림6]에 나타내었다. 작용유도 분자동역학 응용문제들 대부분의 경우 원자들의 정확한 이동경로를 모른다. 따라서 여기서 논의하는 error variable은 정확한 정답에 해당하는 전이경로와 비교했을 때 오류의 정도를 나타내는 양과는 차이가 있음을 지적해야한다. 사실, 정확히 뉴튼의 운동방정식을 만족해도 분산된 Verlet 운동 경로로 표시할 경우 가우시안 형태로 error variable들이 분포함을 보일 수 있다.

그림 4)_ 작용유도 분자동역학 방법의 순서도(flowchart)를 표시했다. 일반으로 주어진 분자구조에 대한 포텐셜 에너지 함수값 계산과 각 원자에 작용하는 힘을 계산하는 부분이 가장 많은 컴퓨터 자원(CPU 시간, 메모리)을 필요로 하는 계산단계이다. 다양한 분자구조에 대해서 포텐셜 에너지 함수값과 각 원자에 작용하는 힘 계산 부분을 병렬처리 한 것이 특징이다. MPI 라이버러리를 이용하여 쉽게 병렬화 할 수 있다. 경로 최적화 과정에서 새로운 시도경로는 LBFGS 루틴의 출력을 통해서 얻을 수 있다. LBFGS 루틴은 원자 이동경로를 입력으로 하여 새로운 원자 이동경로를 출력으로 하는 형식으로 작동한다. 연속적으로 LBFGS 루틴을 호출함으로써 수렴된 원자 이동경로를 얻을 수 있다. 불연속적인 원자들의 위치를 유한한 개수의 시간단계(j=0,1,2,...,P)를 따라서 설정하고 각 단계에서 원자들의 집단운동(concerted motion)을 기술하려는 것이 작용유도 분자동역학이다. 초기, 말기 상태의 분자구조(원자배열)는 각각 시간단계 지수 j=0, j=P에 해당한다. 실질적으로 변화하는 변수들의 수는 (P-1)개의 분자구조들로부터 나온다. 예를 들어 N개의 원자들로 만들어진 분자를 다루는 경우 작용유도 분자동역학 방법으로 결정되어야 할 변수들의 수는 3N(P-1)이다. 계산종료 직전 단계에서 시간단계 지수를 따라서 포텐셜 에너지 함수, 분자구조의 변화, 총에너지 보존의 정도를 출력한다.


error variable들의 분포를 보면 0 에 가까운 평균값을 가지고 있다. 이는 Verlet 운동경로와 상당히 유사한 경로가 작용유도 분자동역학 방법으로 계산되었다는 것을 의미한다. 편차는 총에너지 보존과 관련된 벌칙항 이외에 각 원자들의 운동 에너지 조절과 관련된 벌칙항을 도입함으로 더욱 작아질 수 있음이 알려져 있다. 통상의 분자동역학에서는 다룰 수 없는 원자세계의 희귀사건들을 시간에 대한 경계조건으로 변형시켜서 원자 군의 집단적인 운동(concerted motion)에 기인한 운동 특성을 연구할 수 있도록 설계된 것이 작용유도 분자동역학 계산 방법이다. 계산된 전이경로로부터 원자들이 얼마나 뉴튼의 운동방정식을 충실히 따르는가를 평가할 수 있기 때문에 계산된 전이경로 정확도에 관한 논의가 정량적으로 가능하다. 작용유도 분자동역학 방법에서 선택되는 총시간은 계산에서 사용되는 하나의 입력값이다. 이 값은 특정한 값보다 작을 경우에는 최종적으로 계산된 전이경로는 활성화 에너지가 높게 나올 가능성이 많다. 이는 충분한 시간을 가지지 않고 급격하게 전이경로가 만들어진 경우에 해당한다. 충분히 긴 시간을 입력으로 사용할 경우에는 이러한 문제는 발생하지 않는다. 기준 총에너지값을 활용한 최소-에너지-경로를 찾는 과정에서 전이경로 확립에 필요한 물리적 시간과 반응 전후의 ‘기다리는 시간’의 분리가 자연스럽게 얻어진다.

   
[그림7]에서는 작용유도 분자동역학 방법과 NEB 방법을 비교했다. 탄소60 분자 두 개가 융합하여 하나의 캡슐 모양의 분자 탄소120을 만드는 매우 복잡한 화학반응 문제에 대해서 작용유도 분자동역학 방법과 NEB방법을 각각 적용했다. Tersoff 박사가 제안한 간단한 경험적 포텐셜 에너지 함수를 이용하였다. 이 문제를 통하여 매우 복잡한 화학반응 문제풀이의 경우 NEB 방법상으로는 상대적으로 낮은 활성화 에너지를 가지는 원자 이동경로 또는 전이상태를 찾는 것이 쉽지 않음을 확인했다. 반면, 작용유도 분자동역학 방법에서는 원하는 총에너지 값 근처에서 적극적으로 전이경로를 찾는 것이 가능하다. 물론, 작용유도 분자동역학 방법이 훨씬 많은 포텐셜 에너지 함수의 계산들을 요구한다. 하지만, 매우 복잡한 문제풀이의 경우 NEB 방법을 활용한 문제풀이가 작용유도 분자동역학 방법을 활용할 때보다도 더 강력하지 않음을 알 수 있다. 이는 목적함수를 활용하고 있는 작용유도 분자동역학 방법의 장점이다. 비교적 단순한 문제풀이의 경우 NEB 방법을 활용하면 효율적으로 전이상태를 계산할 수 있다.


그림 5)_ 탄소60 분자 두개가 융합하여 탄소120 갭슐 형태의 분자로 변환하는 과정을 작용유도 분자동역학 방법으로 기술했다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용했다. 64개의 프로세서를 사용할 경우 247초 후에 모든 계산이 종료되었다. 벽걸이 시간(wall-clock time)을 이용하여 병렬 효율성을 표시한 것이다. IBM p690 컴퓨터를 활용하여 테스트를 진행했다. IBM p690에 설치된 hpmcount(V2.4.3) 프로그램을 활용하여 최적화 정도/병렬 효율성을 판단했다. 예를 들면 프로그램 실행 후 64개의 hpmcount관련 화일들이 생긴다. 각 화일에는 wall clock time이 출력되어 있다. 64개중 가장 큰 값을 이용하여 소위 병렬 계산 speedup을 계산했다. 하나의 프로세서를 이용한 경우 wall-clock time 과 비교하여 얼마나 더 빨리 작업이 종료되었는지를 계산했다. 출력은 프로세서의 숫자에 관계없이 정확히 같은 물리적 해를 주고 있으나 실질적인 계산수행에 소요되는 시간상에서 효율이 나타나는 것이 병렬계산의 특징이다. 단일 프로세서 계산에서 1.11 GFLOPS의 계산능력을 보였다. IBM p690 슈퍼컴퓨터가 가진 이론적인 계산처리능력은 6.8 GFLOPS이다.

그림 6)_ 하나의 작용유도 분자동역학 계산에서 얻어진 원자 이동경로들에 ‘정밀도’를 통계적으로 분석 표시했다. 작용유도 분자동역학으로 얻어진 원자들의 이동경로에 Verlet 공식을 그대로 적용하여 Verlet 공식에서 벗어난 정도(error variable)를 수평축의 양으로 두고 그 분포를 수직축에 표시했다. 분포함수가 delta 함수일 경우, 원자들의 이동경로는 정확히 Verlet 경로를 따른다는 뜻이다. 수평축의 단위는 길이 단위이다. 동원된 데이터의 수는 작용유도 분자동역학 문제 풀이의 미지수의 개수와 같이 3N(P-1)이다. 여기서, N, P는 각각 원자들의 수와 시간단계의 수를 나타낸다. 실제 응용계산에서는 N=122, P=300, 시간간격으로 6 fs를 사용하였다. 전이경로의 시작부분과 끝나는 부분에서는 error variable 들의 절대값이 상대적으로 큰 값들을 가진다. 평균값, 0 근처의 값을 가지는 경우는 대부분 전이경로의 중간 부분에 해당하는 경우가 많다. 운동 에너지 조절항을 포함한 확장된 목적함수를 이용할 경우 원자 이동경로의 정확도가 Onsager-Machlup 작용으로 평가했을 때 상당히 개선됨을 보였다. 운동 에너지 조절항을 포함한 목적함수는 원자 이동경로의 '정밀도'을 높이는데 사용될 수 있을 뿐만 아니라 원자 이동경로의 담금질에도 사용될 수 있다. a.u.=0.529177Å


그림 7_ 탄소60 분자들 두 개의 융합과정에서의 포텐셜 에너지의 변화를 작용유도 분자동역학 방법과 nudged elastic band 방법으로 각각 계산했다. 각 분자 구조에 해당하는 에너지 계산을 위해서 Tersoff의 포텐셜 에너지 함수를 이용하였다. 작용유도 분자동역학 방법의 경우 총 에너지가 보존됨을 알 수 있다. 초기, 말기 상태의 분자 모양을 각각 그림 속에서 표시했다. 복잡한 원자들의 이동경로가 존재하는 경우에는 nudged elastic band 방법으로는 낮은 에너지 장벽을 가지는 경로를 찾기가 힘든 반면 작용유도 분자동역학 방법의 경우에는 적절한 target 에너지를 이용하여 허용되는 동역학적 경로를 계산할 수 있다.  nudged elastic band 방법에서는 주로 전이상태의 원자구조와 이에 따른 최소에너지 장벽을 계산하는 것을 목표로 하는 반면 작용유도 분자동역학에서는 원자들의 동적인 이동경로를 계산하는 것을 목표로 한다. 특히, 위에서 취급한 경우와 같이 복잡한 반응경로 계산에서는 전이상태를 계산하는 것이 간단하지 않을 수 있고 초기경로에 민감하게 전이상태가 계산될 수 있다. 실제로 nudged elastic band 방법의 경우 10 가지의 다른 초기경로를 이용한 것 중에서 가장 낮은 에너지 장벽에 관한 경로를 그림에서 표시했다.

조화 푸리에 구슬(HFB)

   조화 푸리에 구슬(harmonic Fourier beads; HFB) 방법은 아주 최근 Brooks 그룹(Scripps Research Institute)에서 개발한 최소 자유에너지 경로(minimum-free-energy path)를 계산할 수 있는 기법이다. 절대 온도 0 K로 계산된 경로는 최소-에너지 경로(minimum-energy path) 계산으로 볼 수 있다. HFB 방법은 NEB 방법에서 이용하는 스프링을 직접적으로 활용하지 않고, 원자가 특정 좌표 근처에 머물게 하는 추가적인 조화 포텐셜, bias potential을 도입함으로써, 전이경로를 계산을 할 수 있게 설계된 것이다. NEB 방법의 경우 스프링에 의한 힘과 분자 내부에서 만들어진 원자에 작용하는 힘, 이들 둘 다 고려함으로써, 전이경로 (전이상태)를 찾아가는 양식이다. 이 때, 특별히 전이경로 방향에 수직인 성분만 활용한다는 것이 특징이다. 조화진동자 형식의 bias potential 이 추가로 부여한 다음 분자동역학 방법을 통하여, 주어진 온도 하에서 분자 구조 샘플링을 시도한다. 특별히, 이러한 추가 포텐셜 하에서 얻어진 평균위치를 계산함으로써 특정 화학 반응과 관련된 자유에너지 변화량을 계산할 수 있다. 이러한 계산 절차는 Kaestner와 Thiel 박사들이 개발한 Umbrella Integration 방법에 기초하고 있다. 조화진동자 형식의 bias potential 함수와 원자위치의 평균값을 계산함으로써 선택된 반응좌표를 따라서 자유에너지 변화량을 계산할 수 있다.

   HFB 방법에서도 NEB 방법, string 방법에서처럼 목적함수가 정의되어 있지 않다. 앞서 지적한 자유에너지 변화량은 자유에너지를 step index에 대한 미분으로 연결시킬 수 있음을 의미한다. 인접한 이미지들이 스프링으로 연결되어 있지 않다. 그야말로 독립적인 에너지 최적화 또는 분자동역학 (MD) 전산모사를 수행함으로써 전이경로를 계산할 수 있다. 다만, HFB 방법에서는 인접한 image들 사이의 거리를 정확히 같게 되도록 만들어주는 별도의 절차를 요구한다. 결국, 이 부분이 NEB 방법에서 “스프링”이 하는 역할을 대신 한다고 볼 수 있다. 원자들의 집단적인 동일 이동 거리를 요구하는 조건을 제외하고는 HFB 방법은 상당히 물리적이라고 볼 수 있다. string 방법에서도 유사한 절차를 거친다. 이러한 작업이 없으면 모든 image들이 '낮은 에너지'쪽으로 쏠리게 될 것이다. 이와 같은 상황을 [그림8]에서 도식적으로 나타내었다. 이렇게 image들이 "쏠리는 현상"을 방지하기 위한 실질적인 방법은 아래와 같다. 분산된 각 원자의 이동경로를 "직선+사인함수들"의 공식을 활용하여, 사실상 연속 함수처럼 표현한다. 이러한 공식을 활용하면, 인접한 image간의 거리계산을 연속함수 적분하듯이 수행할 수 있다. 이러한 수치적 적분을 이용하면 원자들의 집단 이동거리 계산을 할 수 있다. 초기 상태에서 말기 상태까지의 원자들 집단 거리를 수치적으로 계산하여 그 값을 얻는다. 이 값이 L이라고 하면, L/P값이 결국 인접한 image들 사이의 이동거리가 됨을 의미한다. 특정한 매개변수를 선택함으로써 순차적으로 동일한 이동거리를 확보할 수 있다. “직선+사인함수들” 공식을 화용하면 자유 에너지 계산에서도 마찬가지로 일의 양을 연속함수처럼 표현할 수 있다. 일의 양(자유 에너지), 이동 거리를 적분할 때, 연속함수처럼 적분할 수 있다. 해당 적분들은 1차원 적분으로, 심프슨 공식을 적용할 수 있다. [그림9]에서는 HFB 방법의 실질적인 구현에 사용될 수 있는 순서도를 표시했다. NEB 방법, string 방법, ADMD 방법들과 마찬가지로 HFB 방법도 본질적으로 쉽게 병렬화 될 수 있는 형식의 계산들로 이루어져 있다. [그림10]에서는 인접한 image간의 동일한 거리를 요구하는 모듈의 작업 특성을 입력/출력 형식으로 표시했다. 순서적으로 결정된 동일 이동거리를 제시하는 파라미터도 step index의 함수로 동시에 그림에 나타내었다. 매개변수의 결정을 위해서 Brent 알고리듬을 구현한 zeroin 함수를 사용하였다. 매개변수를 순서적으로 구해내는 계산은 하나의 프로세서에서만 수행하였다. 상대적으로 계산양이 많지 않기 때문에 병렬효율성을 저하시키지는 않음을 확인하였다. [그림11]에서는 HFB 방법을 활용한 실제 계산에서 병렬효율성을 테스트해보았다. alanine dipeptide 이성질체의 변환과정에 대한 전산 모사를 수행하였다. 프로세서의 수에 따른 계산 종료 시간을 측정하여 speedup을 구했다. 물론, 동일한 입력에 대해서 동일한 출력을 얻었다. 포트란90 언어와 MPI 라이버러리들을 활용하여 병렬처리 작업을 수행하였다. 계산 수행에 앞서서 예측한 것처럼 HFB 방법을 활용한 자유 에너지 프로파일 계산은 높은 병렬효율성을 보였다. 동일 이동 거리 조건을 요구하는 작업, 자유 에너지 프로파일을 작성하는 작업들은 첫 번째 프로세서에서만 수행해도 전체 계산의 병렬효율성 저하에는 큰 영향을 미치지 않음을 확인할 수 있었다. 이것은 HFB 방법에서 분자동역학 계산 부분이 가장 많은 컴퓨터 자원을 활용하는 부분이라는 것을 의미한다. 이러한 경향은 예제에서 사용한 분자계보다 더 큰 시스템을 취급할 경우 더욱 뚜렷하게 나타날 것이다. 최근 발표된 단순화된 string 방법에서는 조화 bias potential을 활용하지 않고 단순히 image들을 이완시키는 단계를 string 방법 계산의 첫 단계로 한다. 두 번째 단계에서는 동일한 이동 거리를 요구하는 조건을 사용한다. 두 단계로 이루어진 단순화된 string 계산 방식은 기존의 string 방법에 비해서 안정적으로 전이경로를 계산할 수 있게 해 준다. 알고리듬 상으로 이전의 string 방법보다 간결하다. 첫 단계에서 사용하는 힘은 경로의 수직 성분만을 사용하는 방식이 아니다. 모든 성분의 힘을 활용하여 이완 과정을 수행한다. 동일한 이동 거리를 추가로 요구하는 점은 단순화된 string 방법과 HFB 방법에서 공통으로 활용되는 개념이다. 이 개념을 활용하는 것은 전이 경로 계산 알고리듬을 단순화시킬 뿐만 아니라 안정화시킨다.


그림 8_ HFB 방법에서 일어날 수 있는 수렴 과정들(HFB cycle)을 도식적으로 나타내었다. 조화 bias potential 하에서 이완되는 image(분자구조)들은 에너지가 낮은 쪽으로 이완되기 쉽다. 인접한 image들 사이에서 정의되는 원자들의 집단 이동거리를 인위적으로 동일하게 만들어주는 단계를 도입함으로써 인접한 image들이 동시에 에너지적으로 유리한 경로를 만들 수 있다. 각 HFB cycle 마지막 단계에서 위의 조건을 부여한다. HFB 방법에서 동일이동 거리를 요구하는 단계가 NEB 방법에서 원용되고 있는 스프링의 역할과 비교대상이다. 각원자의 이동거리는 매개변수의 연속함수로 생각하고 수치적분을 수행한다. 실제 인접한 image 사이의 원자 집단 거리계산에서는 심프슨 공식을 사용하여 거리를 계산한다. 심프슨 공식을 활용한 계산에서 전체 이동거리(L)를 계산하여 매 단계가 가져야 할 동일한 이동거리(L/P)를 계산한다. 이 계산에서 알아낸 인접한 image 사이의 동일한 이동거리(L/P)를 보장하는 매개변수를 step index 순서를 따라서 결정한다. zeroin이라는 포트란 루틴을 활용하여 특정 매개변수를 찾아낼 수 있다.

 


그림 9_ HFB 방법에 해당하는 순서도를 나타내었다. HFB 방법의 출력물은 최소 자유 에너지 경로(minimum-free-energy path)와 자유에너지 프로파일(free-energy profile)이다. 온도, ‘연결된 image들’이 입력으로 볼 수 있다. 연결된 image들은 HFB cycle을 통해서 계속 변화될 것이다. 하나의 HFB cycle에서는 분자동역학 계산이 수행된다. 이러한 분자동역학 계산으로부터 원자들의 평균 위치가 계산된다. 이 평균 위치가 다음 HFB cylce에 사용되는 ‘연결된 image들’이 된다. 위에서 선택된 입력을 바탕으로 조화 bias potential 하에서 분자동역학을 수행함으로써 각 원자들에 대한 평균위치를 계산할 수 있다. 조화 bias potential의 강도를 결정하는 변수는 HFB 계산에서 활용되는 또 다른 입력이다. 이러한 평균위치로부터 자유 에너지의 변화량을 계산할 수 있다. 얻어진 자유 에너지 변화량은 해석적으로 얻어진 계산 결과는 아니다. 이는 자유 에너지 프로파일 계산에 사용된다. 각 단계의 HFB 이완 단계 사이에는 인위적으로 인접한 image들 사이의 이동거리를 동등하게 만들어 주는 계산 단계가 있다. HFB 방법에서 얻어진 원자들의 평균위치들을 다음 HFB cycle 에서 기준 원자위치가 된다. 이러한 HFB cycle 이 진행되면서 최소-자유에너지 경로가 결정된다. 각 원자들의 경로가 연속함수로 표시되어질 있다. 이를 활용하고 심프슨 공식을 활용하여 인접한 image들 사이의 이동거리를 계산한다. HFB 방법을 활용한 계산은 상대적으로 쉽게 병렬화 할 수 있다. 사이안 색으로 표시된 점선 부분이 병렬화로 처리될 수 있는 부분이다. 독립적인 분자동역학 계산이 수행된다. 각 원자들의 위치뿐만 아니라 각 원자에 작용하는 평균적인 힘, 조화함수 포텐셜에 기인한 힘을 매개변수의 연속함수로 만들어서 적분공식에 사용하였다. 거리 계산 및 자유에너지 프로파일 계산 관련된 1차원 적분은 심프슨 공식을 활용하였다.


그림 10_ HFB 방법에서 사용되는 인접한 이미지들 사이의 동일한 집단 이동거리 조건을 구현한 하나의 예를 보였다. alanine dipeptide의 이성질체 변환과정을 전산 모사 한 예이다. 동일한 이동거리 조건을 부여하기 전과 후의 원자 집단 이동 거리의 분포를 표시했다. ADMD 방법에서 얻어진 이동경로를 입력으로 사용하여 이미지간 원자 집단 이동 거리를 계산한 결과 인접한 이미지간의 거리가 동일하지 않음을 알 수 있다. 동일한 이동거리 조건을 요구하는 분산된 매개변수를 아래쪽 판넬에 나타내었다. 각원자의 이동경로는 해석적으로 얻어지고 일차원 적분은 심프슨 공식을 이용하여 수치적으로 계산하였다. 주어진 경로가 정확하게 동일한 집단이동 경로를 나타낸다면 매개변수들은 정확히 step index를 변수로 하는 1차 함수를 따라 배열될  것이다. 매개변수의 결정을 위해서 Brent 알고리듬을 구현한 zeroin 함수를 사용하였다. 일차원 수치적분을 수행하기 위해서 심프슨 공식을 이용하였다.

 


그림 11_ 한 가지 예제 계산을 통해서 HFB 방법을 활용한 자유에너지 프로파일 계산의 병렬 효율성를 보여주고 있다. alanine dipeptide의 이성질체 변환과정을 전산 모사 한 예이다. wall-clock 시간 대비 speedup을 사용한 프로세서들의 수에 대해서 나타내었다. 단일 프로세서를 사용할 경우 계산 종료에 소요되는 시간은 13.45 시간이다. 주요 계산 시간을 차지하는 hot spot은 HFB cycle 내에서 분자동역학 계산을 수행하는 부분이다. 이들 분자동역학 계산은 각 프로세서들에서 사실상 독립적으로 계산이 수행된다. 예제 계산에서 사용된 분자모델 보다 크기가 큰 분자모델일수록 병렬효율성 자체는 일반으로 증가한다. 동일한 원자집단 이동거리 조건 계산, 자유에너지 프로파일 계산들은 첫 번째 프로세서에서만 수행된다. 이러한 계산 형식이 전체 병렬 효율성 저하에는 큰 영향을 미치지 않음을 알 수 있다. 이러한 현상은 대부분의 컴퓨터 계산시간이 분자동역학 계산들에 할당되기 때문이다. 난수를 활용한 분자동역학 계산이 수행되지만, 계산 결과는 프로세서들의 숫자에 의존하지 않는다.


 응용 연구들

  탄소60(벅민스터 풀러렌)의 형성 과정은 그 동안의 많은 학문적 관심과 실험적 이론적연구에도 불구하고 잘 알려져 있지 않다. 하지만, 탄소동위 원소를 활용한 실험 등 여러 가지 추가적인 실험적 사실로부터 탄소60의 초기 분자모양에 대한 몇 가지 결정적 단서를 확보하고 있는 상황이다. 탄소60은 완전히 분해된 흑연 판상구조로부터 만들어진다. 즉, 흩어져있던 탄소원자들이 하나하나 뭉쳐서 하나의 탄소60 분자를 형성한다. 이는 탄소 동위원소 12 또는 13으로 만들어진 흑연 샘플을 이용한 실험으로부터 알 수 있는 정보이다. 또한, 당초 연구자들의 예측과 달리, 단일 고리구조를 포함하여 다중 고리구조가 탄소60의 유력한 선행구조가 아니라는 것이 밝혀졌다. 선행구조로부터 탄소60이 형성될 때 필요한 활성화 에너지도 C-C 간의 결합에너지 (약 3 eV) 보다 크지 않다는 것이 실험적으로 알려졌다. 탄소60 분자의 형성과 관련된 시간 스케일은 msec 단위에 머물러있다는 실험적 증거도 있다.

   탄소60의 형성과정에서 유리한 선행분자들의 모양은 체인 모양, 6각형, 5각형 고리들을 포함하는 분자구조임을 작용유도 분자동역학 전산모사를 통해서 알 수 있다. 탄소60으로 변형되는데 유리한 조건의 판단 기준으로 두 가지 조건을 활용했다. 첫 번째 조건은 선행분자와 탄소60과의 에너지 차이가 크지 않아야 한다. 두 번째 조건은 선행분자가 탄소60으로 변형될 때 필요한 활성화 에너지가 높지 않아야 한다는 조건이다. 선행 분자구조가 안정된 탄소60 분자에 비해서 에너지적으로 불안할 경우 설령 탄소60으로 형성된다고 하더라도 발열반응의 특성상 너무 많은 에너지가 운동에너지로 변환되어 구조 자체가 불안해 질 수 있다. 활성화 에너지가 3 eV 보다 과도하게 크지 않은 선행 분자모양에 대한 조건을 활용하고 작용유도 분자동역학 방법을 통하여 선행 분자들의 모양과 변환 이동경로를 계산했다. 탄소 원자 60 개로 구성된 부서진 흑연 단일판 구조는 너무 높은 활성화 에너지를 요구함이 작용유도 분자동역학 방법으로 밝혀졌다. 이전에 제안된 여러 가지 선행 분자구조들은 대개 여러 개의 고리모양들을 구성성분으로 가지고 있는 선행 분자들 이였다. 그 이유는 탄소원자들이 적절한 크기의 고리모양을 유지할 때, 에너지적으로 안정하다는 사실을 활용한 것이다. 얽혀진 여러 가지의 고리모양들로 구성된 3차원 구조 모형의 선행 분자들이 탄소60으로 변환되기가 쉽다는 것을 작용유도 분자동역학 방법을 통해서 알았다. 이러한 선행 분자구조들은 전산모사 된 담금질(simulated annealing) 방법을 활용하여 수치적으로 얻었다. 계산된 선행분자들이 탄소60으로 변환될 수 있는 시간적 비율은 온도 영역 1000-2000 K에서 대략 10-3/sec 정도이다. 이는 실험적으로 알려진 사실과 잘 일치한다.

  총 150 단계를 거쳐서 안정된 탄소60의 분자모양을 찾아가도록 작용유도 분자동역학 전산모사를 설정했다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용하였다. 체인형식의 탄소 원자들의 배열은 에너지적으로 유리하다. 체인 형식의 탄소 원자배열은 체인모양이 휘어지는데 필요한 활성화 에너지가 크지 않기 때문에 선행 분자의 모양 변형에서도 매우 유리한 구조라는 것을 알 수 있다. [그림12]에서는 4가지 선행 분자들의 구조들과 각각이 탄소60으로 변형될 때, 여러 단계마다 포텐셜 에너지 함수의 변화를 표시했다. [그림13], [그림14]에서는 서로 다른 두 가지 선행 분자들에 대해서 분자모양의 변화들을 각각 표시했다.

  선행분자 모형에서 탄소60 분자모형으로 분자가 모양을 바꿀 수 있는 방법은 사실 매우 많이 존재한다. 작용유도 분자동역학 방법을 수행하기 전에 두 분자를 잘 중첩시키는 문제를 풀어야 한다. 이는 한 분자의 회전, 병진 운동들에 대한 자유도를 결정하는 문제이다. 두 구조의 RMSD(root mean square deviation)을 목적함수로 해서 한 구조를 병진, 회전시킴으로써 두 구조사이의 최대 중첩 조건을 찾을 수 있다. 선행분자가 탄소60 구조로 바뀔 때에 각 원자들이 움직여야 하는 거리를 크게 잡지 않는 간단한 하나의 방법은 원자순번을 바꾸어서 RMSD값을 최소화하는 방법이다. 이를 수행하기 위해서 simulated annealing 방법을 동원하여 원자순번을 바꾸는 작업을 수행하였다.  


그림 12_ 탄소60으로 변형되어 가는 과정을 따라서 포텐셜 에너지 함수의 변화를 그림으로 나타내었다. 탄소60으로 변형되기 전의 선행 분자구조들을 표시했다. 이들 선행 분자구조들은 각각 국소적 포텐셜 에너지 함수 최소화 작업을 통해서 얻어진 구조들이다. 분자 구조의 변화가 대해서 활성화 에너지는 대략 C-C 결합의 에너지(약 3 eV) 정도와 비슷하거나 그것 보다 작을 수 있음을 확인할 수 있다. 선행 분자들이 탄소60 구조로 변형되는 과정을 전산모사하기 전에 선행 분자와 탄소60 두 분자들 사이의 적절한 포개짐(중첩; superposition)을 위해서 한 분자구조를 병진, 회전 운동을 시켜서 최대한 잘 포개지도록 만들었다. 이러한 작업은 원자순번에 대해서도 수행하였다. simulated annealing 방법을 활용하여 원자순번을 바꾸었는데 조건은 최대한 두 분자들이 잘 포개지는 조건을 활용하였다. RMSD(root mean square deviation)을 활용하여 두 분자들의 포개진(중첩된) 정도를 측정하였다. Tangled polycyclic I, II의 활성화 에너지 장벽은 각각 3.3 eV, 1.0 eV 이다.


그림 13_ 탄소60이 형성되는 과정을 작용유도 분자동역학 방법으로 구했다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용하였다. 선행 분자구조는 simulated annealing 방법을 통해서 얻어진 구조이며, 국소적 포텐셜 에너지 함수 최적화 과정을 거쳐서 얻어진 국소적으로 안정된 구조이다. 총 150 시간단계를 활용하였다. 반응경로 상에서 활성화 에너지가 과도하게 높아지지 않기 위해서는 에너지적으로 안정하면서 휘어지기 쉬운 체인 모양의 탄소 원자 배열이 부분적으로 선행 분자에 존재하는 것이 유리함을 알 수 있다. 선행 분자구조와 탄소60 분자 사이의 포텐셜 에너지 함수값 차이는 41.94 eV이고 활성화 에너지 장벽은 2.0 eV 이다.



그림 14_ 탄소60이 형성되는 과정을 작용유도 분자동역학 방법으로 구했다. 선행 분자구조는 simulated annealing 방법을 통해서 얻어진 구조이며, 국소적 포텐셜 에너지 최적화 과정을 거쳐서 얻어진 국소적으로 안정된 구조이다. 총 150 시간단계를 활용하였다. 선행 분자구조와 탄소60 분자 사이의 포텐셜 에너지 함수값 차이는 39.55 eV이고 활성화 에너지 장벽은 0.4 eV 이다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용하였다.

  탄소60 분자들 두 개가 하나의 캡슐모양 탄소120으로 변환하는 경우도 희귀사건의 하나로 볼 수 있다. 더욱이 탄소나노튜브 속에 정렬해 있는 탄소60 분자들이 상대적으로 높지 않은 800oC의 온도조건에서 탄소나노튜브로 모양을 바꾸는 실험사실을 설명하기는 매우 어렵고 또한 구체적인 원자 이동경로도 알려져 있지 않다. 탄소60 분자들의 변형과정을 연구하기 위하여 작용유도 분자동역학 방법을 활용할 수 있다. 다중벽 탄소나노튜브는 나노튜브들 간 벽사이 간격을 대략 0.34 nm(흑연 구조에서 단일판 사이의 간격)를 띄우고 성장하는 경향이 있다. 바깥쪽 탄소나노튜브의 직경을 알 경우 안쪽에 형성될 탄소나노튜브의 직경은 거의 정해져 있다는 뜻이다. 안쪽에 원하는 직경의 탄소나노튜브를 성장시킬 수 있는 방법으로 활용될 수 있다. 탄소 원자간 결합의 세기에 비해서 외부 온도조건 800oC는 탄소 원자 간의 결합의 파괴와 재생성이 필수 불가결한 캡슐모양의 분자를 생성하기에는 턱없이 낮은 온도조건이다. 하지만 실험적으로는 열처리만으로도 캡슐 또는 탄소나노튜브가 형성된다는 보고가 있다.

  전자빔, 이온빔을 활용하여 탄소나노튜브들을 용접하여 새로운 나노구조물을 만들 수 있음이 잘 알려져 있다. 많은 결공 결함들이 생기고 탄소간의 재결합을 통해서 ‘X’, ‘T’, ‘Y’ 모양의 탄소나노튜브 유도체가 생성될 수 있다고 알려져 있다. 실험적으로 원자의 이동경로를 확인 할 수 없다. 작용유도 분자동역학 방법을 적용한 결과 약 7 eV 정도의 활성화 에너지는 불가결하다는 결론에 도달했다. 이는 탄소-탄소 원자간 결합의 세기를 에너지로 나타냈을 때, 두 개 정도의 원자간 결합이 필수적으로 끊어지는 상황이 필요함을 의미한다. 하지만, 이러한 에너지 장벽은 실험적 사실을 설명하기에는 부적절하다.

  외부에서 추가로 첨가된 탄소 원자가 캡슐 분자를 만드는데 촉매로 활동할 수 있다는 가정 하에, 탄소 원자를 추가로 도입하여 캡슐이 형성될 수 있음을 작용유도  분자동역학 방법으로 확인 할 수 있다. 이 때, 촉매로 사용된 탄소원자는 반응 후 캡슐 밖에 위치할 수도 있지만, 캡슐의 일부분으로 존재할 수 있다. 물론, 추가된 탄소 원자의 수를 두 개로 확장하기도 했다. 이렇게 될 경우, 충분히 낮은 활성화 에너지를 가질 수 있음을 작용유도 분자동역학 방법으로 보일 수 있다. 추가된 탄소 원자가 캡슐형성에 필수 불가결한 탄소 원자 간 결합 파괴 과정에서 발생할 수밖에 없는 활성화 에너지의 증가를 완화시킬 수 있다. 외부에서 추가된 탄소 원자가 캡슐형성 과정에서 발생하는 불완전한 탄소 간 결합의 상태에 직접 관여하여 원자간 결합의 불안정성을 완화시켜 준다. 이러한 이유로 캡슐이 형성되는 과정에서 7 eV 수준의 활성화 에너지를 필요로 하지 않고도 전체 반응이 완료될 수 있음을 보일 수 있다. 자동촉매 반응의 한 예로 해석할 수 있다.

  [그림15]에서는 외부에서 추가된 탄소 원자의 영향 하에서 탄소60 두 개의 분자가 비교적 활성화 에너지가 낮은 방식으로 탄소120 분자로 변화하는 과정을 표시했다. 각 단계마다, 추가된 탄소원자의 결합상태를 국소 상태밀도를 통해서 분석했다. 추가된 탄소 원자가 캡슐형성의 중간단계에 추가적인 탄소 원자간 결합을 유지하면서 적극 개입하는 것을 확인할 수 있다. [그림16]에서는 추가된 탄소원자의 촉매역할에 의한 실질적인 활성화 에너지의 감소를 표시했다. 추가된 탄소 원자의 수가 한 개, 두 개일 때를 각각 추가된 탄소 원자가 없는 경우와 비교할 수 있다. 총 300 단계로 나누어서 탄소60 분자들의 융합과정을 표시했다. 추가된 탄소 원자가 융합된 캡슐의 일부분이 되고 반응전 탄소60에 있었던 하나의 탄소 원자가 캡슐모양 탄소120 밖으로 나오는 상황을 exchange라고 표시했다. 두 개의 탄소 원자가 위의 exchange 과정을 밟을 경우에도 활성화 에너지가 낮아 질 수 있음을 확인했다. 촉매원자의 활약상은 상황에 따라서 많이 달라질 수 있다고 판단된다. 마찬가지로 exchange 과정에서도 계산에서 확인된 방식이외의 다른 방식으로도 탄소120 분자가 만들어질 수 있다.


그림 15_ 두 개의 탄소60 분자들이 융합해서 캡슐모양의 탄소120 분자를 만드는 과정을 작용유도 분자동역학 방법으로 전산 모사했다. 두 개의 촉매작용을 하는 탄소 원자들을 각각 파랑색과 초록색으로 표시했다. 이들이 관여하여 자동촉매 반응을 성사시킨다. 촉매역할을 수행하는 원자 내의 전자들의 국소상태밀도를 각 작용유도 분자동역학 단계마다 표시했다. 총 300시간단계를 활용하였다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용하였다.

 


그림 16_ 탄소60 분자 두 개가 융합하여 캡슐 모양의 탄소120 분자를 만드는 화학반응을 고려했다. 탄소60 분자의 탄소 원자 간 결합의 파괴가 필수 불가결한 상황에서 반응이 시작되어야 한다. 촉매역할을 하는 탄소 원자가 없는 경우는 활성화 에너지가 높은 반면, 촉매역할을 하는 탄소 원자가 한 개 또는 두 개가 추가될 경우 해당 활성화 에너지가 낮아 질 수 있음을 보였다. 추가된 탄소 원자가 캡슐모양의 탄소120의 구성요소로 변환된 경우를 exchange라고 표시했다. 작용유도 분자동역학 방법에서 총 300 시간단계를 활용하였다. 반경험적 tight-binding 포텐셜 에너지 함수를 이용하였다.

  포스트 지노믹시대에 아주 중요한 문제들 중 하나는 단백질 접힘 과정을 원자수준에서 밝히는 것이다. 단백질 접힘 과정은 단백질 공학, 질병 치료의 근원적인 정보를 포함하고 있는 매우 중요한 문제라고 할 수 있다. 많은 이론적 연구자들이 관심을 가지고 단백질 접힘 문제를 다양한 방법으로 풀고 있다. 아미노산 서열로부터 생성된 단백질이 단백질 고유기능을 수행하기 위해서는 단백질 고유의 3차원 형태로 잘 접혀야한다. 단백질 접힘 문제는 아미노산 서열과 단백질 구조, 기능 사이의 근원적인 원리를 찾는 것을 포함한다. 단백질 접힘과 관련된 응용 연구는 질병과 관련된 신약 개발에 응용될 수 있을 것으로 예측된다. 분자 생물학의 근본적인 목표 중의 하나는 단백질의 기능을 원자수준의 정밀도로 이해하는 것이다. 단백질 접힘의 과정은 다양한 변이가 있고 굉장히 복잡한 것으로 알려져 있다. 장시간에 걸쳐서 일어나는 단백질 접힘 현상은 전형적인 희귀사건으로 볼 수 있다. 단백질 접힘의 근본적인 힘들 중 하나는 수소결합이라고 볼 수 있다. 단백질 분자 내에서 여러 개의 수소결합이 일반적으로 가능하다. 극히 간단한 매우 특별한 단백질의 경우를 제외하고 대부분의 단백질의 접힘 문제에서 반응좌표를 설정하기란 매우 어렵다. 수소결합의 숫자, contact들의 숫자, 단백질의 유효반경 등이 반응좌표로 사용되는 경우가 많다.

  단백질 접힘 문제는 이론적으로 계산하기 상당히 까다로운 경우의 응용문제로 분류될 수 있다. 단백질 접힘에 대한 전산모사가 어려운 이유는 단백질 내에 존재하는 원자수가 매우 많고 매우 복잡한 형식으로 장시간에 걸쳐서 접힘 과정이 진행되기 때문이다. 사실, 원자간 결합이 끊어지거나 하는 극단적인 상황이 일어나지 않는다. 대부분의 단백질의 경우 상당히 큰 분자량을 가지고 있다. 실제로 활성화 에너지가 높지 않다는 점을 고려하면 매우 유연한 운동으로 단백질 분자가 접힐 수밖에 없다. 대부분의 경우 많은 평면각들(dihedral angles)의 총체적인 변화들에 의해서 단백질 접힘이 일어난다고 볼 수 있다.

   작용유도 분자동역학 방법을 활용한 응용연구 내용을 살펴본다. 초기상태와 말기상태를 각각 펼쳐진 단백질 분자, 고유의 기능을 수행하는 단백질 구조로 잡고 단백질이 접혀가는 과정을 작용유도 분자동역학 방법으로 전산모사 하였다. 삼차원 xyz 좌표계를 활용하여 원자들의 위치를 표시했다. 제2차 구조의 중요한 두 가지 기본 단위형식으로 나타나는 헬릭스모양과 병풍모양 단백질들의 접힘 과정을 각각 연구했다. 이들 두 가지 단백질 접힘 과정은 이론 및 실험적으로 많이 연구된 것들이다.

  [그림17]에서와 같이 헬릭스(acetyl-Ala10-N-methyl amide) 형성의 경우, 인접한 아미노산 간의 거리가 3인 경우들로부터, contact 들이 먼저 형성된 다음 헬릭스 고유의 수소결합이 형성되는 것을 확인할 수 있었다. 즉, contact 형성에서는 (i,i+3) - (i,i+4)와 같은 순서로 contact가 형성됨을 알 수 있었다. 이는 관련 실험사실, 다른 이론적 계산결과들과 잘 일치하는 내용이다. 아미노산 10개로 구성된 비교적 간단한 단백질 모형이다. 좀 더 복잡한 단백질에 대한 응용도 수행되었다. 작용유도 분자동역학 방법이 단백질 접힘 문제에 적용될 수 있음을 보여준 계산이다.

  [그림18]에서 보인바와 같이 병풍구조(resiudes 41-56 of protein G)를 가지는 베타-헤어핀에 대한 접힘경로를 계산했다. 이 경우에는 단백질의 가운데 부분이 접히면서 머리핀 모양을 취한다. 작용유도 분자동역학 방법으로 접힘경로를 계산한 결과 단백질의 중앙 부위에서 다양한 형식으로 contact들이 형성되고 이 근처에서 첫 번째 수소결합이 만들어진다는 것을 알았다. 나머지 수소결합들은 중앙으로부터 바깥쪽으로 향하는 순서로 만들어짐을 알아내었다. 이는 실험적 사실과 잘 일치하는 결과이다. 다른 일부 이론계산에서와는 차이를 보이는 경우가 있음을 알았다.
 

그림 17_ 알파 헤릭스 구조(acetyl-Ala10-N-methyl amide)를 가지는 간단한 단백질의 접힘 과정을 작용유도 분자동역학 방법으로 전산모사 했다. 접힘 과정을 따라가면서 단백질 자체의 포텐셜 에너지 함수의 변화를 나타내었다. 실질적인 물분자들을 전산모사에서 사용하지 않고 묵시적인 유효 물분자 효과를 포함하는 포텐셜 함수를 사용하였다. 총 400 시간단계를 이용하였다. 경험적 포텐셜 에너지 함수(tinker 패키지에 있는 amber 포텐셜 함수)를 이용하였다. 수소결합들이 만들어지기 전에 형성되는 contact 들의 형성 순서가 실험 사실과 잘 일치한다.

 

그림 18_ 베타 헤어핀 구조(resiudes 41-56 of protein G)를 가지는 단백질의 접힘 과정을 작용유도 분자동역학 방법으로 연구했다. 접힘 과정을 따라가면서 단백질 자체의 포텐셜 에너지 함수의 변화를 나타내었다. 실질적인 물분자들을 전산모사에서 사용하지 않고 묵시적인 유효 물분자 효과를 포함하는 포텐셜 에너지 함수를 사용하였다. 총 1000 시간단계를 이용하였다. 작용유도 분자동역학 계산에서는 반응좌표를 미리 정의하지 않고 단백질의 접힘 과정을 전산 모사 한다는 특징이 있다. 분자 모형 그림에서 단백질의 펼쳐진 상태(j=0)와 고유의 기능을 수행하는 접힌 상태(j=1000)뿐만 아니라 접힘 과정의 중간 단계를 시간단계 순으로 나타내고 있다. 경험적 포텐셜 에너지 함수(tinker 패키지에 있는 amber 포텐셜 함수)를 이용하였다. 첫 번째 수소결합이 단백질 중간 부분에서 형성됨을 알 수 있다. 나머지 수소결합들은 중간부분에서 바깥쪽 방향으로 나가는 순서로 형성되었다. 이러한 수소결합들의 형성과 관련된 공간적 분포와 시간적 순서는 실험적 사실과 잘 일치하는 것이다.

  단백질 접힘 연구에 적용된 작용유도 분자동역학 방법은 기존의 분자동역학 방법의 응용에서와 달리 아래에서 지적된 몇 가지 가정들로부터 자유롭다. 고온 전산모사 방법의 가정이 없다. 반응좌표의 선택이 없다. 접힌 단백질에 대한 인위적인 판단 기준을 사용하지 않는다. 지수함수 운동학적 확률론을 사용하지 않는다.

  금속 내에 존재하는 전위(dislocation) 결함의 이동도 전형적인 희귀사건의 한 형태로 볼 수 있다. 결정체 금속 물질들의 기계적 성질들은 전위 결함들의 성질과 이동에 의해서 결정되는 경우가 많다. 소성 변형 등은 전위 결함의 운동들에 밀접하게 연관되어 있음이 알려져 있다. 구리 금속 내에 존재하는 cross-slip 문제를 작용유도 분자동역학 방법으로 연구했다. fcc 결정 내에서의 cross-slip은 screw dislocation이 glide 평면을 (111)에서 다른 평면으로 바꾸는 과정을 말하는 것이다. 반도체나 금속의 소성의 이해에 필수적인 원자들의 움직임이라고 볼 수 있다. 원자모델을 사용하지 않고서는 cross-slip 문제를 제대로 풀 수가 없다. 반면, 모든 연구자들이 동의하는 관련 원자모델이 정립되어 있지 않다. 현재 까지 알려진 기존 모형들에 비해서 활성화 에너지가 더 낮은 전위 이동경로를 찾을 수 있었다. 기존의 Friedel-Escaig 미케니즘과 Fleischer 미케니즘의 중간 형식으로 cross-slip 이 이동함을 확인했다. 전위 core에서의 원자들의 집단적 움직임으로 활성화 에너지를 더 낮출 수 있음을 알 수 있다. [그림19]에서 원자모형과 포텐셜 에너지 함수의 변화를 표시했다.



그림 19_ 구리 내의 전위의 모양을 표시했다. 구리 내의 전위의 모양과 전위의 움직임에 따른 활성화 에너지 장벽을 표시했다. 중간 단계의 core 구조들을 동시에 표시했다. core 부분에서의 집단적인 원자들의 움직임으로 발표된 활성화 에너지보다 약 0.5 eV 정도 더 낮은 활성화 에너지 장벽을 작용유도 분자동역학 방법으로 얻었다. 경험적 포텐셜 에너지 함수를 활용하였다. 총 40 시간단계를 이용하였다. 1.5x107 자유도에 대한 최적화 과정을 거친 계산 결과이다.

    

그림 20_ 최근 시간적, 공간적 정밀도의 순서적 체계를 가지는 통합된 계산방식들이 많이 제안되고 있다. 작용유도 분자동역학 방법과 준연속 방법을 혼용한 혼성 모델링 방법의 한 예를 표시했다. 작용유도 분자동역학 방법은 시간, 공간적으로 확장된 다중스케일 계산방식으로도 활용가능하다.

  전이경로 계산의 미래

   분자세계에서 일어나는 희귀사건을 전산모사 할 수 있는 다양한 방법들을 소개했다. 실제 실험 상황에서 나타나는 화학반응들과 연관된 원자들의 집단운동(concerted motion)들은 수 밀리초(ms) 이상의 시간 간격에 걸쳐서 일어나기도 한다. 이들은 원자세계의 희귀사건이라고 할 수 있으며 새로운 계산방법 개발 없이는 다루기 매우 까다로운 문제들로 분류될 수 있다. 최근 이러한 희귀사건들을 기술할 수 있는 계산방법들이 개발되어 응용되기 시작하였다. 여러 연구자들의 노력으로 과거에는 다룰 수 없었던 화학반응의 구체적인 원자 재배열 단계들을 이론적으로 설명할 수 있게 되었다. 또한, 과거에는 다룰 수 없었던 장시간에 걸쳐 일어나는 여러 가지 실험적 사실들을 이론적으로 계산할 수 있음을 보였다.

  새로운 계산방법들 중에는 기존의 분자동역학 방법을 원용한 것도 있으며  수학적으로 새롭게 정립된 것들도 있다. 예를 들면, hyperdynamics 방법, parallel replica 방법, temperature accelerated dynamics 방법, NEB 방법, string 방법, 작용유도 분자동역학 방법, 조화 푸리에 구슬 방법 등을 들 수 있다. 앞에 기술한 세 가지 방법은 분자동역학 방법에서의 초기조건 문제 특성을 그대로 가지면서도 특정한 희귀사건들을 효율적으로 계산할 수 있다. 희귀사건을 짧은 전산모사 시간 안에서 관찰할 수 있도록 특수하게 처리된 방법이다. 반면 뒤에 열거한 네 가지 방법들은 초기상태와 말기상태로 대표되는 특정 희귀사건이 일어나는 주요 경로를 찾아내는 방법이다. NEB 방법, string 방법, Dimer 방법, ART 방법들은 원자들의 구체적인 이동경로보다는 전이상태를 찾는 것에 초점을 두고 있으며 가장 낮은 에너지 장벽을 계산할 때 유용하게 사용할 수 있는 방법이다. 매우 복잡한 화학반응의 경우 작용유도 분자동역학 방법을 활용하면 에너지적으로 가능한 낮은 활성화 에너지와 관련된 전이경로를 적극적으로 찾을 수 있다. 초기경로에 덜 민감하게 실현가능성이 높은 전이경로를 찾는데 있어서 유리하다. Dimer 방법 또는 ART 방법과 결합된 kinetic Monte Carlo 방법은 미리 정해둔 사건테이블을 사용하지 않고 운동학적 원자들의 재배열을 시간의 함수로 연구할 수 있다는 점에서 매우 의미 있는 이론적 발전으로 해석된다. 앞으로 많은 응용문제 풀이가 뒤따를 것으로 예측된다.

  본 논설에서는 특별히 작용유도 분자동역학 방법이 다양한 분야에서 응용될 수 있음을 보였다. 분산된 작용함수를 대체하기 위해서 분산된 원자 이동경로를 활용한 작용함수에 벌칙항들(총에너지보존, 운동에너지 조절)을 도입함으로써 수치적으로 편리하게 원자들의 이동경로를 얻을 수 있다. 작용유도 분자동역학에서 얻어진 원자 이동경로들은 Verlet 경로와 비교하여 그 정밀도를 정량적으로 평가할 수도 있다. 작용유도 분자동역학 방법에서는 통상의 분자동역학 방법과 마찬가지로 포텐셜 에너지 함수 그리고 원자에 작용하는 힘을 계산함으로써 전이경로의 결정이 가능하다. 병렬효율성이 높은 특징을 가지고 있는 작용유도 분자동역학 방법은 다양한 분야에서 거대계와 관련된 응용문제 풀이에 적용될 수 있음을 보였다.

  [그림20]에서 볼 수 있듯이 작용유도 분자동역학 방법은 공간적, 시간적으로 확장된 다중 스케일 계산방식으로도 확장될 수 있을 것으로 예측된다. 작용유도 분자동역학 방법의 최대 특징은 화학반응에 대해서 미리 반응좌표를 선택하지 않고도 원자 재배열과 관련된 전이경로를 계산할 수 있다는 점이다. 일반적인 원자 재배열 과정 연구에 적용될 수 있을 것으로 기대된다.  유한온도 string 방법 또는 조화 푸리에 구슬 방법을 활용할 경우 최소-자유에너지 경로와 자유에너지 프로파일을 계산할 있다. 이들 방법을 활용한 많은 응용연구가 있을 것으로 판단된다. 이것은 'chain of state' 형식의 모델로 최소-에너지 경로, 최소-자유에너지 경로 모델 연구가 가능함을 의미하는 것이다.

  작용유도 분자동역학 방법을 이용하여 탄소60의 형성, 자동촉매 작용으로 형성 되는 탄소60의 융합, 촉매원자의 역할, 단백질의 접힘, 금속 내에서의 전위의 움직임 등의 응용문제를 소개했다. 작용유도 분자동역학은 기존의 분자동역학 방법이 그 전산모사 시간의 한계에 부딪히는 경우에 광범위하게 적용될 수 있는 계산방법이다. 기존의 NEB 방법, string 방법, Dimer 방법 등과의 연관성 등에 대해서도 토의했다. 향후, 작용유도 분자동역학 방법을 활용한 응용연구들이 나노물질들의 생성, 변형, 그리고 기능을 이해하는데 큰 도움을 줄 것으로 기대한다.

  분자세계로 가면 다양한 형식의 전이경로들이 있다. 뿐만 아니라 전이경로들을 통해서 연결되는 새로운 나노물질들과 새로운 원자재배열과 연관된 프로세스들이 있을 수 있다. 파인만 박사는 일찍이 저 바닥으로 가면 새로운 세상이 있음을 알고 있었다. 사실 따지고 보면, 분자세계로 가면 엄청나게 많은 사건들이 즐비하다. 원자들은 초당 1012 번 이상 진동 할 수도 있다. 분자세계는 그 자체가 엄청난 가능성의 시험무대이다. 분자세계에는 파인만 박사가 지적한 공간뿐만 아니라 엄청난 시간들이 있다고 할 수 있다. 분자세계의 전이경로를 찾아서 여행을 떠날 때가 되었다.


pdf 형식으로 정리된것:
 http://incredible.egloos.com/3947906

cf. 작용유도 분자동역학 방법 (action-derived molecular dynamics) :
http://incredible.egloos.com/372174

cf. 분자동역학 (molecular dynamics) 입문 :
http://incredible.egloos.com/1679232
cf. 자유에너지 프로파일 계산 :
http://incredible.egloos.com/2342691
cf. HFB 방법에서도 동일한 이동 거리 조건을 활용하였다 : Ilja V. Khavrutskii, C. L. Brooks III 의 웹페이지와 프리프린트:
http://mccammon.ucsd.edu/~ikhavru/
cf. Harmonic Fourier Beads 방법, gradient augmented HFB 방법:
http://incredible.egloos.com/3068649

Action-CSA 방법:




핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax