분자역학 및 분자동역학 첫걸음

이 글은 제가 석사 1학년 때인 2000년 2학기 때, 박수영 교수님의 <고분자 전자재료> 과목 기말 레포트 용으로 제출했던 글 가운데 일부를 발췌하여, 화학부 물리화학실험과 일반화학실험 가운데 전산화학과 관련된 내용을 소개할 때 사용했던 글입니다. 잘못된 부분도 있고, 이제는 필요없는 부분도 있지만 더이상 볼 사람도 없고 지금와서 뜯어 고치는 게 별 의미도 없다고 생각해서 방치하고 있습니다. 우연히라도 방문하신 분 가운데, 글에서 눈에 거슬리는 내용이 보이시면 부담갖지 마시고 직접 고쳐 주시면 됩니다. 그러라고 위키로 페이지를 만든 거니까요 :-)

1 들어가며

컴퓨터의 계산 속도가 이제는 글자 그대로 하루 하루가 다르게 빨라지고 있다. 가격 대 성능 비 역시 거의 매일 상승하고 있다 해도 과언이 아닌 바, 개인용 컴퓨터와 워크스테이션의 경계가 무너진 지도 이미 오래이다. 이러한 실질적인 주변 여건이 보다 많은 사람들로 하여금 컴퓨터를 이용한 연구에 눈을 돌리게 하고 있다. 이런 경향은 화학이나 분자 생물학 분야에서 특히 두드러지고 있는데, 이것은 우연이라기보다는 이들 학문이 다루고 있는 계의 특성이 컴퓨터를 활용하기에 적당한 점들이 많기 때문이다. 예를 들어, 고분자가 용액 속에서 서로 꼬이는 현상이라든지, 단백질이 물 속에 들어가서 보이는 구조적 변화 등의 문제를 다룬다고 생각해보자. 이들 문제는 적절한 가정과 직관, 그리고 그에 근거한 근사법을 이용해 단순하고 반복적인 덧셈・뺄셈 또는 미적분의 연속으로 변형할 수 있다. 이렇게 변형된 문제는 계산 량이 손으로 처리할 수 있는 수준을 훨씬 넘어서는 방대한 것이기는 하지만, 결국은 컴퓨터를 이용한 작업에 최적화된 종류의 것이다. 그리하여, 컴퓨터를 이용한 분자 모델링 및 모의 실험은 단순 반복적인 계산을 효율적으로 수행하기 위한 알고리듬의 개발과 적절한 가정 및 근사법의 개발, 그리고 응용이라는 측면에서 발전되어 왔다. 알고리듬 개발의 경우는 계산 물리학이 한창 발전을 하던 60년대에서 80년대까지 이루어진 성과물들이 이미 이 분야의 고전으로 확립된 상황이기는 하지만, 컴퓨터의 성능 향상과 더불어 다룰 수 있는 계의 범위가 확장되면서 그에 적합한 메모리 관리나 새로운 개념을 도입한 프로그래밍 기법의 개발로 이어지고 있다. 흔히 사용되는 수학 라이브러리들도 이 분야의 성과물이라 하겠다. 하지만, 이는 상당 부분 전산과학의 영역이기 때문에 전산 전공자가 아닌 사람들이 이 분야에서 보이는 활약은 그리 두드러지지 못하다. 사실, 많은 사람들은 다루고자 하는 계에 적합한 가설 및 근사법의 개발이나 또는 기존의 물리・화학적 개념들을 컴퓨터를 이용한 모의 실험에 적합하게 변용 하는 일에 노력을 경주하고 있다. 또한, 이미 개발된 근사법과 그것을 구현한 도구들을 어떠한 시스템에 어떻게 응용할 수 있을 것인가에 대해서도 많은 연구를 진행하고 있다.

물론, 분자 모델링 및 모의 실험(혹은 통칭하여 계산화학)이라는 영역이 아직은 실험 및 실험 데이터와 독립하여 존재할 수 있을 정도로 내적 완결성을 갖추지는 못한 것이 사실이다. 그러나, 이러한 이론적 작업을 통해 실험으로부터 얻은 결과물을 우리가 이해할 수 있는 물리・화학적 언어로 번역할 수 있게 되고 나아가 이를 기반으로 실험적으로 얻기 힘든 결과를 예측할 수도 있다. 게다가, 간단한 경험 법칙과는 비교할 수 없을 정도로 정확하게, 실험 이전에 이미 어느 정도의 물성을 예측할 수도 있게 되었다. 그리하여, 간단한 이론적 도구들은 이제 실험을 하는 사람들에게도 역시 필수적인 도구라고 하겠다. 이러한 맥락에서, 이 글에서는 알고리듬 개발이 아닌, 실질적인 화학계(chemical system)에 계산 화학을 응용하는 일과 그에 관련된 문헌을 읽기 위해 필요한 기본 개념들, 그 가운데에서도 분자 역학(및 분자 동역학)과 관련된 기본 용어들을 간략하게 설명하고 이들 방법이 실제로 응용되고 있는 측면에 관해 간단히 소개하고자 한다.

2 분자 모델링 및 모의 실험

우선, 이 모델(model)이라는 말이, 다른 영역에서 일하는 사람에게는 다른 방식으로 이해될 수도 있겠다. 예를 들어, 기계공학에서 ‘물’을 다룰 때 사용하는 모델은 나비에-스토크스 방정식을 만족하는 일종의 연속체일 것이나, 분자 모델링에서 사용하는 ‘물’은 TIP3P나 TIP4P와 같이 특정 모양과 쌍극자 모멘트와 같은 특성을 갖고 서로서로 떨어져 있는 구조물이다. 이 글에서 사용하게 될 “분자 모델링”은 분자라는 말에서도 알 수 있듯이 원자나 분자 수준에서 대상을 모형화하는 쪽임을 쉽게 알 수 있을 것이다. 하지만, 고분자가 물 속에 들어가 있는 시스템과 같은 경우에는 물을 연속적인 유체로 보고(mesoscale) 고분자는 개별 분자로 모델링하는(microscale) 등, 두 영역을 혼용하는 경우도 많다. 직감적으로 받아 들이면 될 듯 싶다.

여느 학문 분야와 마찬가지로, 분자 모델링 및 모의 실험에서도 다루고자 하는 시스템과 보고 싶은 성질에 따라 크게 구분되는 방법론이 존재한다.

이 두가지 커다란 분류체계는 Lipkowitz, Applications of Computational Chemistry to the Study of Cyclodextrins, Chemical Reviews, 1998, vol. 98, No. 5, p. 1830에서 언급된 내용을 참고한 것이다.

그 하나는 QSPR(Quantitative Structure Property Relationships)이나 QSAR(Quantitative Structure Activity Relationships)에서처럼, 분자 구조와 물리화학적 성질 혹은 생물학적 반응간의 연계를 기존의 데이터에서부터 외삽이나 내삽 등의 통계적 방법론을 이용해 새로운 물질의 성질을 유추하여 찾아가는 방법으로, 이는 요즈음 관심의 대상이 되고 있는 신약 개발이나 bioinformatics의 기본 아이디어를 이루고 있다. 다른 하나는, 우리가 알고 있는 물리 이론들을 처음부터 적용해 나가며 관심 있는 시스템을 다루는 방법이다. 여기에서 우리가 알고 있는 물리 이론이란 곧 전자기학을 포함한, 고전 역학과 양자역학으로 대별되는 물리학이다.(이 글에서는 이 후자에 관해서만 간단히 정리해 보기로 한다.)

이 후자의 영역도 크게 두 분야로 나누어 발전해 왔다.

이하 이어지는 문단은 정확한 문헌을 통해 확인하지는 않은 지극히 개인적인 판단이나, 나름대로 유효하다고 생각되어 전체적인 글도 이러한 생각에 근거하여 기술할 것이다.

한 영역은 관심 있는 (분자)계가 갖는 구조와 그에 수반되는 전기적 혹은 구조적 특성처럼 시간 개념이 명확하게 개입되지 않는 (분자)계의 특성, 즉 평형 성질을 연구하는 것이고 다른 한 영역은 시간에 따른 계의 변화를 연구하는 것이다. 이들 각 영역은 양자 역학이나 고전 역학 가운데 관심 있는 계에 적합한 방법을 택해 나름의 방법론을 발전시켜 왔다. 여기에서 ‘나름의 방법론’이란, 대상을 바라보는 물리학의 체계와는 별개의 관점에서, 내적 기술을 위해 자신의 영역을 구축해 온 분자 역학(MM), 분자 동역학(MD), 몬테 카를로 방법(MC) 등의 다양한 방법론을 의미한다. 이러한 방법론적 분류는 단순화시켜 얘기하자면 위상 공간(phase space) 속에서 계를 기술하는 궤적(trajectory)을 움직이게 하는 방식에 따른 분류라고 하겠다. 그런데, MM・MD・MC 등이 대부분 고전 역학에 근거한 PES(Potential Energy Surface)를 중심으로 초기 발전과 응용이 이루어 졌기 때문에 많은 경우 이들 방법은 종종 고전 역학에 근거한 방법의 하위 개념으로 뭉뚱그려져 이야기된다.

