오픈 소스화된 「MOPAC2016」을 Python으로 최대한 활용하기 (제1회: 도입 및 기초 계산 편)

서론: 계산 화학의 민주화와 MOPAC의 귀환

과거에 양자 화학 계산(분자 궤도법)은 고가의 워크스테이션과 더욱 비싼 상용 소프트웨어 라이선스를 보유한 일부 연구실만의 전유물이었습니다. 하지만 그 상식은 극적으로 변하고 있습니다.

그 상징적인 사건 중 하나가 반경험적 분자 궤도법 소프트웨어의 대명사인 「MOPAC [1]」의 오픈 소스화입니다.

오랫동안 상용 소프트웨어로 제공되었던 최신 버전 「MOPAC2016」이 Apache License 2.0 하에 공개되어, 누구나 무료로 이용할 수 있게 되었습니다 [2].

본 연재(총 3회)에서는 개발자에게 경의와 감사를 표하며, 이 개방된 강력한 계산 엔진을 Python에서 최대한 활용하는 방법을 소개합니다.

제1회 (이번 회차): MOPAC의 도입과 Python 스크립트를 통한 자동 계산 및 데이터 추출의 기초 -
제2회: 머신러닝 (ML) 응용을 겨냥한 대량 분자 데이터 생성 및 기술자 (Descriptor) 산출 -
제3회: PySide6를 이용한 GUI 앱 개발 (MOPAC 머신러닝 데이터셋 생성 앱 자작)

제1회인 이번 회차에서는 환경 구축부터 시작하여, Python의 화학 라이브러리인 「RDKit」과 연계해 분자의 3차원 구조 생성부터 HOMO-LUMO 에너지 산출까지 일괄적으로 수행하는 스크립트를 작성합니다.

「커맨드 라인은 조금 어렵다」는 화학 계열 분들이나, 「Python은 쓸 줄 알지만 화학 계산은 처음이다」라는 엔지니어 분들도 이해할 수 있도록 배경 지식을 포함하여 해설해 나가겠습니다.

제1장: 왜 지금, MOPAC인가?

반경험적 분자 궤도법이란

양자 화학 계산에는 크게 「ab initio (제1원리) 계산」과 「반경험적 (Semi-empirical) 계산」이 있습니다.

Gaussian이나 GAMESS로 대표되는 ab initio 계산은 슈뢰딩거 방정식을 엄격하게 풀려고 하기 때문에 정밀도가 높은 반면, 계산 비용이 방대합니다. 분자가 커지면 계산 시간은 지수 함수적으로 증가합니다.

반면 MOPAC가 채택하고 있는 반경험적 분자 궤도법은 일부 적분 계산을 실험값에 기반한 파라미터로 대체함으로써 계산 속도를 극적으로 가속화한 기법입니다.

「정밀도가 낮지 않을까?」라고 생각하실 수도 있지만, 유기 화학의 반응 설계나 대량의 분자 스크리닝, 혹은 머신러닝을 위한 초기 기술자 생성과 같은 용도에서는 현재에도 매우 강력한 도구입니다.

MOPAC2016과 PM7 메서드

이번에 사용하는 MOPAC2016에는 **PM7 (Parametric Method 7)**이라는 해밀토니안 (Hamiltonian, 계산 기법)이 탑재되어 있습니다.

이전 버전 (PM6 등)에 비해 수소 결합이나 분산력 (Van der Waals force)의 기술이 대폭 개선되어, 생체 분자나 유기 결정과 같이 비공유 결합성 상호작용이 중요한 계에서도 더욱 신뢰성 높은 결과를 얻을 수 있습니다.

이 정도 수준의 계산 엔진을 conda install 한 번으로 손에 넣을 수 있게 되었습니다. 이를 사용하지 않을 이유가 없습니다.

제2장: 개발 환경 구축

먼저 계산 환경을 갖추겠습니다. 이번에는 크로스 플랫폼 환경 구축 도구인 Conda (Miniconda 또는 Anaconda)를 이용합니다. Windows, macOS, Linux 어디에서든 동일한 절차로 구축이 가능합니다.

1. Python 환경 준비

아직 Conda를 도입하지 않았다면, Miniconda 공식 사이트 등을 참조하여 설치해 주세요.

설치 후, 터미널 (Windows의 경우 Anaconda Prompt 또는 PowerShell)을 열고 본 연재를 위한 가상 환경을 생성합니다.

# mopac-env라는 이름으로 Python 3.11 환경을 생성
conda create -n mopac-env python=3.11
conda activate mopac-env

2. MOPAC 설치

이 부분이 가장 핵심적인 하이라이트입니다. 이전에는 웹사이트에서 라이선스를 등록하고, 비밀번호를 받고, 바이너리를 다운로드하여 경로(Path)를 설정하는... 등의 번거로운 절차가 필요했지만, 현재는 conda-forge 채널에서 다음 명령어 하나로 완료됩니다.

