오픈 소스화된 「MOPAC2016」을 Python으로 최대한 활용하기 (제2회: MI용 데이터셋 작성)

서론: 점으로부터 선, 그리고 면으로

[지난 회차(제1회)]에서는 MOPAC2016의 설치부터 Python을 사용하여 가장 기본적인 계산(HOMO-LUMO 에너지 산출)을 수행하는 방법까지를 해설했습니다. 자신의 PC 안에서 명령어 하나로 양자 화학 계산(Quantum Chemistry Calculation)이 돌아가는 모습에 감동하신 분들도 계실 것입니다.

하지만 현대의 화학 연구, 특히 머티리얼즈 인포매틱스 (Materials Informatics, MI) 문맥에서 단 하나의 물성치(Descriptor)만으로 충분한 경우는 드뭅니다.

머신러닝 (Machine Learning) 모델에 '분자의 개성'을 올바르게 가르치기 위해서는 전자적 성질, 입체적 성질, 열역학적 안정성 등 다각적인 정보(Feature Vector)가 필요합니다.

제2회인 이번 테마는 「기술자(Descriptor)의 망라적인 자동 산출」 입니다.

MOPAC2016은 적절한 키워드를 지정하고 출력 파일을 깊이 있게 분석함으로써 놀라울 정도로 많은 물리량을 제공합니다. 이번에는 지난 회차의 코드를 대폭 확장하여, 다음과 같은 '보물산'을 자동으로 CSV 파일에 담아내는 데이터 산출기를 구축합니다.

**용매 효과 (COSMO 법)**를 고려한 에너지와 표면적·부피 -
정적 분극률 (Static Polarizability) 및 쌍극자 모멘트 -
열역학적 물성 (엔탈피, 엔트로피, 열용량, 생성열) -
HOMO/LUMO 근방의 궤도 준위 (HOMO-5 〜 LUMO+5)

이것이 가능해지면, 나머지는 분자 리스트를 루프(Loop)시키는 것만으로 당신만의 '고품질 화학 데이터베이스'가 완성됩니다. 그럼 Python과 MOPAC의 심연한 연계의 세계로 발을 들여놓겠습니다.

제1장: MI 기술자로서의 물리량과 MOPAC 키워드

코드 해설에 들어가기에 앞서, "왜 그 계산이 필요한가?"라는 화학적 전략을 세웁니다. 이번에 타겟으로 삼을 물리량과 그것을 MOPAC에서 출력시키기 위한 키워드에 대해 해설하겠습니다.

1. 용매 효과의 도입 (EPS, COSMO)

유기 반응의 대부분은 용액 중에서 일어납니다. 진공(기상)에서의 계산 결과만으로는 극성 분자의 거동을 오판할 가능성이 있습니다.

MOPAC에는 **COSMO 법 (Conductor-like Screening Model)**이라는 용매 효과를 근사적으로 도입하는 기능이 갖춰져 있습니다.

키워드: EPS=78.4

• 이것은 유전율 ($\epsilon$)을 지정하는 것입니다. 78.4는 '물'의 값입니다. 이것을 지정하는 것만으로 수용액 중에서의 안정 구조나 에너지를 계산해 줍니다. 유기 용매라면 아세토니트릴(37.5)이나 DMSO(46.7) 등을 지정합니다.

2. 정적 분극률 (STATIC)

분극률 (Polarizability)은 외부 전기장에 대한 전자 구름의 '왜곡되기 쉬운 정도'를 나타내는 지표입니다. 이는 분자 간 상호작용(분산력)이나 광학적 특성(굴절률 등), 나아가 화학적 반응성(친핵 공격을 받기 쉬운 정도)과도 상관관계가 있는 중요한 물성입니다.

통상적인 구조 최적화만으로는 출력되지 않기 때문에 전용 키워드를 추가합니다.

키워드: STATIC

• 이것을 지정하면 정적 분극률 텐서와 그 평균값을 계산합니다.

3. 열역학적 물성과 진동 해석 (FORCE, THERMO)

분자가 '존재할 수 있는지'를 알기 위해서는 에너지만으로는 불충분하며 자유 에너지가 필요합니다. 진동 해석을 수행함으로써 영점 진동 에너지 (Zero-Point Energy, ZPE) 보정을 포함한 생성 엔탈피나 엔트로피를 산출할 수 있습니다.

키워드: FORCE THERMO

• FORCE: 힘의 상수 행렬 (Hessian)을 계산하여 진동수를 구합니다. -
• THERMO: 진동수를 바탕으로 지정 온도(표준은 298K)에서의 열역학 함수를 산출합니다.

4. 연결된 작업 (Linked Job)의 재설계

이것들을 한 번에 수행하기 위해 이번에도 연결된 작업 (Linked Job)을 사용합니다. 다만, 제1회보다 고도화된 구성입니다.

Step 1: PM7 PRECISE EF STATIC EPS=78.4

• 수중에서 구조 최적화를 수행하고, 그 구조에서의 분극률 (STATIC)도 계산한다.

Step 2: PM7 PRECISE FORCE THERMO EPS=78.4 OLDGEO