이러한 방법들을 이용하면 우리는 진공 상태에 있는 분자 하나의 성질을 자세히 알 수 있다. 예를 들면, 생성열과 같은 열역학적 성질, 쌍극자 모멘트 값을 포함한 전기적 모멘트 값 등을 구할 수가 있다. 하지만, 신약 개발이나 여러 가지 물성 연구에서처럼 분자 하나의 성질보다는 이들 분자가 모여서 어떤 집합적 성질을 보이는 가도 매우 흥미 있는 영역이다. 우리가 관심을 갖는 대상의 특성이 성분 분자 개개의 고유한 특성들에 의존하는 것이 사실이지만, 동시에 이들 분자가 어떤 식으로 모여서 배열하게 되는 가도 물성에 지대한 영향을 미치기 때문이다. 즉, 계산 화학의 목적은 개별 분자의 특성을 이론적으로 연구하는 동시에 개별 분자가 여러 개 모여 있을 때 이들이 보여줄 수 있는 성질을 선험적으로 예측하는 데 있다고 하겠다. 보트(Bott)는 이를 “원하는 물성을 지닌 물질을 분자 수준에서 설계(design at the molecular level)하는 것”이라 표현했다.

2.1 양자 역학에 근거한 방법

양자 역학에 근거한 방법은 주로 (분자)계의 전기적 구조를 정확하게 계산하는 일에 사용된다.

물론, 양자역학을 이용한 방법이 분자의 전기적 구조를 연구하는 일에만 사용되는 것은 아니다. 90년대 들어와서 관심의 초점이 되고 있는 First Principle Molecular Dynamics나 QM/MM은 그 좋은 예라고 할 수 있다. 하지만, 아직까지는 양자역학을 이용한 방법의 대부분의 활용 예가 주어진 분자 시스템의 구조를 최적화하거나 주어진 구조가 갖는 에너지 혹은 그와 관련된 전기적 성질을 연구하는 쪽이므로 여기에서도 간단하게 이 부분만을 이야기한다.

여기에는 사용되는 근사의 정도에 따라 ab initio법과 반 경험적 방법(semi-empirical method)으로 대별되는 방법론이 존재한다. 그러나, 기본적인 아이디어는 다음과 같이 동일한 곳에서 출발한다 : 양자 역학을 이용하면 관심 있는 계의 전자・핵의 공간적 위치를 기술할 수 있다. 대개의 경우 이 작업은 본-오펜하이머 근사 하에서 분자 궤도 함수법을 사용하여 구현된다. 즉, 본-오펜하이머 근사를 따라 고정된 핵의 전자기장 속에서 전자의 움직임이 자기 모순 없는 장(SCF : self-consistent field)을 만족할 때까지 전자에 관련된 슈뢰딩거 방정식을 계산한 다음 핵을 약간 움직인 후 다시 전체 에너지가 자기 모순 없는 장을 만족하도록 반복하여 계산하고, 또다시 전자에 대한 과정부터 반복하는 것이다. 이렇게 해서 더 이상 분자 시스템의 에너지가 감소하지 않으면 우리는 주어진 분자가 가질 수 있는 가장 안정한 전기적 구조를 얻었다고 얘기할 수 있다. 여기에서, ab initio법은 기본적인 물리 상수와 계산을 간단히 하기 위한 이론적 근사 몇 가지만으로 처음부터 끝까지 수치적 혹은 해석적 계산을 통해 원하는 값을 얻는 방식인데 반해, 반 경험적 방법은 ab initio법의 이론 전개 구조를 따라 가면서 나타나는 수많은 적분들을 실험에서 얻을 수 있는 수치들로 파라메터화함으로서 계산 량을 대폭 줄이는 방식으로 구현된다.

주로 사용되는 반경험적 방법인 CNDO, INDO, MNDO, AM1, PM3 등의 구별은 결국 특정 시스템에 대한 적분치 파라메터들을 어떠한 방식으로 최적화했는가에 따라 구별되는 것이다.

그리하여, ab initio법을 이용하면 분자의 구조나 생성열 외에도 전이 상태나 들뜬 상태에 대한 정보도 얻을 수가 있지만 계산 량이 엄청나게 많아지는 데 반해(∝n4), 반 경험적 방법을 사용하면 일정 수준의 정보를 적은 계산 량(∝n3)으로 얻을 수가 있다.

여기서의 n은 기저 함수의 수이다.

2.2 고전 역학에 근거한 방법

고전역학을 이용한 분자 모델링 및 모의 실험에서는, 계가 고전적 입자들이 고전 역학적 포텐셜에 의해 모여서 이루어진 것으로 간주되며 각각의 입자는 고전 역학으로부터 계산되는 힘에 따라 움직이도록 모형화 되어있다. 여기에는 계의 정적인 성질을 연구하는 데 초점이 맞추어진 분자 역학 외에도 시간에 따른 계의 움직임을 연구할 수도 있는 분자 동역학이 매우 큰 부분을 차지하고 있다.

그리고 몬테 카를로 방법도 가장 많이 사용되는 방법론 가운데 하나이다. 최근에 와서는 양자 역학 대 고전 역학이라는 분류 기준과, 분자 동역학 대 몬테 카를로 방법이라는 분류 기준이 서로 다른 관점이라는 사실이 명확해 지고 있지만, 앞에서 언급한 내용과 마찬가지로, 그 초기 발전 양상에서 분자 동역학이나 몬테 카를로 방법이 주로 고전 역학을 이용한 분자 모델링으로부터 발전되어 온 관계로 고전 역학에 근거한 방법론을 이야기할 때 분자 동역학과 몬테 카를로 방법을 그 하위 분류 개념으로 사용하는 경향이 있다. 이 글에서도 MD를 고전 역학에 근거한 모델링의 방법에 포함시켰지만, 이는 순전히 편의상의 분류이다. 분자 동역학의 경우는 따로 장을 두고 기본 개념들을 이야기할 것이나, 몬테 카를로 방법은 발표 주제에서 벗어나므로 이 글에서 다루지는 않을 것이다.

2.2.1 분자 역학 : MM : Molecular Mechanics

분자 역학(MM)은 고전 물리학을 이용하여 관심 있는 (분자)계의 구조, 에너지 등의 성질을 계산하는 방법이다. 이 방법에서는 경험적 힘장(EFF ; Empirical Force Field)이라는 개념이 등장하게 되는데, 이는 분자의 PES를 재현하기 위해 여러 가지 간단한 수학적 함수들을 합해서 표현되는 일종의 포텐셜 에너지 함수이다. 분자 역학에서는 분자를 (물리적으로는, 전자에 기인하는) 일종의 탄성력으로 묶여 있는 입자들의 집합으로 간주하는데, 이 탄성력은 흔히 내부 좌표계의 함수로 표현된 EFF의 미분으로부터 계산된다. 이렇게 입자들의 위치와 각각의 입자에 미치는 힘이 계산되면, 입자들은 구조 최적화 즉 에너지 최소화 과정을 거쳐 좀 더 안정한 공간적 배치를 찾아가게 되는 것이다. 또한, FF의 함수 꼴로 나뉘어 표현된 각각의 에너지 항들이 그 구조에 어떤 정도로 기여하고 있는 지도 분석할 수 있다. 이는 곧 관심 있는 계의 현상을 주도하는 물리적 연원을 밝히는 일인데, 이를 통해 원하는 물성을 얻기 위한 ‘분자 수준에서의 설계’의 첫걸음을 뗄 수 있는 것이다. 하지만, 전형적으로 사용되는 MM의 EFF는 전자에 관련된 정보를 명확한 함수 형태로 가지고 있지 않으므로 화학 결합의 생성과 분해처럼 전자의 움직임과 밀접한 관계가 있는 경우를 기술하기 위해서는 다루고자 하는 계에 적합한 특수한 함수 형태를 포함하는 EFF를 도입해야 한다.

사실 Molecular Mechanics라는 방법론의 활용보다는, Molecular Mechanics Approximation이라는 개념이 더 유용/중요하다고 생각한다. 우리가 사용하는 고전역학 모의실험 프로그램들 거의 대부분에서 접하는 고전적 힘장은, 개념적으로 서로 구분지을 수 있는 분자 내부의 운동들을 모사하는 각각의 함수 항들 - bond/vibration term, angle term, dihedral term, improper term, vdW term, Coulomb term, 등등 - 을 독립적으로 구축하고, 계의 에너지를 이들 함수의 합으로 표현하고 있다. 이 자체가 일종의 근사일 수밖에 없는데, 이렇게 분자계의 에너지를 분자 내부 운동의 normal mode와 energy scale에 근거해 독립적으로 구축한 항들의 합으로 표현하는 방식을 분자역학근사라고 한다. 이는 매우 자연스럽고 당연한 것으로 보이지만, 조금만 생각해보면 이게 정말 보기만큼 당연한 것인지 확신하기가 쉽지 않다 -_-;; (적어도 본인은 이런 근사가 아주 작위적이라고 생각해서, 물리화학 학점이 개판 5분전이다 -_-;;)

2.2.2 분자 동역학 : MD : Molecular Dynamics

