[LNP] Step4_create_top_hexagon

너무 힘들어 ㅠㅠ......

그래도 해야지.....

 

강한 여성은 해낼 것이다....!

 

 

 

## 이전 글 ##

2026.10.30 - [무한도전!!/MD 시뮬레이션 돌리자!!] - [LNP] Step3_hexagon_channels

[LNP] Step3_hexagon_channels

 

 

---

 

이 노트북은 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 객체를 돌려준다.

1. read_itp는 파일 경로가 아니라 줄 리스트를 입력으로 받는다.
2. 분자 정의를 담을 빈 상자이다.
3. 반환값이 없고, 넘겨준 forcefield를 직접 수정한다. ITP의 [ moleculetype ] 섹션마다 Block을 하나씩 만들어 forcefield.blocks 에 등록한다.
4. 다 채워진 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만 다시 실행해서 목록만 교체할 수 있다.

 

 

 

 

 

여기까지 진행해서 토폴로지 파일을 만들었다..!

이제 다음 노트북으로 넘어가보자....