[C.Chem] Gaussian 16을 이용한 QM 계산의 기본 절차

Chemical reaction의 reactivity를 계산하는 procedure

지난 학기 대학원 수업에서 계산법 소개용으로 작성했던 프레젠테이션인데, 그냥 버리기는 아까워 여기에 적어볼까 한다. 사실 계산화학을 하는 사람이 대한민국에 얼마나 있는지 잘 모르기 때문에 큰 효용성이 있을지는 의문이지만..

실험실마다 방법이 다르겠지만, 내가 속한 실험실에서의 계산 진행의 절차는 대략 다음과 같다.

1. Draw overview first.

항상 모든 계산의 시작은 '어디로부터' 시작해서 '무엇을' 위해 하는지를 생각하는 것이다. 하나의 방향성을 가지고 있지 않으면 계산 중간중간에 나타날 수많은 난관들을 극복하기가 매우 까다롭다.

간단한 overview를 논문을 읽거나 랩미팅, 저널클럽 등에서 스터디한 후 'Chemdraw'나 'Chemsketch'를 이용하여 반응의 전반적인 mechanism을 러프하게 그리면 된다.

 

2. Take the structure & Structure optimization

실험하는 연구실과 co-work를 진행하는 경우, 해당 연구실에서 crystal structure나 NMR 자료를 가지고 있는 경우가 있다. 이러한 경우에는 일이 쉽게 풀릴 수 있는데, PDB와 같은 파일 포맷을 직접 가져와서 계산을 진행할 수 있기 때문이다. 물론 실험실에서 가지고 있는 물질의 부피가 많이 bulky 하다면 일부를 트리밍해서 계산에 적합한 부분만 남겨놓을 필요가 있다. 만약 가져올 구조가 없다면 프로그램을 이용하여 초기 구조를 직접 그려야 한다. 

다음으로 가져온 구조로부터 최적화된 구조를 찾는다. 이 과정에서 사용하는 키워드는 'opt'이다. opt는 'optimization'의 약자로 구조를 최적화하라는 명령어이다. 연관된 키워드에는 'rcfc' 'readfc' 'calcfc' 'recalcfc=N' 'calcall' 등을 주로 사용한다. 나열한 키워드들은 모두 'force constant'와 관련이 있는 keyword들이다. 내가 찾고자 하는 구조의 force constant를 가지고 있으면 원하는 최적의 구조를 찾기까지의 시간이 짧아질 수 있기 때문이다.

 

3. Scan

최적화된 구조를 얻었다면, 그 구조의 에너지가 Reference RC가 된다. 이제 이 RC로부터 원하는 생성물을 얻기 위한 전이상태(Transition state)의 구조를 찾아야 하고, 그를 위해 특정한 두 원소 사이의 거리를 스캔한다. 키워드는 다음과 같다.

---

# opt=modredun(or modredundant)

Title

charge multiplicity

X Y Z

A B S N (-)n

---

'modredun'을 통해 우리가 원하는 A와 B 원소의 거리를 고정한다. 이후 구조를 최적화하고, n Å 만큼의 거리를 늘린다(좁힌다). 그리고 이를 N회 반복한다.

맨 처음에 얻은 RC는 TS보다 낮은 에너지를 가지므로, 해당 스캔의 에너지는 우상향한다. 그러다가 어느 순간에 에너지가 감소하는 구간이 발생하는데 (정확히는 한번 꺾이는 구간), 이 부분을 'Approxiate TS'라고 한다. Approximate가 붙는 이유는 아직 확실히 TS라고 보기 어렵기 때문.

물론 이후 계산을 진행하기 전에 현재 찾아낸 Approximate TS의 에너지가 타당한지의 여부를 판단해야 한다.

 

4. TS calculation

Approximate TS의 구조를 가져와서 frequency calculation을 진행한다. 이를 통해 진동수를 구할수가 있는데, 만약 단 하나의 '가상 진동수(imagninary frequency)'가 나타난다면, 해당 구조가 TS일 가능성이 매우 높아진다. 이 구조로부터, TS 구조의 최적화를 시작하게 된다. 일반적으로 내 실험실에서 사용하는 키워드는 다음과 같다.

---