분자 동역학에서는 MM에서의 EFF가 만드는 위상 공간 속에서 (분자)계가 (뉴턴) 운동 방정식을 따라 어떤 궤적을 그리는 지를 관찰하게 된다. 물론 여기에서도 입자를 움직이게 하는 힘은 EFF의 미분으로부터 얻어진다. 그러나 MM과는 달리 MD에서는 운동 방정식을 적분함으로써 입자들의 다음 위치와 속도를 계산한다. 그렇게 얻은 위상 공간 속에서의 궤적(혹은 단순히 입자들의 좌표와 운동량)을 통계 역학적 지식을 동원해 해석함으로써 계에 대한 여러 가지 정보를 얻을 수 있다. 또한, 입자들의 위치와 속도에 관한 정보를 구함으로써 QM이나MM과는 달리, 포텐셜 에너지뿐만 아니라 운동 에너지까지도 같이 생각할 수 있게 된다. 따라서, 운동 에너지가 존재하기 때문에 우리는 (분자)계가 만든 PES를 시작 위치의 국지적 에너지 최소치 부근만이 아니라 좀 더 다양한 부분의 PES를 돌아다닐 수 있게 되므로, 단순한 MM과 비교해 볼 때 좀 더 안정한 계의 상태/구조를 얻을 수가 있게 된다. (물론, 고전적인 MD만으로 안정한 구조를 찾는 일은 매우 비효율적이다.) 따라서 단순히 분자 한 두 개로 이루어진 계라면 MM을 이용해 최적화된 구조를 찾을 수도 있겠지만, 여러 분자들로 이루어진 계의 평형 구조를 보고 싶을 때에는 MD를 사용하기도 한다. 근자에는 평형 상태에서 계가 보여주는 구조를 찾는 작업보다는 시간에 따른 계의 변화 자체를 관찰하는 일에 MD를 더 많이 사용한다.

이제까지 MM과 MD에 대해 간략하게 정리해 보았다. 이제 이들 각각의 방법론에서 자주 사용되는 개념과 방법들에 대해 좀 더 자세히 살펴보고, 실제로 이들 방법이 사용되는 몇 가지 예들을 소개하고자 한다.

3 분자 역학 (Molecular Mechanics)

이 부분은 Goodman, 『Chemical Applications of Molecular Modelling』과 Gelin, 『Molecular modeling of Polymer Structures and Properties』을 중심으로, Jensen, 『Introduction to Computational Chemistry』의 내용 가운데 일부를 첨가하여 정리한 것이다. 그 외, 개인적 기억에 의존하지 않고 특별히 인용한 자료는 각주로 처리하겠다.

3.1 MM의 실제

흔히 MM은 0 K에서의 MD라는 식으로 간단히 넘어가는 경우가 많은데 하는데, 사실 MM의 발전되어 온 주된 동기/방향은 분자의 정확한 구조와 에너지를 선험적으로 계산하는 것이다. 이 방법론에서는 오로지 핵의 움직임만이 관심 대상이 되며, 그 핵은 경험적 힘장(이하 FF)에 의해 그 움직임을 지배받는다. 이는 우리가 양자 화학에서 가장 먼저 배우게 되는 본-오펜하이머 근사에 의한 것으로, 아주 들뜬 상태만 아니라면 전자의 움직임이 핵의 움직임보다 아주 빠르기 때문에 핵은 언제나 평형 상태에 이른 전자의 평균적인 영향력 하에서만 움직이게 된다. 그러므로, 대부분의 경우 FF에서는 전자가 미치는 영향이 명확한 함수 형태로 나타나기보다는, FF를 표현하는 여러 함수 형태에 간접적으로만 고려되는데 이것이 결국 평형 상태에서 시간에 대한 평균값에 해당한다.

MM은 기본적으로 FF의 이식성(transferability)을 하나의 가정으로 받아들인다. 여기서 이식성이라 함은, 일관적으로 만들어진 파라메터와 FF의 함수 형태가 FF를 구축하기 위해 데이터를 제공한 분자들뿐만 아니라 전혀 다른 분자 계에도 만족스럽게 적용된다는 이야기이다. 사실, 이는 MOT와 같은 양자 화학적 도구에 의해 검증해 볼 때 그리 자명한 것은 아니다. 양자 화학적인 관점에서 볼 때, 화학 결합이란 분자 내 각각의 원자로부터 전자가 분자 전체에 걸쳐 재 분포되며 발생하는 현상이기 때문이다. 따라서, 각각의 원자는 분자에서 마다 다를 수밖에 없으며 따라서 원자나 기능기를 단위로 하여 분자의 성질을 분석하는 것은 타당하지 않다는 결론을 얻을 수 있다. 하지만, 핵의 움직임과 관련된 사항만을 생각할 경우에는 FF의 이식성이라는 가정도 나름대로 타당성을 가질 수가 있다. 실험 결과도 이러한 접근 방법이 상당히 타당함을 암시한다. 예를 들어, 우리는 적외선 분광학을 통해 분자의 고유한 지문 영역을 관찰할 수 있고 이 사실은 개별 분자의 유일성을 암시한다. 하지만 동시에 적외선 스펙트럼에서 특정한 기능기들이 나타나는 영역은 여느 분자에서나 어느 정도 일정하게 나타난다. NMR에서도 마찬가지로 생각할 수 있다.

위에서 언급한 이식성이라는 가정을 전제하면, 이제 우리는 실제로 MM을 사용하여 분자의 에너지를 구하게 된다. MM에서는 운동 에너지라는 개념이 도입되지 않기 때문에 우리는 오직 FF에 의한 분자 계의 에너지만을 구하게 되는데, 우리는 이것을 분자 역학 에너지(E_MM, 혹은 steric energy)라고 부른다. 그런데, 이 값은 실험하는 사람들이 흥미를 갖고 있는 자유 에너지가 아니다.

일반적으로 자유 에너지(G)는 분자 모델링이나 모의 실험을 통해 구하기가 쉽지 않은데, 이는 분자의 엔트로피를 계산하는 데에 구조적으로 많은 어려움이 존재하기 때문이다. 하지만, 실험하는 사람의 입장에서는 자유 에너지가 바로 평형상수와 직결된 값이므로 매우 중요하다.

물론, E_MM은 내부 에너지와도 값이 다르지만, 이는 영점을 어떻게 설정하는가의 문제이고 실질적으로 Del_E_MM과 Del_U의 값은 동일하다. 그렇다면 자유 에너지를 구하기 위해 어떠한 방법을 사용해야 할 것인가? 물론 FEP나 MD를 이용하여 이론적으로 구할 수도 있지만, 많은 경우에는 E_MM을 자유 에너지에 대한 근사값으로 사용한다. 이것은 타당한 근사가 될 수 있을까? 이제 H = U + PV라는 간단한 열역학적 관계식을 생각해보자. 일정한 압력 하에서라면 이 식은 Del_H = Del_U + P*Del_V가 된다. 그런데, 간단한 MM에서는 외부 압력이 0이므로 결국 우리는 Del_H ~ Del_U ~ Del_E_MM라는 관계식을 얻을 수 있다. 그러므로 열역학의 기본 관계식인 Del_G = Del_H - T*Del_S 에서 오직 엔트로피 변화만이 문제가 되고 있다. 하지만, 많은 MM 계산에 있어 분자 계의 엔트로피 변화는 엔탈피 변화에 비해 그 크기가 매우 적다. 따라서, Del_S ~ 0 이라는 가정이 타당한 계에 대해서라면 다음과 같이 근사를 한다고 해도 크게 문제되지 않을 것이다 : Del_G ~ Del_H ~ Del_U ~ Del_E_MM 그러므로, MM의 대략적인 과정은 EFF의 구축그에 따른 에너지 계산결과의 분석 및 해석이라는 과정으로 이루어진다고 할 수 있을 것이다.

그런데, 고립계라면 계의 총 에너지는 불변이어야 하므로 포텐셜 에너지가 적은 분자들은 운동 에너지 값이 증가하게 된다. 그렇다면, 실제로 개별 분자가 좀 더 낮은 포텐셜 에너지의 구조를 선호하는 것을 어떻게 설명할 수 있을까? 그것은 계의 총 에너지에는 변화가 없지만 엔트로피가 증가하기 때문이다. 즉, 계는 가장 큰 엔트로피에 해당하는 구조를 취하게 되고 그 때 분자는 낮은 을 갖는데 이는 곧 분자의 움직이는 정도가 큰 상황에 해당한다. (물론 지금 이야기하고 있는 에너지는, 앞서 언급했듯이, 계의 자유 에너지가 아니다.) 즉, 실제 분자는 가장 낮은 포텐셜 에너지에 해당하는 구조를 취할 확률이 가장 높다고 말할 수 있다.

주어진 PES 상에서 최소점을 찾는 기법들은 이미 수학자들에 의해 다양하게 개발되어 있으므로 우리는 FF로 주어지는 계의 PES를 따라 움직이며 최소점에 해당하는 분자의 모양을 알아내기만 하면 된다. 그러나, FF와 같은 어떤 다변수 함수의 전체 최소점을 찾는 일은 거의 불가능한 일임이 알려져 있다. 그러므로, 우리가 PES의 어느 점에서부터 출발하는가에 따라 국소적 최소치에 해당하는 값이 다를 것이므로 가능하면 PES 상의 여러 점에서 출발하여 분자의 구조와 에너지를 조사할 필요가 있다. 이에 대해서는 발표 시간에 간단한 예를 들 것이다. 일반적으로, 조그마한 분자(약 60개의 원자로 구성된 분자)의 전체적 최소 에너지에 해당하는 구조를 계산으로 주장하기 위해서는 약 20000개에서 50000개의 다른 점들에서부터 에너지 최소화를 시도해서 그 중 가장 낮은 에너지 값을 가진 구조를 가장 안정한 분자 구조라고 주장하는 것이 통례였다. 그리고 최적화된 에너지 구조를 구하는 문제는 아직까지 매우 활발한 연구가 진행되고 있는 분야로 지금 사용되는 기법들도 그 종류가 매우 다양하다. MD를 이용하여 최적화된 구조를 찾는 경우에는, 온도를 충분히 높게 올렸다가 서서히 떨어뜨리면서 최적화된 구조를 찾는 방법이 가장 흔하게 사용되고 있다. 이는 온도가 높아지면 분자 구조를 결정하는 데 가장 큰 영향을 미치는 비틀림 각이 만드는 에너지 장벽(보통 상온에서 2~3RT 정도로 생각한다.)을 쉽게 넘나들 수 있기 때문에 단순한 MM보다는 좀 더 안정한 구조를 찾을 수 있는 가능성이 증가하게 되는 것이다. 이 외에도 임의의 여러 지점에서 출발을 해서 구조를 찾는 MC나 MD와 MC를 서로 혼합한 방법들도 많이 사용된다.