• Step 1의 구조 (OLDGEO)를 사용하여...

)를 이어받아, 수중에서 진동 해석(FORCE)과 열역학 계산(THERMO)을 수행한다.

• Step 1의 구조(OLDGEO)를 사용하여

이로써 구조, 전자, 열, 응답 특성을 한 번에 취득할 수 있습니다.

제2장: Python을 이용한 「정규 표현식 (Regular Expression)」 처리

MOPAC의 출력 파일(.out)은 사람이 읽기 위한 리포트 형식(고정 길이 텍스트에 가까운 형식)으로 작성되어 있습니다. 이를 머신러닝용 CSV(구조화 데이터)로 변환하려면, **정규 표현식 (Regular Expression)**을 구사하여 필요한 수치를 추출해야 합니다.

이번 코드에서는 Python의 re 모듈을 풀 활용하고 있습니다. 몇 가지 테크닉을 소개하겠습니다.

1. 부동 소수점 수의 캡처

MOPAC의 수치는 "공백"으로 구분되어 있습니다.

([-+]?\d+\.\d+)

이는 "부호(+ 또는 -)가 있어도 되고 없어도 되며, 숫자가 이어지고, 점(dot)이 있으며, 다시 숫자가 이어지는" 패턴, 즉 실수를 포착하기 위한 정형구입니다.

2. 여러 줄에 걸친 블록의 탐색

예를 들어 진동수나 궤도 에너지는 여러 줄에 걸쳐 출력됩니다.

re.DOTALL 옵션을 사용하면 . (도트)가 개행 코드(newline)에도 매치되도록 하여,

r"EIGENVALUES\s*\n(.*?)(?=\n\s*[A-Z])"

와 같이 작성함으로써, "EIGENVALUES라는 헤더부터 다음 대문자 헤더가 나오기 전까지의 모든 텍스트"를 통째로 가져올 수 있습니다.

3. 표 형식 데이터의 특정

"쌍극자 모멘트 (Dipole Moment)"나 "분극률 (Polarizability)"은 유사한 단어가 여러 번 등장하기 때문에, 단순한 검색으로는 "초기 구조의 값"이나 "최적화 도중의 값"을 잘못 가져올 리스크가 있습니다.

따라서 "표의 경계선"이나 "단위 (Debye, ANG.**3)"를 행 끝에 가지는 행 등, 유니크한 특징을 가진 주변 문자열을 앵커(Anchor, 닻)로 삼아 목표로 하는 행을 핀포인트로 특정합니다.

제3장: 구현 코드의 상세 해설

그럼 이번 코드(run_mopac_mi.py)를 살펴보겠습니다.

페놀 (Phenol)을 예로 들어, 앞서 소개한 모든 기능을 탑재한 스크립트입니다.

1. 라이브러리 임포트 및 설정

이번에는 데이터 처리의 주역으로서 pandas를 도입했습니다 (설치되어 있지 않다면 pip install pandas 명령어로 설치해 두시기 바랍니다). 추출한 데이터를 딕셔너리(dict)에 저장하고, 마지막에 DataFrame으로 변환하여 CSV로 만드는 것이 현대적인 정석입니다.

import sys
import os
import subprocess
...

2. RDKit을 이용한 초기 구조 생성 (전회와 동일)

이 부분은 제1회와 마찬가지로 RDKit의 MMFF (Force Field)를 사용하여 타당한 초기 배좌(conformation)를 찾아냅니다. 이 부분이 적절하지 않으면 MOPAC 계산이 수렴하지 않거나 국소해 (Local Minimum, 이상한 구조)에 빠질 수 있기 때문에, MI용 데이터셋 작성 시에는 소홀히 할 수 없는 부분입니다.

(※ 코드의 "1. RDKit을 이용한..." 부분은 제1회 참조를 위해 생략)

3. 고기능 인풋 파일 작성

이 부분이 이번의 하이라이트 중 하나입니다. 문자열 결합으로 입력 파일을 만들지만, 키워드 설정에 주목해 주세요.

# ==========================================
# 2. MOPAC 인풋 파일 작성 (2 Step Linked Job)
# ==========================================
...

이로써 "수중의 페놀"을 정밀 계산할 준비가 완료되었습니다.

4. 출력 파일의 철저한 분석

계산 실행 후 (subprocess.run 이후), .out 파일의 내용을 분석합니다. 여기서부터의 코드는 그대로 라이브러리로 사용할 수 있을 정도의 완성도를 갖추고 있습니다.

A. 기본 데이터와 쌍극자 모멘트

쌍극자 모멘트 (Dipole)는 출력 형식이 경우에 따라 달라지기 때문에 (점군 대칭성이 있는 경우 등 표 형식이 되는 경우), 2단계로 나누어 취득합니다.

# Dipole
dip_table_match = re.findall(r"SUM\s+[-+\d\.]+\s+[-+\d\.]+\s+[-+\d\.]+\s+([-+\d\.]+)", text)
if dip_table_match:
...

B. 궤도 에너지 주변 탐색 (HOMO/LUMO ± 5)