conda install -c conda-forge mopac

설치가 끝났다면 정상적으로 도입되었는지 확인해 봅시다.

mopac --version

버전 정보가 표시되면 성공입니다. 이것으로 당신의 PC 안에 계산 화학 연구실이 세워졌습니다.

3. RDKit 설치

이번에는 계산 대상이 되는 분자의 초기 구조를 생성하기 위해, 화학 정보학용 라이브러리인 RDKit을 사용합니다. SMILES 표기법(문자열)으로부터 3차원 분자 구조를 생성하기 위해 필수적인 도구입니다.

conda install -c conda-forge rdkit

이것으로 필요한 도구는 모두 갖춰졌습니다.

나머지는 여러분의 실행 환경에 따라 JupyterLab 등을 적절히 설치하시면 됩니다.

제3장: Python을 이용한 자동 계산 전략

여기서부터는 실제로 Python 코드를 작성해 보겠습니다.

이번 목표는 "SMILES(문자열)를 입력하는 것만으로, 가장 안정적인 입체 형태(Conformation)를 찾아내어, MOPAC로 구조 최적화와 궤도 에너지 계산을 수행하고, 결과를 표시하는" 일련의 흐름을 자동화하는 것입니다.

수동으로 진행할 경우, "Avogadro 등의 소프트웨어로 분자를 그리기" $\rightarrow$ "MOPAC용 입력 파일 만들기" $\rightarrow$ "명령어 실행하기" $\rightarrow$ "텍스트 에디터로 로그 파일(Log file)을 열어 수치 찾기"라는 절차를 거치게 되지만, 이 모든 것을 Python이 대신하게 합니다.

스크립트 구성은 다음의 4단계입니다.

RDKit을 이용한 3D 구조 생성 및 형태 탐색 $\rightarrow$ MOPAC 입력 파일(.mop) 생성 $\rightarrow$ MOPAC 실행 (subprocess) $\rightarrow$ 출력 파일(.out) 분석 및 결과 표시

각 단계에 대한 기술적인 포인트를 설명하겠습니다.

Step 1: RDKit을 이용한 형태 탐색의 중요성

분자, 특히 유기 분자는 단일 결합의 회전에 따라 다양한 형태(Conformation)를 취합니다. MOPAC 등의 구조 최적화는 어디까지나 "입력된 구조로부터 가장 가까운 안정점(극소점, Local Minimum)"을 찾는 작업입니다.

만약 뒤틀린 불안정한 구조를 입력하면, 불안정한 상태 그대로 계산이 수렴해 버려 올바른 물성치를 얻을 수 없습니다.

따라서 이번에는 RDKit을 사용하여 여러 가지 형태(Conformer)를 생성하고, 힘의 장(Force field) 계산(MMFF)으로 대략적인 에너지를 추정하여, "가장 안정할 것으로 예상되는 구조"를 선별한 뒤 MOPAC에 전달하는 처리를 구현합니다.

코드 내에서는 AllChem.EmbedMultipleConfs와 AllChem.MMFFOptimizeMoleculeConfs를 사용합니다. 이를 통해 계산의 성공률과 신뢰성을 격상시킬 수 있습니다.

Step 2: 연결 작업 (Linked Job)의 활용

MOPAC에는 **연결 작업 (Linked Job)**이라는 편리한 기능이 있습니다. 이는 하나의 입력 파일로 여러 계산을 연속해서 수행하게 하는 기능입니다.

이번에는 다음과 같은 2단계 계산을 한 번에 실행합니다.

• 제1 작업: 구조 최적화 (Geometry Optimization)

  • 키워드: PM7 PRECISE STATIC EF
  • 의미: PM7 방법을 사용하여 고정밀(PRECISE) 및 정적(STATIC)으로, 고유 벡터 추적법(EF)을 통해 구조를 최적화한다.

• 제2 작업: 전자 상태 및 진동 계산 (Electronic/Vibrational Analysis)

  • 키워드: PM7 PRECISE FORCE OLDGEO
  • 의미: PM7 방법을 사용하여, 제1 작업에서 최적화된 구조(OLDGEO)를 이어받아, 힘의 상수(Force constant)나 진동(FORCE)을 계산한다.

입력 파일 내에서 &라는 행을 삽입함으로써 이러한 작업들을 연결할 수 있습니다. 이를 통해 "최적화된 좌표를 복사하여 다시 계산"하는 번거로움을 줄일 수 있습니다.

Step 3 & 4: 실행 및 분석

실행은 Python 표준 모듈인 subprocess로 수행합니다.