이 보고서에서는 이제 앞에서 이야기 한 세 과정 가운데서도 가장 핵심이 되는 부분인 경험적 힘장에 대해서만 간단히 살펴보기로 한다.

3.2 (경험적) 힘장 : (Empirical) Force Field : (E)FF

MM의 구조에서 가장 핵심적인 부분을 차지하는 것이 바로 FF이다. 일반적인 경우, PES를 얻기 위해서는 슈뢰딩거 방정식을 풀어야 하는데 반해, FF는 실험으로부터 얻은 자료 등을 활용해 PES를 재현하도록 만들어 낸 함수들의 집합이다. 따라서, 원래는 양자 현상인 화학 결합이나 분해를 FF를 이용하여 재현할 수 있는 이유는 우리가 다루는 계가 고전 역학적으로 행동하기 때문이 아니라 관찰되는 현상을 재현하기 위해 만들어낸 FF 속에 이미 평균화된 양자적 효과들이 고려되었기 때문이다. 따라서 광자의 흡수를 통한 전자의 전이나 전자 전달 현상, 혹은 양성자 전달 현상처럼 전적으로 양자적인 현상을 기술하는 일에는 FF를 기반으로 한 전통적인 MM을 활용할 수 없다.

이러한 FF의 이론적 연원은 진동 분광학자들이 진동 스펙트럼을 설명하기 위해 도입한 해석(analytical)식들에 있다. 이 식들은 대개의 경우 내부 좌표계를 사용하여 표시되어 왔다. 이는 데카르트 좌표계에서 분자계를 표현하는 것보다 식이 간단하게 표현되기 때문이지만, 또한 이것이 물리・화학적 입장에서 보았을 때 직감적으로 훨씬 자명하고 분자계가 갖고 있는 내재적 대칭성을 활용할 수도 있기 때문이기도 하다. 우선, 국소적이거나 전체적인 에너지 최소치에서 에너지 E가 E=E(Si), 즉 일반화된 좌표계인 {Si}의 함수로 표현된다고 하자. 그러면, 이 최소값 근처에서는 작은 정도의 변형에 대해 다음과 같은 테일러 전개식이 성립한다 :

E = E_o + SUM{G_i*(Del_S_i)} + SUM{K_i*(Del_S_i)^2} + SUM{F_ij*(Del_S_i)(Del_S_j)} + ...

이제, G_i = round dE / round dS_i 라고 하면 최소값에서 이 값들은 모두 0이 된다. 즉, E = E_o 가 되는데, 이 값을 분자에 내재하는 잔여’긴장(strain) 에너지라고 해석할 수 있다. 이제 최소점에서 점점 벗어날수록 두 번째 항은 PES에 대한 선형 근사가 된다. 세 번째와 네 번째 항은 모두 Del_S_i의 2차 전개 항이다. 그러나, 세 번째 항은 하나의 좌표계 S_i만을 포함하는 대각(diagonal)항들의 합이고 네 번째 항은 교차항(cross term)들의 합으로 이루어진 항이다. 여기에 조화 근사(harmonic approximation)를 사용하면 이것이 초창기의 FF에서 채택한 함수 꼴의 근간이 된다. 물론, 조화 근사 하에서 분자를 다룰 경우 Del_S_i가 커지게 되면 필연적으로 비 조화 항들까지 고려해 주어야 한다. 이렇게 표현된 (주로) 분자 내부 상호 작용에 의한 에너지 항에, 반 데르 발스 힘과 같은 원거리 작용 분산력과 전자기적 상호 작용 항과 같은 분자 간 상호 작용에 의한 에너지 항, 그리고 다루고자 하는 계에 특별히 존재하는 기타 항들을 함께 고려해 주어서 힘장(FF)을 완성하게 된다. 여기에서, 분자 간 상호 작용에 포함되는 에너지 항들은 그 꼴이 물리학적 이론에 의해 거의 해석적인 형태로 주어지게 되므로 수치 연산이나 계산 편의를 위한 식의 변형을 제외한다면 개념적으로는 명확하나, 이를 FF에서 구현하기 위해서는 사용하고자 하는 FF에 사용된 수식에 맞도록 적당하게 파라메터화할 필요가 있다. 그러나, 일상적으로 사용하는 분산력이 대개가 레너드 존스 포텐셜이나 모스 포텐셜이거나 또는 이를 약간 변형한 것이기는 하지만, 개개의 경우마다 파라메터를 개발해 주는 것 필요하다. 그러므로, 많은 경우 FF의 개발은 분자 내부 상호 작용을 정확히 묘사하는 FF의 수식 꼴의 개발과 계에 특별하게 존재하는 상호 작용을 고려하는 수식을 개발하는 일에 주된 발전 목표를 두고 있는 동시에 정확한 파라메터의 개발도 비슷한 비중으로 이루어지고 있다고 해야 할 것이다.

예를 들면, 용액 내 용질 분자간 상호 작용, 용매 분자와 용질 분자간의 상호 작용, 그리고 용매 분자간의 상호작용을 다룰 수 있는 함수꼴의 개발 등이 그것이다. 인터넷을 이용할 경우, 숭실대학교 분자설계연구센터의 교육 교재로 사용되고 있는 자료를 참조하는 것이 큰 도움이 될 것이라 생각된다. 특히, <MMㅡMD-MC-Vib-기초>를 참조하라. 이외에도, 양자 계산의 기초나 chemometrics 개론, 그리고 분자 모델링에 사용되는 좌표계 등에 관한 자료도 올라와 있다. 홈페이지 주소는 http://camd.soongsil.ac.kr/sites01.html 이다.

FF를 구성하는 공통적인 요소에는 원자 형식(atom type), 원자 형식에 포함되지 않은 전하, 원자 형식을 부여하는 규칙들, 에너지를 나타내는 함수 꼴들, 각 함수 항들에 대한 파라메터 등이 있다. 여기에서 말하는 원자 형식이란, 한 원자가 주변 환경에 따라 취할 수 있는 여러 가지 다른 혼성화, 전하 등이 각각의 FF가 목적하고 있는 바에 의거해 서로 구별 가능하도록 원자의 종류를 세분화시킨 원자의 모델이다. 예를 들면, 주기율표에는 탄소가 C 하나 뿐이지만, FF의 구성 목적에 따라 sp・sp2・sp3 혼성, 탄소 양이온, 탄소 음이온 등으로 세분하여 C1, C2, C3, C+, C- 등의 내부 기호로 구별하게 되는 다수의 C가 존재하게 된다. 그리고, 각각의 원자 형식에 따라 FF를 구성하는 함수 꼴에 포함되는 계수와 상수 부분이 일종의 파라메터 표로 저장이 되어 있다. 예를 들어, 탄소 분자간의 결합을 나타내는 함수 꼴이 후크의 법칙에 근거한 것이라면, 여기에는 탄소의 원자 형식에 따라 서로 다른 힘상수와 평형 거리가 결합에 필요한 파라메터 표로 저장되게 된다.

그러므로 이제 FF를 이해하기 위한 가장 핵심적인 단계는 S_i와 K_i, F_ij를 결정하는 일이다. 그런데, FF를 구성하는 함수들의 형식은 우리에게 매우 익숙한데, 왜냐하면 FF에서 채택한 내부 좌표계는 결합 길이, 결합 각, 비틀림 각(torsional angle), 그리고 나머지 비결합 자유도로 표현되기 때문이다. 그런데, 현재 사용하고 있는 FF들은 거의 대부분이 3N-6개의 자유도를 갖는 함수 꼴을 사용하지 않고, 더욱 많은 정보를 담기 위해 공유 결합 구조를 설명하는 일에 5N, 그리고 비 결합 상호 작용을 나타내는 일에 N^2/2 개의 자유도를 사용하고 있다. 따라서, 이들 함수 꼴을 구성하는 데 사용된 정보가 서로 모순되지 않도록 일관성을 유지하는 일이 무엇보다도 중요하며, 또한 서로 다른 FF 간에 함수 꼴이나 파라메터를 교환하는 일이 거의 불가능하다는 것이 자명하다. 그러므로, FF를 구성하는 함수의 형태는 그 FF의 응용 영역을 결정하는 중요한 요인이 되는 것이다. 그리고, FF를 구성하는 함수를 선택하는 일에는 컴퓨터의 계산에 적합한 형태의 함수를 선택하는 일이 중요하다. 일례로, 분자간 상호 작용을 다룰 때 모스 포텐셜이 좀 더 좋은 결과를 줌에도 불구하고 레너드 존스 포텐셜을 사용하는 이유는 후자가 계산하기가 훨씬 용이하기 때문이다. 이는 지수 함수의 계산에 시간이 걸리기도 하지만 흔히 사용되는 12-6 레너드 존스 포텐셜은 r^-6을 한 번 계산하면 r^-12는 제곱이라는 연산을 한 번만 취해주면 되기 때문에 알고리듬적으로도 매우 적합하다고 할 수 있는 것이다. 한가지 더 주의할 점은 가능하면 FF로 이루어진 PES가 연속일 수 있도록, 나아가 C1이 되도록 구성하는 것이 좋다. 왜냐하면 우리가 보고 싶어하는 성질이나 그것을 계산하기 위해 필요한 항 가운데 PES의 기울기 이상의 미분이 필요한 경우가 허다하기 때문이다. 그러한 이유에서, 우리는 주어진 함수를 통해 선험적으로 PES의 모양을 알아 낼 수는 없기 때문에 가능하면 다항식의 전개항 꼴로 FF를 구성하는 것이다.