MI(Machine Learning)에서는 HOMO/LUMO의 에너지 차이(Gap)뿐만 아니라, "그 주변에 궤도가 얼마나 밀집되어 있는가(DOS, Density of States)"가 중요한 정보가 될 수 있습니다.

따라서 HOMO를 중심으로 전후 5개의 궤도를 추출하는 로직을 구현했습니다.

# 점유 궤도의 수(NO. OF FILLED LEVELS)를 가져와서, 그곳을 기준(HOMO)으로 삼음
filled_match = re.findall(r"NO\. OF FILLED LEVELS\s*=\s*(\d+)", text)
# 모든 고유값 리스트를 가져옴
...

이 루프 처리를 통해 데이터셋의 열(Column)로서 Energy_HOMO뿐만 아니라 Energy_HOMO-1 등이 자동으로 생성됩니다.

C. 정적 분극률 (Polarizability)

STATIC 키워드 특유의 출력을 포착합니다. 이 부분은 매우 혼란스러운 출력이 많은 구간이므로, Average Polarizability from:이라는 문자열 블록을 특정하고 나서 그 안의 ANG.**3 행을 찾는 신중한 방식으로 구현했습니다.

# 블록을 특정
static_block_match = re.search(r"Average Polarizability from:.*?ANG\.\*\*3", text, re.DOTALL)
if static_block_match:
...

D & E. COSMO 부피와 열역학 데이터

COSMO법에 의한 분자 표면적(Area)과 부피(Volume), 그리고 THERMO 계산에 의한 열역학량을 가져옵니다.

특히 생성열(Heat of Formation)은 원자의 결합 에너지와 비교하여 분자의 안정성을 논의하는 데 있어 가장 기본적인 값입니다.

# 열역학 데이터 (298K)
thermo_298_block = re.search(r"^\s*298\.00.*?\n\s+TOT\.\s+([^\n]+)", text, re.DOTALL | re.MULTILINE)
if thermo_298_block:
...

F. 진동수 체크

마지막으로 진동수(Frequencies)를 가져옵니다.

구조 최적화가 올바르게 수행되어 극소점(안정 구조)에 도달했다면, **모든 진동수는 실수(양수)**가 됩니다. 만약 음수(MOPAC에서는 출력상 마이너스가 붙는 허진동)가 있다면, 그것은 전이 상태(Saddle point)임을 시사합니다.

여기서는 최소 진동수(Min_Freq_cm-1)를 기록해 둠으로써, 나중에 "이 데이터는 정말 안정 구조인가?"를 스크리닝할 수 있도록 합니다.

5. CSV 저장 및 정렬

마지막으로 수집한 데이터를 Pandas DataFrame으로 변환하고, 보기 쉬운 순서로 정렬하여 CSV로 저장합니다.

여기서는 orb_sort_key라는 함수를 정의하여, 궤도 에너지 열(HOMO-5, ..., HOMO, LUMO, ...)이 에너지 순서대로 나열되도록 고안했습니다.

제4장: 실행 및 결과 해석

실행 방법

이 스크립트를 실행하면 디렉토리에 phenol_data.csv가 생성됩니다.

python run_mopac_mi.py

생성되는 데이터 예시 (발췌)

다음과 같은 데이터가 출력됩니다 (Jupyter 환경에서의 아웃풋).

단 하나의 분자(페놀)로부터 이토록 풍부한 정보를 추출할 수 있었습니다.

이것이 머신러닝에서의 "설명 변수(X)"가 됩니다.

요약: 데이터셋 생성 엔진의 완성

이번에는 MOPAC2016과 Python을 조합하여 실용적인 MI용 기술자(Descriptor)를 자동 산출하는 시스템을 구축했습니다.

• 키워드 선정: EPS, STATIC, THERMO를 조합함으로써 정적인 에너지뿐만 아니라 환경 응답 및 열물성까지 취득할 수 있다.
• 정규 표현식 해석: 복잡한 출력 파일도 패턴 매칭을 구사하면 구조화된 데이터로 변환할 수 있다.
• 확장성: 이 코드의 SMILES_STR을 리스트화하여 루프를 돌리면, 방대한 수의 분자 리스트 데이터라도 자는 동안 데이터셋이 완성된다.

이번에 작성한 코드는 말하자면 "분자 데이터 생성 엔진"의 핵심 부품입니다.

다음 회, 최종 회가 될 제3회에서는 이 핵심 부품을 통합하여, GUI (Graphical User Interface) 애플리케이션화에 도전합니다.

검은 화면(터미널)을 어려워하는 동료나 실험 화학자들을 위해, 분자 리스트 파일을 선택하기만 하면 이 고도의 계산이 백그라운드에서 실행되고 결과가 CSV로 툭 튀어나오는——그런 친절한 설계의 도구를 「PySide6」를 사용하여 제작합니다.

Python을 이용한 화학 도구 개발, 드디어 집대성입니다. 기대해 주세요!

부록: 이번에 사용한 소스 코드 전체

import sys
import os
import subprocess
...

Discussion