# opt=(ts,nofreeze,rcfc,noeigentest)

---

먼저 frequency calculation을 진행한 후에 계산을 했으므로, 'rcfc' 키워드를 통해 force constant를 가져오게 된다. 이를 통해 초기의 force constant를 가져와 사용할 수 있고, 계산 속도의 증대가 이루어진다. 'nofreeze' 키워드는 이전의 chk 파일에 'modredun'에 의한 원소들의 구조 고정이 있었던 경우, 이를 모두 해제하라는 의미이다. 마지막으로 'noeigentest'는 eigentest를 진행하지 말라는 의미이다. 이 또한 계산 속도의 상승을 목적으로 작성한다.

TS 구조 계산이 끝나면 다시 frequency calculation을 진행한다. 이 과정에서 단 하나의 가상적인 진동수가 나오고, 에너지 역시 타당하다면 우리는 이를 일차적으로 'TS Structure'라고 볼 수 있게 된다.

 

5. IRC Calculations

TS 구조를 찾았으므로, 우리는 이로부터 실제 Reactant(RC)와 Product(PC)를 찾아야 한다. 이를 해결해주는 것이 IRC Calculation이다. IRC calculation은 TS로부터 시작해서 몇 포인트 앞, 몇 포인트 뒤의 구조를 계산한다. 그리고 이로부터 구조 최적화를 진행하면, 우리가 원하는 RC와 PC(아닐 수도 있다)가 나오게 되는 것이다. IRC Calculation의 키워드는 다음과 같다.

---

# irc=(rcfc,forward) / irc=(rcfc,reverse)

---

이전의 TS 구조에서 force constant를 'rcfc'를 통해 가져오게 되고, forward는 그로부터 앞방향, reverse는 그로부터 뒷방향으로 IRC를 진행하라는 의미이다. 일반적으로 forward에서 PC, reverse에서 RC가 나오나, 정방향 역방향이 그 자체로 '반응물'과 '생성물'을 의미하지는 않는다. 반드시 구조를 직접 눈으로 확인해야 한다.

만일 IRC가 잘 진행되지 않는다면, 'LQA'알고리즘을 사용하거나 'Stepsize=N'으로 스텝사이즈를 임의로 조정하여 문제를 해결할 수 있다. 일반적으로 Stepsize는 25정도로 조정하고 줄여가는 방식으로 진행한다.

IRC가 끝나게 되면, 해당 구조들로부터 구조 최적화를 시작한다. 이는 IRC의 최종 결과물을 알 수 있게 해주며, 우리는 이 결과를 통해 'RC'와 'PC'를 정확히 구분할 수 있게 된다.

6. Basis set Effect

모든 계산은 함수(Functional)와 그에 대한 기저 함수(basis set)로 이루어져 있다. 매우 다양한 종류의 함수와 기저함수가 존재하며, 각각의 공식 또한 천차만별이라 여기에 다 적기에는 여백이 부족하다. 일반적으로 내 실험실에서 사용하는 함수는 Unrestricted 'B3LYP'이고, 기저 함수로 'Def2-SVP'와 'Def2-TZVPP'를 사용한다.

이 기저 함수, 즉 Basis set에서 문제가 발생하는데, 그것은 바로 계산 시간의 문제이다. 높은 수준의 basis set은 에너지의 정밀도가 높지만, 계산의 시간이 매우 오래 걸리며, 낮은 수준의 basis set은 계산의 시간이 매우 빠르지만, 정밀도가 떨어지게 된다.

따라서, 우리는 구조 최적화 과정에서 낮은 수준의 basis set을 활용하여 계산 시간을 절약하고, 완성된 최종 구조물들의 single-point calculation만을 높은 수준의 basis set을 사용하여 계산할 수 있다. 이를 통해 상당히 많은 시간을 절약할 수 있다.

한편 이 결과 에너지의 차이가 발생하게 된다. 낮은 수준의 basis set과 높은 수준의 basis set의 에너지 차이를 Basis set effect라고 한다. 이 또한 반응의 기전에 따라 그 정도가 달라지므로, 반드시 자신의 반응에 적합한 함수를 찾고자 노력해야 한다.

 

Edited by 푸른삿포로
2022.01.14