[LNP] Step4_create_top_hexagon
너무 힘들어 ㅠㅠ......
그래도 해야지.....
강한 여성은 해낼 것이다....!
## 이전 글 ##
2026.10.30 - [무한도전!!/MD 시뮬레이션 돌리자!!] - [LNP] Step3_hexagon_channels
[LNP] Step3_hexagon_channels
이제 드디어 Step3에 도달했다.... 또 열심히 달려보자... Step3는 Step2에서 만든 채널 구조(m.gro)를 복제한 뒤, 그 안에서 육각형 배열을 이루는 채널 묶음만 골라내는 단계를 진행한다.이는 크게 세
www.cheonggwan.com
이 노트북은 Step3에서 잘라낸 Hexagon.pdb를 읽어서 GROMACS 토폴로지 파일 hexagon.top을 만든다.
아래 import 6개는 각각 이 과정의 한 부분씩을 담당하고 있다.
# 라이브러리 불러오기
import MDAnalysis as mda
import numpy as np
import re
from vermouth.forcefield import ForceField
from vermouth.gmx.itp_read import read_itp
import networkx as nx
import MDAnalysis as mda 는 구조 파일 읽기 위한 import이다.
분자동역학(MD) 구조와 trajectory 파일(PDB, GRO, XTC 등)을 파이썬 객체로 불러온다.
from vermouth.forcefield import ForceField, from vermout.gmx.itp_read import read_itp 는 ITP를 파싱하기 위한 라이브러리이다. vermouth는 Martini 개발 그룹이 만든 라이브러리로, martinize2와 polyply 내부에서 쓰인다.
이는 GROMACS 토폴로지(.itp)를 파이썬 구조로 다룰 수 있게 해준다.
- ForceField : 여러 분자 정의(block)를 담는 컨테이너. forcefield.blocks 는 {분자이름: Block} 형태의 딕셔너리다.
- read_itp(lines, forcefield) : .itp 파일의 줄 목록을 파싱해서 [ moleculetype ] 하나하나를 block으로 만들어 forcefield에 넣어준다.
# Step 1: Generate topology
해당 단계는 Hexagon.pdb에 맞는 GROMACS 토폴로지 파일(hexagon.top)을 만드는 과정이다.
1. itp_to_forcefield() 함수 정의
def itp_to_forcefield(itp):
"""
Returns a forcefield with all the molecule in the itp added.
The itp should specify the path to the ITP file.
"""
path = itp
with open(path, 'r') as f:
lines = f.readlines() # 1. 파일을 줄 단위 리스트로 읽음
forcefield = ForceField(name='custom') # 2. 빈 forcefield 생성
read_itp(lines, forcefield) # 3. 파싱해서 forccefield 에 채워 넣음
return forcefield # 4. 반환
목적: ITP 파일 경로를 받아서, 그 안의 모든 분자가 들어 있는 vermouth ForceField 객체를 돌려준다.
- read_itp는 파일 경로가 아니라 줄 리스트를 입력으로 받는다.
- 분자 정의를 담을 빈 상자이다.
- 반환값이 없고, 넘겨준 forcefield를 직접 수정한다. ITP의 [ moleculetype ] 섹션마다 Block을 하나씩 만들어 forcefield.blocks 에 등록한다.
- 다 채워진 forcefield 를 반환한다.
분자 하나 = 그래프 하나 / bead는 노드, 결합은 edge가 된다.
2. ITP 6개를 읽고 하나의 딕셔너리로 합치기
# 2-1 Import the itps which are sub dependencies
ff_rna = itp_to_forcefield('ITPs/RNA.itp')
ff_lipids = itp_to_forcefield('ITPs/martini_v3.0.0_phospholipids_v1.itp')
ff_lisbeth = itp_to_forcefield('ITPs/ILs.itp')
ff_sterols = itp_to_forcefield('ITPs/martini_v3.0_sterols_v1.0.itp')
ff_solvents = itp_to_forcefield('ITPs/martini_v3.0.0_solvents_v1_BMHB.itp')
ff_ions = itp_to_forcefield('ITPs/martini_v3.0.0_ions_v1.itp')
# 2-2 Combine all the blocks in one big dictionary (a block is a molecule in our case)
all_blocks = {}
for forcefield in [ff_rna, ff_lipids, ff_lisbeth, ff_sterols, ff_solvents, ff_ions]:
for name, block in forcefield.blocks.items():
if name in all_blocks.keys():
print(f'A double definition was found for {name}, this often results in issues!')
all_blocks[name] = block
#all_blocks['RNA'].nodes(data=True)
2-1 파일별로 forcefield 만들기
ff_rna = itp_to_forcefield('ITPs/RNA.itp')
ff_lipids = itp_to_forcefield('ITPs/martini_v3.0.0_phospholipids_v1.itp')
ff_lisbeth = itp_to_forcefield('ITPs/ILs.itp')
ff_sterols = itp_to_forcefield('ITPs/martini_v3.0_sterols_v1.0.itp')
ff_solvents = itp_to_forcefield('ITPs/martini_v3.0.0_solvents_v1_BMHB.itp')
ff_ions = itp_to_forcefield('ITPs/martini_v3.0.0_ions_v1.itp')
각 파일에 들어있는 분자 수( [ moleculetype ] 개수 )는 다음과 같다.
| 변수 | 파일 | 분자 수 | 내용 |
| ff_rna | RNA.itp | 1 | RNA |
| ff_lipids | phospholipids_v1 | 113 | DOPC, DSPC, DOPE 등 인지질 |
| ff_lisbeth | ILs.itp | 36 | 이온화 지질 |
| ff_sterols | sterols_v1.0 | 2 | CHOL, CHOX |
| ff_solvents | solvents_v1_BMHB | 59 | 물(W), 유기용매 등 |
| ff_ions | ions_v1 | 20 | NA, CL, 버퍼 bead 등 |
2-2 모든 Block을 하나로 합치기
all_blocks = {}
for forcefield in [ff_rna, ff_lipids, ff_lisbeth, ff_sterols, ff_solvents, ff_ions]:
for name, block in forcefield.blocks.items():
if name in all_blocks.keys():
print(f'A double definition was found for {name}, this often results in issues!')
all_blocks[name] = block
#all_blocks['RNA'].nodes(data=True)- 바깥 루프: forcefield 6개를 하나씩 돈다.
- 안쪽 루프: forcefield.blocks는 {분자이름: Block} 딕셔너리다. .items() 로 이름과 Block을 쌍으로 꺼낸다.
결과 : all_blocks는 약 231개 분자가 담긴 딕셔너리다.
마지막 줄에 주석 처리된 코드는 디버깅용으로, 주석을 풀면 RNA의 bead 목록을 볼 수 있다.
3. 분자별 bead 목록 (molecule_identifiers) 만들기
molecule_identifiers = {}
for molecule in all_blocks.keys():
molecule_identifiers[molecule] = []
for key, values in all_blocks[molecule].nodes(data=True):
molecule_identifiers[molecule].append((values['resname'], values['atomname']))
#molecule_identifiers
목적은 복잡한 Block 그래프에서 매칭에 필요한 정보 두 가지 (resname, atomname) 만 뽑아서 단순한 리스트로 만든다.
- 핵심 문법: block.nodes(data=True)
- vermouth의 Block은 networkx.Graph 를 상속한 클래스라서, networkx의 메서드를 그대로 쓸 수 있다.
- data=True 를 주면 (노드키, 속성딕셔너리) 쌍을 돌려준다. 여기서는 노드키(key)는 쓰지 않고 속성(values)만 쓴다.
- 노드는 ITP [ atoms ] 에 적힌 순서대로 나온다 (networkx는 노드를 추가한 순서를 유지한다). 그래서 리스트 순서가 PDB의 원자 순서와 맞다.
4. PDB 읽기
# Read in the GRO file
universe = mda.Universe('Hexagon.pdb') ## This pdb you generate from step3, using VMd or similar tool, to manually select the desired core- VMD 같은 도구로 원하는 코어 영역만 잘라 저장한 PDB를 MDAnalysis Universe 로 불러온다.
- 주석에는 "GRO file" 이라고 적혀있지만 실제로는 PDB를 읽는다.
- 이후 셀에서는 universe.residues (residue를 파일에 나온 순서대로 순회), universe.atoms(전체 bead 수) 만 사용한다.
5. residue를 분자로 매칭하기
# Start matching the resnames
# This is not a very generic approach, but it can be expanded. The issue
# is that there is no true solution to this problem. For certain molecules
# might be subgroups of other molecules, meaning it is to a certain degree
# an ambigous problem. Here we map small segments over long segments for
# simplicity, this might or might not be a problem for you system.
# This can largely be circumvented by making sure that all resnames are
# unique for every molecule type.
matches = []
active_molecule = []
in_RNA = False
for residue in universe.residues:
# Matching by residue name if the molecule
# is a residue.
if residue.resname in all_blocks.keys():
in_RNA = False
matches.append(residue.resname)
# The naming of IONS is weird as they are separated
# molecules which have the same resname
elif residue.resname == 'ION':
in_RNA = False
matches.append(residue.atoms.names[0])
# The RNA or any polymer is an issue for
# their residues are not molecules. Here
# I wrote a basic detection for RNA. But
# this might need to be changed for other
# polymers. I only match by length, this
# might be lazy but for now I think it
# will work.
elif residue.resname in ['A', 'U']:
if in_RNA == False:
in_RNA = len(residue.atoms)
else:
in_RNA += len(residue.atoms)
if in_RNA == len(molecule_identifiers['RNA']):
in_RNA = False
matches.append('RNA')
else:
in_RNA = False
print(f'UNKNOWN RESIDUE {residue.resname}')
break
#print(matches)
PDB는 residue 단위 정보만 있고, 그게 어떤 분자(moleculetype)인지는 적혀 있지 않다.
대부분은 residue 하나가 분자 하나지만, RNA와 이온의 경우에는 그렇지 않다.
RNA는 A, U, ... 와 같은 residue 여러 개가 모여 분자 하나가 되고, 이온은 resname은 모두 ION 인데, 실제 분자 이름은 NA, CL 처럼 다르다.
- 변수 준비
matches = [] # 결과: 분자 이름을 PDB 순서대로 쌓는 리스트
active_molecule = [] # ⚠️ 선언만 하고 쓰지 않음 (남은 코드)
in_RNA = False # RNA를 읽는 중인지 표시하는 플래그 겸, 지금까지 센 bead 수
- 4가지 분기
residue를 하나씩 보며아래 순서대로 검사한다.
| 분기 | 조건 | 동작 | Hexagon 예시 |
| ① 일반 분자 | resname이 all_blocks에 있음 | resname을 그대로 분자 이름으로 추가 | MC3H, CHOL, W |
| ② 이온 | resname == 'ION' | 첫 원자 이름을 분자 이름으로 추가 | 원자 이름 CL → 'CL' |
| ③ RNA | resname이 'A' 또는 'U' | bead 수를 누적하다가 RNA 한 분자 크기에 도달하면 'RNA' 추가 | 아래 설명 |
| ④ 모름 | 위에 모두 해당 안 됨 | UNKNOWN RESIDUE 출력 후 루프 중단 | - |
* RNA 판별 (자세히)
elif residue.resname in ['A', 'U']:
if in_RNA == False:
in_RNA = len(residue.atoms) # 첫 뉴클레오타이드: 카운터 시작
else:
in_RNA += len(residue.atoms) # 이어지는 뉴클레오타이드: bead 수 누적
if in_RNA == len(molecule_identifiers['RNA']): # RNA ITP의 bead 수(357)와 같으면
in_RNA = False
matches.append('RNA') # RNA 1개 완성
* 한계
- bead 수만 비교한다. 서열은 확인하지 않기 때문에, 길이만 같으면 다른 RNA도 통과한다.
- 'A', 'U' 만 인식한다. G나 C가 들어간 RNA를 쓰면 ④로 분류되어 멈춘다.
- RNA가 잘려서 bead 수가 357에 못 미치면 'RNA' 가 추가되지 않는다.
- in_RNA 가 False(플래그)와 정수(카운터) 두 역할을 동시에 한다.
6. 연속된 같은 분자 묶어 세기
successive_elements = []
previous_match = None
counter = None
for match in matches:
if match == previous_match: # 직전과 같은 분자 → 개수 +1
counter += 1
else: # 다른 분자가 나옴 →
if previous_match is not None:
successive_elements.append((previous_match, counter)) # 직전 묶음 저장
counter = 1 # 새 묶음 시작
previous_match = match
successive_elements.append((match, counter)) # 루프가 끝난 뒤 마지막 묶음 저장
긴 matches 리스트를 (분자 이름, 연속 개수) 쌍으로 압축하는 것이 목적이다.
"전체 개수"가 아니라 "연속 개수"인 이유는 [ molecules ] 는 좌표 파일의 순서를 그대로 따라야 하기 때문이다.
그래서 같은 분자라도 중간에 다른 분자가 끼어 있으면 따로 적어야 한다. 예를 들어 MC3H 980 한 줄로 합치면 안되고, 채널마다 MC3H 140을 7번 적어야 한다.
루프 밖에는 마지막 append가 있는데 이는 루프 안에서 "다음 분자가 바뀔 때"의 직전에 묶음을 저장하기 때문이다.
그래서 맨 마지막 묶음은 루프 안에서 저장될 기회가 없어서 따로 저장한다.
7. 원자 수로 검증하기
# Quick sanity check to see if the amount of atoms
# based on the segment counting is indeed the equal
# to the total amount of atoms in the PDB/GRO.
amount_of_atoms = 0
for element in successive_elements:
amount_of_atoms += len(molecule_identifiers[element[0]])*element[1]
all_fine = amount_of_atoms == len(universe.atoms)
print(f'The amount of atoms based on the segments in correct: {all_fine}')
토폴로지대로 계산한 bead 수와 PDB의 실제 bead 수가 같은지 확인하는 것이 목적이다.
이 검증으로 잡을 수 있는 문제들은 다음과 같다.
- 처음 보는 residue를 만나 cell이 중간에 멈춘 경우, 계산한 수가 모자라게 된다.
- RNA가 잘려서 'RNA'로 인식되지 않은 경우, 그 bead들이 빠진다.
- PDB에서 분자를 자를 때, 일부 bead 만 남은 경우
하지만, 이 검증은 총합만 비교한다는 한계가 있다.
원자 이름이나 순서까지는 확인하지 않아서, 오류끼리 우연히 상쇄되면 통과할 수 있다.
all_fine 값은 다음 과정에서 Hexagon_mol_counting.top 을 쓸지 말지 결정하는 데 쓰인다.
8. Hexagon_mol_counting.top 쓰기
# Writing the mol_counting file which can be copied or imported in the top.
if all_fine:
with open('Hexagon_mol_counting.top', 'w') as f:
for successive_element in successive_elements:
f.write(f'{successive_element[0]}\t{successive_element[1]}\n')
print('The mol_counting.top as been written succesfully.')
else:
print('No ouput has been generated as the amount of atoms does not match.')
이전의 all_fine이 True 일 때만 파일을 쓴다.
원자 수가 안 맞는 토폴로지는 애초에 grompp에서 에러가 나기 때문에 처음부터 만들지 않는다.
파일 내용은 다음과 같다.
successive_elements의 (이름, 개수) 쌍을 한 줄씩, 탭으로 구분해서 쓴다.
GROMACS는 공백과 탭을 구분하지 않아서 탭을 써도 괜찮다.
이 파일에는 [ molecules ] 헤더가 없고 본문만 있다. 헤더는 이후 과정에서 생성되는 hexagon.top에 들어간다.
9. 최종 토폴로지 hexagon.top 쓰기
topol = open('hexagon.top', 'w')
topol.write(f'#include "ITPs/martini_v3.0.0.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0_ffbonded.itp"\n')
topol.write(f'#include "ITPs/RNA.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_phospholipids_v1.itp"\n')
topol.write(f'#include "ITPs/ILs.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0_sterols_v1.0.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_solvents_v1.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_ions_v1.itp"\n')
topol.write(f'[ system ]\n')
topol.write(f'LNP\n')
topol.write(f'\n')
topol.write(f'[ molecules ]\n')
topol.write(f'#include "Hexagon_mol_counting.top"\n')
topol.close()
만들어지는 파일의 구조는 다음과 같다.
#include "ITPs/martini_v3.0.0.itp" ─┐ ① force field 정의
#include "ITPs/martini_v3.0_ffbonded.itp" ─┘
#include "ITPs/RNA.itp" ─┐
#include "ITPs/martini_v3.0.0_phospholipids_v1.itp" │
#include "ITPs/ILs.itp" │ ② 분자 정의 (moleculetype)
#include "ITPs/martini_v3.0_sterols_v1.0.itp" │
#include "ITPs/martini_v3.0.0_solvents_v1.itp" │
#include "ITPs/martini_v3.0.0_ions_v1.itp" ─┘
[ system ] ─┐ ③ 시스템 이름
LNP ─┘
[ molecules ] ─┐ ④ 분자 목록
#include "Hexagon_mol_counting.top" ─┘ (cell 9에서 만든 파일)
- Force field 정의
Force field 정의는 반드시 맨 앞에 와야한다.
| 파일 | 내용 |
| martini_v3.0.0.itp | [ defaults ], [ atomtypes ] (P1, SC4, TQ5 같은 bead 타입), [ nonbond_params ] (bead 쌍마다의 LJ 상호작용) |
| martini_v3.0_ffbonded.itp | b_C1_C1, a_C1_C4_C1 같은 결합과 각도 파라미터 매크로(#define) |
GROMACS는 위에서 아래로 읽기 때문에, bead 타입과 매크로가 분자 정의보다 먼저 정의되어 있어야 한다.
실제로 ILs.itp 에서 ffbonded.itp 의 매크로를 쓰고 있어서, 이 순서를 바꾸면 에러가 난다.
- 분자 정의
cell 3에서 vermouth로 읽은 ITP와 같은 파일들이다.
이 파일 안의 [ moleculetype ] 이름 (MC3H, CHOL, W, CL 등)이 ④에서 쓰는 이름과 정확히 일치해야 한다.
- [ molecules ] + #include
GROMACS의 #include는 C 전처리기처럼 해당 위치에 파일 내용을 그대로 붙여넣는다.
분자 목록을 별도 파일로 분리해두면, 구조가 바뀌었을 때 cell 9만 다시 실행해서 목록만 교체할 수 있다.
여기까지 진행해서 토폴로지 파일을 만들었다..!
이제 다음 노트북으로 넘어가보자....