그리고, 이렇게 FF를 각각의 내부 좌표계에 대한 전개식으로 나타냄으로서 얻는 중요한 이점 가운데 하나가 바로 ab initio나 기타 실험 자료를 이용하여 FF를 구축할 수 있다는 점이다. ab initio를 통해 직접 원하는 계의 PES를 구축한다는 것은 개념적으로도 결코 간단한 일이 아니다. 왜냐하면, 우리가 상식적으로 상정하게 되는 데카르트 좌표계(혹은 기타 직교 좌표계)를 기준으로 PES를 생각하게 되면 독립변수가 어떤 것들인지 확인하기도 힘들 정도이다. 그러니, 우선 분자 간 상호 작용과 분자 내부 상호 작용을 분리하여 생각하자. (물론, 나중에 분자 간 상호 작용과 분자 내부 상호 작용의 결합(coupled)항들을 넣어 줌으로서 이 분리로 인한 오차를 보정할 수 있다.) 그러면, 우리는 결합 길이, 결합 각, 비틀림 각 등을 한 종류의 독립 변수로 보고 이 독립 변수 한 종류에 대해서만 전개식에 나타나는 계수나 상수 파라메터를 ab initio 계산 등을 통해 구하게 된다. 즉, 이렇게 결합 길이에 대한 식과 관련 파라메터를 따로 구축하고, 그리고 결합 각, 비틀림 각 등에 대해서도 같은 일을 해주면 대강의 PES를 만들 수가 있는 것이다. 물론, 결합 각의 변화에 따라 결합 길이도 변하는 등의 효과를 고려하기 위해 결합항들을 고려해 주면 좀 더 세련된 묘사가 가능하다. 이 경우, 주의해야 할 것은 ab initio 계산 등의 결과를 이용해 FF를 구성하는 각각의 함수 꼴에 대한 파라메터를 조정해 나갈 때, 전반적으로 일치하는 모양을 얻어야 하는 것이 아니라 최소값 부근의 값들을 가장 정확하게 맞추는 일을 우선 해야한다. 왜냐하면, 에너지 최소화 과정에서 가장 많은 정보를 이용하게 되는 부분이 바로 에너지 최소값 근처 부근이기 때문이다.

3.3 많이 사용되는 (E)FF의 예

가장 널리 사용되는 FF들에 대하여 간략하게 정리해 보면 다음과 같다.

이는 Goodman, 『Chemical Applications of Molecular Modelling』, pp. 190~192에 나온 부분을 중심으로 하여 Gelin, 『Molecular modeling of Polymer Structures and Properties』, pp. 32~33, 및 『Forcefield-Based Simulations : General Theory & Methodology』, pp. 18~60에 나온 내용에 몇몇 기본적인 참고 사항을 첨가한 것이다. 특히, 마지막 책에서는 개별 힘장들이 취한 각각의 함수꼴들을 간단명료하게 정리해 놓고 있어, 힘장마다 가지고 있는 독특한 구조를 파악하는데 많은 도움을 얻을 수 있다.

AMBER(Assisted Model Building with Energy Refinement)  : 콜만(Kollman)이 펩타이드와 핵산의 모의실험을 위해 개발한 FF로, 주로 그 분야에 널리 사용되고 있으나 다른 경우에 있어서도 널리 사용되는 범용 FF 가운데 하나이기도 하다. 단, 재료 과학 분야의 연구에는 사용하지 않도록 권장되고 있다. AMBER와 관련된 주요 논문은 P. K. Weiner & P. A. Kollman, Journal of Computational Chemistry, 2, pp. 287~303, 1981과 S. J. Weiner, P. A. Kollman, et al, Journal of Computational Chemistry, 7, pp. 230~252, 1986, 그리고 S. J. Weiner, P. A. Kollman 등이 쓴 JACS, 106, pp. 765~784, 1984가 있다. 데카르트 좌표계를 사용하며, 함수는 주로 2차 항까지 전개된 함수 형태를 사용한다. 기하 구조 데이터와 진동 스펙트럼 데이터 모두를 이용하여 FF를 만들었다. diagonal한 함수 꼴 외에 torsion과 비결합 함수 형태도 사용하며 다양한 원자 형식을 보유하고 있어 원래는 all atom FF이지만 수소를 명확히 나타낼 수도 있고 그렇게 하지 않을 수도 있다. 홈페이지는 http://www.amber.ucsf.edu/amber/이다.

CHARMm(Chemistry at Harvard Macromolecular Mechanics)  : FF의 이름으로보다 단백질 및 생체 분자를 위한 범용 MM/MD(/MC) 프로그램으로 더 유명하다. 하버드의 카플러스 그룹이 만든 포트란 코드로 고분자량의 단백질에 관련된 계산에 최적화되어 있다. 주요 논문은 B. R. Brooks, et al and M. Karplus, Journal of Computational Chemistry, 4, pp. 187~217, 1983과 L. Nilsson & M. Karplus, Journal of Computational Chemistry, 7, pp. 591~616, 1986이 있다. 채택한 좌표계나 함수 형태는 AMBER와 비슷하나, 에너지 최소값들을 고려함에 있어 좀 더 다양한 경우의 데이터를 활용해 FF를 만들었다. 최근에 ab initio 계산을 이용하여 CHARMm FF를 개선하려는 움직임이 있으며, C(++) 코드로의 포팅 프로젝트도 진행 중이다. 한양대 화학과의 원영도 교수님이 카플러스 그룹에서 박사 후 연구원 과정에 계시면서 CHARMm 개발에 기여하였다.

CFF93  : 1979년에 리프슨 등에 의해 개발된 CFF를 모태로 한 FF로, 해글러(Hagler) 등에 의해 개발되었다. CFF는 H, C, N, O, S, P와 할로겐 원자나 이온, 알칼리 금속 양이온, 그리고 생화학에서 중요하게 사용되는 2가 금속 양이온 등을 포함한 광범위한 유기 화합물의 실험 자료로부터 구축한 힘장이다. 또한, 조그만 분자와 커다란 분자 모두를 정확하게 정의하기 위해 비조화 항들과 coupling 항들을 이용하여, 개발자들 스스로 CFF93을 기존의 FF와는 정확도에서 차이를 보인다는 함의로 Class II라고 분류하였다. CFF93은 매우 정밀한 대신에 함수 꼴이 무척 복잡하고 엄청난 양의 파라메터가 사용되었다. 그리고 스케일링을 위해 특별히 7개의 파라메터를 도입하였다. 이 파라메터들의 이식성에 관해서는 M. Hwang, et al, JACS, 116, pp. 2515~2525, 1994를 참조하라. CFF 힘장은 이후 PCFF, CFF95 등의 개발로 이어진다. 그리고, CFF93 이후의 힘장은 AMBER나 기타 고전적인 힘장들보다 좋은 결과를 준다는 것이 알려져 있다.

DREIDING  : 매이요 등에 의해 개발된 FF로, 유기・생 분자 및 다수 무기 분자들의 구조 예측과 동역학을 관찰하는 데 사용할 수 있는 범용 FF이다. 채택한 좌표계나 함수의 형식은 AMBER나 과 유사하나 원자 형식이 최소화되어 있으며, 기하구조와 torsion barrier 데이터를 이용해 FF를 구축하였다. 이 FF의 목적은 글자 그대로의 ‘범용’ FF 구현에 있고 지금까지의 발전 방향도 그 기본 철학에 충실하고 있다. 금속 화합물을 포함하는 화합물에 대한 계산 역시 가능하도록 설계되었다. 새로운 화합물이나 실험 데이터가 부족한 분자 계에 대해 좋은 기준이 되는 결과를 제공한다고 알려져 있다. 유기, 생 그리고 main-group 무기 분자들의 구조 예측과 동역학을 모의 실험하는 데 많이 사용된다. 주요 논문으로는, S. L. Mayo, et al & W. A. Goddard III, JPC, 94, pp. 8897~8909, 1990이 있다.