분석은 정규 표현식(re 모듈)을 사용합니다. MOPAC의 출력 파일은 역사적 경위로 인해 고정 길이 포맷(Fixed-length format)에 가까운 텍스트 형식이지만, 버전에 따라 공백의 수가 미세하게 달라질 수 있으므로 유연한 정규 표현식으로 수치를 포착하는 방법을 채택했습니다.

특히 HOMO/LUMO 에너지는 출력 파일 끝부분에 있는 요약 섹션(HOMO LUMO ENERGIES)에서 가져오지만, 경우에 따라 이온화 포텐셜(IONIZATION POTENTIAL)로부터 추정하는 로직도 포함함으로써 분석의 견고성(Robustness)을 높입니다.

제4장: 구현 코드 완전판

그럼, 실제로 동작하는 코드를 소개하겠습니다. 다음 코드를 run_mopac.py 등의 이름으로 저장해 주세요.

이번에는 예시로 페놀(Phenol) 계산을 수행합니다.

import sys
import os
import subprocess
...

제5장: 실행 결과와 해설

실행 방법

위의 코드를 저장한 디렉토리에서 다음 명령어를 실행합니다.

python run_mopac.py

출력 결과 예시

다음과 같은 로그가 터미널에 표시되면 성공입니다.

[Oc1ccccc1]의 구조 생성과 MMFF 컨포머 탐색(Conformer Search)을 시작합니다...
최안정 컨포머의 MMFF 에너지: 8.2341 kcal/mol
MOPAC 인풋 파일(Input File)을 생성했습니다: phenol-0.mop (연결 작업)
...

결과의 화학적 의미

여기서 얻은 **HOMO(최고 점유 궤도, Highest Occupied Molecular Orbital)**와 **LUMO(최저 비점유 궤도, Lowest Unoccupied Molecular Orbital)**는 분자의 반응성을 이해하는 데 있어 가장 중요한 지표입니다.

• HOMO 에너지가 높음: 전자를 주기 쉬워 산화되기 쉽다 (친핵성(Nucleophilicity)이 높다).
• LUMO 에너지가 낮음: 전자를 받기 쉬워 환원되기 쉽다 (친전자성(Electrophilicity)이 높다).
• HOMO-LUMO Gap: 이 값이 작을수록 전자가 들뜨기(Excitation) 쉬우며, 전도성이나 광흡수 특성과 관련이 있습니다 (색이 나타나기 쉬운 등).

페놀의 결과를 보면, 벤젠 고리에 OH기(전자 공여체, Electron Donating Group)가 붙어 있기 때문에 벤젠 단독일 때보다 HOMO가 높아져, 친전자성 치환 반응(니트로화 등)이 일어나기 쉽다는 점 등을 수치로부터 읽어낼 수 있습니다.

생성된 파일 확인하기

스크립트를 실행하면 디렉토리 내에 phenol-0.mop (입력), phenol-0.out (출력), 그리고 phenol-0.arc (아카이브) 등이 생성됩니다.

.out 파일을 텍스트 에디터로 열어보세요. 이번에 Python으로 추출한 정보 외에도 원자의 전하(Charge), 쌍극자 모멘트(Dipole), 열역학적 생성열(Heat of Formation) 등 보물창고와 같은 데이터가 기록되어 있는 것을 알 수 있을 것입니다.

요약 및 다음 예고

이번에는 MOPAC2016의 오픈 소스화를 계기로, Python을 사용하여 페놀의 분자 궤도 계산을 자동화하는 방법을 소개했습니다.

포인트를 복습해 보겠습니다.

MOPAC는 Conda로 쉽게 설치할 수 있다. RDKit으로 컨포머 탐색을 수행한 뒤 MOPAC에 전달하면 안전하다. 연결 작업(&)을 사용하면 최적화와 물성 계산을 한 번에 수행할 수 있다. 출력 파일은 정규 표현식(Regular Expression)으로 해석할 수 있다.

단 수십 줄의 코드이지만, 여기에는 계산 화학 자동화의 정수가 담겨 있습니다. 이 SMILES_STR 변수를 바꾸는 것만으로 모든 유기 분자의 계산이 가능합니다.

하지만 데이터 사이언스에 활용하기에는 궤도 에너지 데이터만으로는 부족하겠죠?

다음 제2회에서는 이 스크립트를 확장하여, "하나의 분자로부터 많은 물성치를 산출하고 머신러닝용 데이터셋(CSV)을 작성하는" 방법에 대해 해설하겠습니다.

계산 화학과 인포매틱스(Informatics)의 융합, 이제부터가 진짜 시작입니다. 기대해 주세요!

(저자 주: 본 기사의 코드는 MOPAC2016 및 RDKit의 특정 버전에서 동작 확인을 거쳤습니다. 환경에 따라 수치의 미세한 차이나 경고(Warning)가 발생할 수 있습니다)

Discussion