MM2/MM3/MM4  : MM2는 1977년에 앨링거(Allinger)가 개발한 FF로, 원래는 조그만 분자들을 주된 대상으로 개발하였으나 상당히 광범위한 분자들에 대해서도 좋은 결과를 주었다. 그리하여, 새로운 FF는 MM2와 비교해서 어떤 점에서 더 좋고 어디가 문제인지를 지적해야 하는 일종의 표준으로 자리 매김을 했다. MM2와 관련한 주요 논문은 N. L. Allinger, JACS, 99, pp. 8127~8134, 1977가 있다. MM3는 MM2에 이어 앨링거가 1989년에 개발한 FF로 당연히 MM2의 몇몇 결점을 보완하는 데 그 목적을 둔 FF이다. 특히, MM2는 진동 주파수를 잘 재현하지 못했었는데 MM3에 이르러 이 문제를 어느 수준까지는 해결했다. MM3에 관련된 주요 논문으로는 L. Allinger, et al, JACS, 111, pp. 8551~8556, 1989가 있다. MM 시리즈는 cross term들을 명확한 함수 형태로 사용하고 기하 구조와 진동 스펙트럼 외에도 생성열을 이용해 FF를 구축해 왔기 때문에 크지 않은 분자에 대해서는 가장 정확한 FF로 군림해 왔다. 그리고 1996년에 이르러 MM4를 내놓으며 이제까지 정확하게 계산하기 힘들었던 진동 주파수나 회전 장벽 에너지도 만족스럽게 재현할 수 있게 되었다. MM4에 관련된 주요 논문으로는 L. Allinger, et al, Journal of Computational Chemistry, 17, pp. 641~668, 1996이 있다.

MMFF(Merck Molecular Force Field)  : 머크의 계산 화학자인 할그렌(Halgren) 등에 의해 개발된 FF로, 작은 분자들뿐만 아니라 단백질에 관한 계산에서도 상당히 정확한 결과를 주는 것으로 알려져 있다. FF 제작에 필요한 데이터 중 많은 부분을 상당히 높은 수준의 기저 함수를 이용한 ab initio 방법에 의존한 것으로 알려져 있다. 이 FF는 MD에도 상당히 적합하다. 주요 논문으로는 T. A. Halgren, Journal of Computational Chemistry, 17, pp. 490~519, 1996이 있다. MMFF93의 주된 용도는 단백질이나 핵산과 같은 거대한 수용체(receptor)와 다양한 종류의 리간드 간에 존재하는 상호작용을 연구하는 데 있다.

UFF(Universal Force Field)  : 1992년에 라페(Rappé) 등에 의해 개발된 FF로, 다른 FF에서는 기능기에 근거한 원자 형식을 사용한 것과 대조적으로 오로지 원소 하나 하나와 혼성화, 그리고 결합 양식에만 근거해 이를 파라메터화한 FF이다. 이 FF의 목표는 주기율표 전체에 걸친 모든 원자들을 간단하게나마 모두 다룰 수 있는 FF를 만드는 것이다. 당연히 금속 화합물에 대해서도 계산이 가능하다. 발상은 매우 신선하나, 각각의 상황에서 이 FF가 얼마나 신뢰할 만한가를 어떻게 판단할 것인가에 대해서는 아직도 일치된 바가 없다. 하지만, main-group 화합물과 유기 화합물, 그리고 다양한 금속 착화합물에 대해 어느 정도 만족스러운 결과를 줌이 보고된 바 있다.

『Forcefield-Based Simulations : General Theory & Methodology』, pp. 18~60. UFF에 관련된 주요 논문으로는 A. K. Rappé, et al, JACS, 114, pp. 10024~10035, 1992가 있다.

Reax FF(First Principles Force Fields)  : 칼텍의 Goddard III group에서 만든 FF로, 체계적인 함수꼴을 이론적으로 제안하고, 그 구현에 필요한 값들을 양자 계산(주로 DFT-GGA)을 통해 얻어 새롭게 만든 힘장이다. 주된 이용 목적은 Mesoscale modeling을 하기 위함이나, 그 정확도와 결과의 만족도는 왠만한 양자 계산을 무색케 한다해도 과언이 아니다. 10년 가까이 꾸준히 계발 보완 발전되고 있으며, 각 영역별로 그 성능에 관한 발표가 보고되고 있어 그 신뢰도는 날이 갈수록 높아지리라 사료된다.

이 외에도 주기 경계 조간 하에서 무기 산화 유리에 대한 계산을 위해 개발된 Glassff라든지, 유기금속화합물을 다룰 수 있게 설계 된 SHAPES 등 특수한 목적을 위해 만들어진 수많은 FF들이 존재한다. 이렇게 다양한 FF에 관해 어느 정도의 정보를 가지고 있을 떼에야 우리는 좀 더 광범위한 계를 보다 정확하게 다룰 수 있을 뿐만 아니라, 계산 결과 가운데 개별 힘장에 고유한 인위적 기여도가 어느 정도인지 판단할 수도 있고, 나아가 새로운 힘장을 개발하는 일도 시도할 수 있다. 그리고 무엇보다도, 아직까지는 어느 상황에서나 좋은 결과를 주는 하나의 FF는 존재하지 않으므로 MM에서는 다루고자 하는 계에 가장 최적화된 FF를 선택하는 일이 무엇보다도 중요하고 이를 위해서는 반드시 다양한 힘장들에 대한 기본적인 지식을 가지고 있을 필요가 있다.

4 분자 동역학(Molecular Dynamics)

이 부분은 Frenkel & Smit, 『Understanding Molecular Simulation』의 4장을 중심으로 몇몇 다른 자료들을 참고한 것이다.

컴퓨터를 이용한 모의 실험을 통해 우리는 다체계의 특성을 연구할 수 있다. 그러나, 컴퓨터 모의 실험을 통해 원하는 모든 성질을 직접적으로 측정할 수 있는 것은 아니다. 또, 모의 실험을 통해 얻은 양이 실제 실험에서 측정한 성질에 대응하는 것도 아니다. 즉, 모의 실험 과정을 통해 얻는 데이터들과 우리가 실험적으로 얻는 계의 특성들이 일대일 대응 관계를 보이지 않는다. 왜냐하면 우리가 실제 실험을 통해 얻는 자료들은 원칙적으로 특정한 평균값이기 때문이다. 그러므로, 모의 실험을 할 때에는 어떤 종류의 평균들을 계산해야 할지 미리 판단하여야 한다. 이 지점에서 통계 역학이 개입한다. 컴퓨터 모의 실험과 실제 실험이 어떠한 관계에 있는 지는 우리가 모의 실험 결과로부터 어떤 평균값을 취해야 할 것인가에 대해 판단을 내린 후에야 가능한 것이다. 그렇다면, 이러한 결정이 내려진 후 우리는 과연 무엇을 할 것인가? MD를 통해 물리적 가관측량을 측정하기 위해서는 우선 원하는 양을 계를 구성하는 입자의 위치와 속도의 함수로 표현할 수 있어야 한다. (이는 우리가 다루는 물리학의 체계에서 비롯된 것이므로 다른 해석 체계를 가지고 있다면 다른 변수를 이용할 수도 있겠다.) 그리고, 계를 구성하는 입자의 위치와 속도 자료가 시간에 따라 축적되면 이제부터는 우리가 배운 고전 혹은 양자 통계역학이라는 해석 도구를 이용해 원하는 물리량을 얻게 된다. 이 때, 수업 시간에는 흔히 간과하고 넘어가던 자료의 오차와 관련된 양(예를 들면, fluctuation)이 매우 중요한 정보원이 된다. 왜냐하면, 우리가 실제로 계산을 통해 얻은 궤적에 관련된 자료들이 바로 들쭉날쭉한 숫자들의 집합이기 때문이다. 이제 MD의 내부 구조를 간략하게 살펴봄으로써 이를 이해해보고자 한다.

고전적 의미에서의 MD란 고전적 다체계의 평형 성질과 운반 성질을 계산하는 기술로, 여기에서 고전적이라 함은 핵의 운동을 고전 역학에 따라 기술한다는 의미이다. 하지만, MD의 용도는 사용하는 사람들에 따라 확장될 수 있을 것인데, 지금까지는 주로 형태 분석(conformation search), 통계 역학적 앙상블의 생성과 그를 이용한 계의 분석, 그리고 분자 운동의 연구 등에 사용되고 있다. 여기에서, MD를 통한 구조 분석에서는 MM을 통한 구조 분석에서보다 훨씬 많은 구조에 관해 데이터를 모을 수 있다는 장점 외에도 우리가 실제 관찰하게 되는 시간에 대해 평균된 결합 길이 등의 정보를 직접적으로 얻을 수 있다는 장점이 있다. 아무튼, 이러한 MD는 컴퓨터 모의 실험이기는 하지만 개념적으로는 많은 부분에서 실제 실험과 매우 유사하다.

① 샘플을 준비하는 과정 : N개 입자로 구성된 모델 시스템을 선택한다.
② 평형을 기다리는 과정 : 뉴턴 방정식을 계의 성질이 시간에 대해서 일정할 때까지 푼다.
③ 측정을 실행하는 과정 : 평형에 도달한 계의 성질을 실제로 측정한다.

이제, 실제 MD 프로그램의 내부를 간단하게 살펴보자.

① 필요한 인수를 읽는다.
② 계를 초기화한다 : 계의 초기 위치와 속도를 설정한다.
③ 모든 입자 각각에 작용하는 힘을 모두 계산한다.
④ 뉴턴 방정식을 적분한다 : 원하는 시간까지 계산한다.
⑤ 측정을 원하는 양의 평균값을 얻는다.

e. g. generic algorithm

       program md        call init                               ! initialization        t=0        do while (t .lt. tmax)          ! MD loop
               call force(f, en)               ! determine the forces                call integrate(f, en)           ! integrate the equation of motion                t=t+dt                  call sample             ! sample averages
       enddo           stop                    end

4.1 초기화

모의실험을 시작하려면 우선 각각의 입자에 초기 위치와 속도를 배정해 주어야 한다. 이 때 입자들끼리 서로 겹치거나 해서는 안 된다. 그럴 경우, 계의 에너지가 물리적으로 무의미하게 커지기 때문이다. 운동 에너지가 존재하므로, 계의 온도 역시 저의 가 가능하고 따라서 계의 초기 온도 역시 설정해 주어야 해야 한다. 모의실험 상에서의 온도는 통계 역학의 결과인 등분배 법칙(equi-partition law)을 이용해 k_B*T(t) = SUM_i{mv_a,i(t)^2/N_f}와 같이 정의한다. 따라서 순간 온도 T(t)를 원하는 온도인와 맞추기 위해서는 모든 속도를 (T/T(t)^(1/2)로 스케일링 해주면 된다. 그러면, 볼츠만 분포나 가우시안 분포를 따르는 임의의 수를 발생시켜 각각의 분자들에 적당한 속도를 부여하게 된다. 하지만, 이것은 온도의 초기화에 사용되는 기법일 뿐, 정준 앙상블 모의 실험을 하기 위해 온도를 일정하게 유지하는 방법으로 사용될 수는 없다.

4.2 힘의 계산과 운동 방정식의 적분

많은 경우 알고리듬 상에서 속도 그 자체를 이용해서 뉴턴 방정식을 풀지는 않는다. 힘과, 현재와 전 단계에서의 입자들 위치를 이용해 다음 단계에서의 위치를 구하는 것이 일반적이다. 그런데, 대개의 MD에서 가장 많은 시간이 소요되는 단계가 바로 이 힘을 계산하는 단계이다. 입자 수(N)의 승으로 시간이 필요한 알고리듬도 많이 사용되지만, 논의의 명확성을 위해 가장 간단히 N의 제곱 승으로 시간이 필요한 알고리듬을 생각해 보자.

① 입자 i와 J간의 x, y, z 성분에 대해 각각의 현재 거리를 구한다. 물론, 우선은 피타고라스 정리에 따른 거리 제곱 값을 구하는 것이다. 왜냐하면, 제곱근을 계산하는 것이 알고리듬에 따라 불필요한 경우도 많거니와 그 계산에 무척이나 많은 시간이 소요되기 때문이다. 그리고, 대개의 경우에는 단순 정육면체 주기 조건에서 i와 j의 최 근접 주기 상과의 상호 작용을 생각하게 된다.

② 이제, 충분히 가까운 거리 내에 있는 입자들간에 작용하는 힘을 계산한다. 이제 우리는 모든 힘을 다 알고 있으므로 뉴턴의 운동 방정식을 적분할 수 있다. 예전에 가장 널리 사용되던 알고리듬이 바로 베를레이 알고리듬이다. 이 알고리듬은 간단하면서도 대개의 경우 아주 좋은 결과를 준다. 고전적인 베를레이 알고리듬에서는 새로운 위치를 계산하는 일에 속도를 사용하지 않는다.

정확한 속도의 계산은 무척이나 중요한데, 왜냐하면 이것이 운동 에너지와 직결되어 있기 때문이다. 또한, 운동 에너지는 순간 온도와 직결되어 있다. 그리하여 베를레이 알고리듬과 동일한 궤적을 주는 베를레이-유사 알고리듬들을 사용한다.

속도가 필요할 경우에는 기초적인 유한 차분법을 이용해 계산할 수가 있다. 이 경우, 새로운 위치는 Del_t의 4승의 오차를 갖지만, 속도는 Del_t의 2승의 오차를 갖게 된다.

MD에서의 운동 방정식을 적분하는 알고리듬들은 사실 수치해석에서 사용해온 ODE solver들 내지는 거기에 바탕을 둔 알고리듬이기 때문에, 구현하는데 있어 별다른 어려운 점은 없다. 하지만, 매우 오랜 회수를 반복하며 적분을 해야하기 때문에 MD에서 실제 사용할 integrator를 선택하는 것에는 몇 가지 기준이 존재한다. 우선 드는 생각은 적분하는데 걸리는 시간이 적어야 할 것 같은데, 사실 운동 방정식을 적분하는 일에는 그다지 많은 시간이 요구되지 않기 때문에 시간이라는 요소를 희생하더라도 다음의 요건을 만족하는 알고리듬을 사용하는 것이 매우 중요하다.

① 큰 시간 간격을 사용해도 정확한 새 위치를 구할 수 있어야 한다 : 이는 고차 미분 계수를 저장해 놓고 그 값을 사용하는 것으로 구현 가능하다. 이것은 거의 당연한 사실이지만, 계산 시간과 정확도 사이에서 어떤 절충 지점을 찾을 수 있을 지는 다루고자 하는 계의 성질과 보고자 하는 성질을 명확하게 파악하고 있을 때에야 가능할 것이다.

② 에너지 보존 법칙을 잘 재현할 수 있어야 한다 : 에너지 보존 법칙은 단시간 에너지 보존 법칙과 장시간 에너지 보존 법칙을 나누어 생각할 필요가 있다. 고차 미분계수를 사용하는 세련된 알고리듬의 경우에는 짧은 시간의 모의 실험에서는 아주 훌륭한 에너지 보존 법칙을 재현하지만, 장시간의 걸친 모의 실험에서는 전체 에너지를 잘 보존하지 못하는 경향이 있다. 이와는 대조적으로, 베를레이 류의 알고리듬은 짧은 시간에 걸친 모의 실험에서는 에너지 보존 법칙을 재현해 내는 정도가 썩 좋다고는 할 수 없으나, 장시간에 걸친 모의실험에 대해서는 에너지 보존 법칙을 잘 재현해 낸다. 그러나 아직까지는 이 두 종류의 에너지 보존 법칙을 모두 만족하는 알고리듬이 존재하지 않는다. 따라서 리아프노프 불안정성과 에르고딕 성질을 잘 고려해 계에 적합한 적분 알고리듬을 택할 필요가 있다.

③ 뉴톤 방정식이 그런 것처럼 알고리듬도 시간 가역적이어야 한다 : 하지만, 알고리듬에서 미래와 현재의 위상 공간 좌표가 대칭적인 역할을 하지 않는 한 이것을 컴퓨터 모의 실험에서는 거의 재현할 수 없다. 하지만, 알고리듬 선택에서 시간 가역성보다 더 중요한 것은 알고리듬이 심플렉틱한가의 여부이다.

대개의 경우는, 그라스만 대수의 꺾쇠곱(wedge product)을 보존하는가의 여부에 관한 것으로 해석해도 된다.

대개의 시간 비가역적 알고리듬은 심플렉틱하지 않다. 심플렉틱 하지 않은 알고리듬은 위상 공간에서 우리 시스템의 부피가 점점 팽창하게 만들 것이고 이는 에너지 보존 법칙에 위배된다! 그래서, 많은 시간 비가역적 알고리듬에서 심각한 장기간 에너지 보존 법칙 위배가 발생하게 되는 것이다. 현재와 미래의 위상 공간 좌표 변환과 관련된 야코비안을 통해 우리는 알고리듬이 위상 공간에서 부피를 보존하는가 아닌가를 확인해 볼 수 있다. 가장 널리 사용되어 온 베를레이 알고리듬을 이 기준에 비춰 평가한다면 다음과 같은 얘기를 할 수 있다.

① 빠르다 : 하지만 컴퓨터의 성능이 날로 향상되고 있는 요즈음에는 이것이 최우선적으로 고려되어야 할 점이라고 하기는 어렵다.

② 시간 간격이 클 경우에는 정확하지 못하다 : 따라서, 우리는 힘을 자주 계산해야 한다.

③ 적은 메모리를 필요로 한다 : 사실, 이것도 그 중요성이 많이 퇴색하고 있다.

④ 단기간 에너지 보존도 그럭저럭 잘 재현해내는 데다가, 장기간 에너지 보존 법칙을 매우 훌륭하게 재현한다고 볼 수 있다 : 베를레이 알고리듬은 시간 가역적이고 위상공간에서 부피를 보존하는 것으로 알려져 있다. 또한, 베를레이 알고리듬은 계의 총 에너지를 정확하게 보존하지는 못하지만, 무한소 시간에서 실제 헤밀토니안으로 수렴하는 유사 헤밀토니안을 보존함이 알려져 있다. 베를레이 알고리듬을 통해 구한 궤적 자체는 그리 훌륭한 것이 못되지만, 이러한 장점 때문에 우리는 베를레이 알고리듬을 오랫동안 이용해 온 것이다.

대부분의 MD에서는 베를레이 유사 알고리듬(특히, 속도-베를레이 알고리듬)이면 대부분 만족할 만한 결과를 얻을 수 있다. 하지만, 단기간 정확도를 유지하면서 큰 시간 간격을 사용하려고 하거나 혹은 주어진 시간 간격에 대해 좀 더 정확한 결과를 얻고 싶을 때에는 예측자-수정자 알고리듬과 같은 고차 알고리듬을 사용하기도 한다. 그러나, 앞에서도 말했듯이 고차 미분 계수를 사용하게 되면 시간 가역적이지도, 부피 보존적이지도 않게 된다. 따라서 원하는 성질을 확실하게 설정하고 모의 실험을 해야 목적에 부합하는 결과를 얻을 수 있는 것이다.

이제 많이 사용되고 있는 몇몇 다른 알고리듬을 생각해 보자.

① leapfrog 알고리듬 : 반 정수 시간 간격에서의 속도를 구하고 이를 이용해 새 위치를 구하는 알고리듬으로 베를레이 알고리듬을 약간 변형하여 쉽게 얻을 수 있다. leapfrog 알고리듬은 베를레이 알고리듬에서 바로 유도된 것이므로 베를레이 알고리듬과 동일한 궤적을 준다. 그러나, 여기서 사용되는 속도와 위치는 알고리듬 상에서 정의되는 시간이 서로 다르다. 이는 결국 운동 에너지와 위치 에너지가 동시에 정의될 수 없음을 의미하고, 따라서 leapfrog 알고리듬을 따를 경우 우리는 전체 에너지를 직접 구할 수가 없다. 물론, 적당하게 수식을 변형하여 운동 에너지나 포텐셜 에너지 가운데 한 쪽은 평균값을 사용하는 등의 방법도 있으나, 그럴 경우 자료의 요동(fluctuation)에 숨겨진 많은 정보를 모두 상실하게 된다.

② 속도-베를레이 알고리듬 : 동시에 정의된 속도와 위치를 이용하도록 베를레이 알고리듬을 다듬어 정리하면, 새로운 위치를 구한 후에 새로운 속도를 구하고 이를 이용해 새로운 힘을 구하는 속도-베를레이 알고리듬을 얻을 수 있다. 이 알고리듬이 장시간의 MD에는 가장 널리 사용되는 알고리듬이다.

③ 베만 알고리듬 : 원래의 베를레이 알고리듬에서보다 정확한 속도를 구할 수 있다. 따라서 좀 더 훌륭하게 전체 에너지가 보존되는 것처럼 보인다. 하지만, 베만 알고리듬은 시간 가역적이지 못하다.

④ 속도 수정 베를레이 알고리듬 : 이 알고리듬은 위치와 속도 모두 오차가 Del_t의 4승으로 주어진다. 그러나 이 알고리듬에서는 현재 속도가 다음 시간 단계에서 계산되어야 하므로 v(t)를 알기 위해서 t + Del_t에서의 위치와 힘을 미리 알아야만 한다.

이 외에도 거의 모든 상미분 방정식에 응용이 가능한 룽게-쿠타 4승 알고리듬이 사용되기도 한다: 룽게-쿠타 알고리듬은 매 단계마다 에너지를 구하여야 하고 시간 간격도 매우 작아야 한다는 제약이 있어서 긴 시간에 걸친 동역학을 보는 일에는 잘 사용되지 않는다. 하지만, PES 상에서 움직이는 분자들의 반응동역학 등을 연구하기 위한 분자동역학 모의실험에서는 정확도가 아주 중요하기 때문에 대부분 룽게-쿠타 알고리듬이나 그에 기반한 알고리듬을 사용하고 있다. (한동안 김명수 교수님 연구실에서 publish 했던 논문들 가운데 이런 예가 상당히 많이 보인다.)

4.3 고전 역학의 한계를 넘어서

앞서 이야기하였듯이, 고전적인 분자 동역학은 입자(원자, 분자)들이 적당한 다차원 힘장에 의해 고전 역학의 운동 법칙을 따라 움직인다고 가정하고 있다. 그러므로, 이 방법의 한계점은 다소간 명확하다 할 수 있다. 우선, 정확한 힘장을 구축하는 것이 매우 힘들며, 많은 수의 입자에 대해서 오랜 시간에 걸쳐 적분을 해야한다는 점도 당연히 문제가 된다. 그리고 정확한 통계 역학적 앙상블을 얻는 일 역시 힘들다. 하지만, 무엇보다도 양자 역학적 움직임이 모두 무시되고 하나의 PES가 입자의 운동을 지배하게 되므로, 수소 원자의 움직임이 중요시되는 저에너지 과정에서의 터널링이나 간섭 현상 등을 기술하기란 거의 불가능했다.

이에 대해서는 Tully, Molecular dynamics with electronic transition, Journal of Chemical Physics, vol. 93, pp. 1061~1071을 참조하라.

그리하여 90년대 들어와서는 이러한 고전적 의미의 분자 동역학의 한계를 극복하기 위한 노력들이 본격적으로 가시화되기 시작했다. 그 가운데 가장 초기에 이루어진 발전을 가운데 하나는 전자 전이를 고려해 줌으로서 단일한 PES에서만 움직일 수 있어서 야기되는 한계점을 극복하려는 노력이다. 또한, 카와 파리넬로가 1985년에 제시한 CPMD는 경험적 힘장이 아니라 양자 역학적 효과를 고려한 독특한 라그랑지안에 기반하여, 전자구조 계산 방법을 사용하여 구한 포텐셜 에너지로부터 분자들의 움직임을 계산해 나가는 방법으로, 최근 컴퓨터의 성능이 향상되고 병렬화가 쉬워짐에 따라 매우 각광을 받고 있는 방법이다. 또한 파인만 경로 적분을 이용해 양자 효과가 고려된 분배함수를 구하는 PIMD 방법도 최근들어 그 사용례가 늘어나고 있다. PIMD가 양자 효과를 고려한 계의 정적 성질을 관찰하는 데 그 목적을 두고 있다면, 파인만 경로 적분에서 나타나는 센트로이드(centroid)를 이용해 양자 효과를 고려한 동적인 계의 성질을 탐구하는 방법(특히, Voth Group의 Cao가 개발한 CMD)도 개발되었다. 즉, 이와 같이 분자 동역학이 갖는 매력과 그 한계를 극복하려는 노력은 바로 이 시점에서도 현재 진행형인 것이다.

5 맺으며

얼마 전까지만 해도 병렬 프로그래밍이나 PC 클러스터링은 전산 관련학과 사람들의 전유물이었다. 하지만, 이제는 상업용이 아닌 분자 모델링 프로그램들 가운데에서도 다수가 병렬 옵션을 지원하며 나아가 권장하고 있기까지 하다. 이는 가격 대 성능 비라는 관점에서 앞으로도 불가피하게 지속될 경향이라 하겠다. 그러므로 우리가 컴퓨터 모의 실험을 이용해 다룰 수 있는 계의 크기는 더욱 커질 것이고 그 정확도도 지금보다 훨씬 높아질 것이다. 그리하여, 지금은 거의 대부분의 ab initio 계산을 위한 프로그램들이 ‘단추만 누르면 되는’ 팩키지의 형태로 보급된 것과 마찬가지로, 얼마 지나지 않아 다양한 목적의 유사한 프로그램들이 팩키지화되어 춘추전국시대를 방불케 하는 시장 쟁탈전을 벌일 것임에 틀림없다. 그렇지만, 호프만이 이야기한 것처럼 프로그램이 주는 것은 방대한 숫자들로 이루어진 데이터일 뿐 그로부터 유의미한 정보를 얻어내는 일은, 적어도 당분간은 사람의 몫이다. 그리고, 택할 수 있는 방법론이 다양한 만큼 이제는 다루고자 하는 계에 가장 적합한 방법론을 결정하는 일이 많은 부분을 차지하게 될 것이다. 그러기 위해서는, 계산 분야에서 자주 사용되는 단어들을 어느 정도는 익숙해질 필요가 있다 하겠다.

물론, 이 보고서는 이제는 누구나 사용하는 일상적인 단어 몇 개를 나열하는 수준에서 그쳤지만 이는 보고자의 능력 한계와 보고서의 성격에 기인한 것이지 분야의 협소함 때문은 아니다. 또한, MM에서 기본적인 설명 없이 흔히 등장하는 단어들만을 중심적으로 다루다보니, 주기 경계 조건하에서 쿨롱 힘에 의한 상호 작용을 고려하는 방식 가운데 하나인 Ewald summation이나, MD에서 소정준 앙상블 외의 앙상블을 구현하는 방법 등 거의 모든 계산 화학 책에서 기본적으로 다루고 있는 자세한 내용들은 모두 생략하였다. 이러한 것은 소위 말하는 기본서들 거의 대부분이 다루고 있는 내용이기 때문에, 반드시 소개를 할 필요가 있으나 보고서의 성격상 과감히 생략하였다. 대신, 참고한 도서는 거의 모두 각주로 처리하였으므로, Numerical Recipes 등의 수치 해석 관련 서적을 제외한 책들 가운데 가장 많이 언급되는 초급 ~ 중급 참고서적들을 간단히 소개하며 보고서를 끝내려 한다. 얼마 되지 않은 기간이지만, 개인적으로는 좋은 책을 추천 받는 일이 얼마나 중요한 가를 절감하였으므로 한 사람만이라도 계산화학참고문헌 목록 가운데 도움을 받을 책을 건질 수 있다면 그것만으로도 충분히 만족스러운 일이라 생각한다.

6 p.s. NAMD / VMD

지금은 University of Pennsylvania의 Saven Group에 계시는 이원선 박사님께서 연재하시고 계신 NAMD 기본 사용법을 NAMDandVMD에서 읽어 볼 수도 있다.

7 같이 보기

• 분자역학
• 분자동역학
• 결합 쌍극자 모멘트