무한도전!!/MD 시뮬레이션 돌리자!!
[LNP] Step7.1_create_top_hexagon
청관(淸觀)
2026. 11. 9. 13:58
Step 7.1 노트북을 실행한 뒤, 다시 Step 7_Solvate 노트북으로 돌아간다.
이 노트북의 목적은 hexagon_coat_solv.gro 를 읽어서 "어떤 분자가 어떤 순서로 몇 개 있는지"를 GROMACS 용 .top 파일로 적는 것이다. 이를 위해 구조 파일 - MDAnalysis 을 읽는 과정과 분자 정의(.itp) - vermouth 를 읽는 과정이 필요하다.
# 라이브러리 불러오기
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- MDAnalysis 는 분자 구조/궤적 분석 라이브러리이다. universe.residues 로 잔기를 파일 순서대로 하나씩 훑으면서, 이름이 MC3H, CHOL, WN 인지 등을 확인한다. 순서가 중요한 이유는 .top의 [ molecules ] 순서가 .gro의 순서와 정확히 같아야 하기 때문이다.
- ForceField 는 분자 정의들을 담는 그릇 역할의 클래스이고, read_itp 는 GROMACS .itp 파일을 해석하는 함수이다. Block은 분자 하나의 정의를 의미하는데, 이 정의로 분자별 비드 수를 알 수 있어서 RNA 인식과 검산이 가능해진다.
# Step 1: Generate topology for the solvated system
Step 1은 물을 채운 시스템의 분자 목록(.top)를 새로 만드는 과정이다.
노트북 이름은 토폴로지 생성이지만, 실제로 하는 일은 .gro 에 있는 분자를 순서대로 세어서 적는 것이다.
GROMACS는 .gro - 각 비드의 좌표, .top - 그 비드들이 어떤 분자인지(.itp 정의 + [ molecules ] 목록) 두 파일을 함께 읽는다.
두 파일은 비드 수와 순서가 정확히 일치해야 한다.
insane이 WN를 추가했는데, Step 6 노트북의 .top 에는 WN이 없어서 grompp가 원자 수 불일치 에러를 낸다.
그래서 물이 포함된 목록이 새로 필요하다.
# 1. .itp 파일 하나를 읽는 함수
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()
forcefield = ForceField(name='custom')
read_itp(lines, forcefield)
return forcefield
.itp 파일 하나를 읽어 분자 정의로 바꾸는 함수들을 만든다.
- read_itp 는 값을 돌려주는 대신 넘겨받은 forcefield 를 직접 채운다.
- 결과에서 forcefield.blocks['CHOL'] 처럼 분자 이름으로 정의를 꺼낼 수 있다.
# 2. .itp 6개를 읽어 all_blocks로 합치기
# 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')
# 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)
RNA, 인지질, 이온화 지질, 스테롤, 용매, 이온 정의를 all_blocks 하나로 합친다.
같은 이름이 두 번 나오면 나중 정의가 앞의 것을 덮어쓰므로 경고를 출력한다.
# 3. 분자별 비드 목록 만들기
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은 비드를 노드로 갖는 그래프라서, .nodes(data=True) 로 비드 번호(key)와 속성(values)을 함께 얻는다.
- 속성 중 잔기명과 비드명만 뽑아서 정리한다.
- 이후 셀에서 이 목록의 길이(비드 수) 만 사용한다. RNA 인식과 원자 수 검산에 쓰인다.
# 4. 구조 파일 불러오기
# Read in the GRO file
#!vmd ../hexagon/whole/whole.pdb
universe = mda.Universe('hexagon_coat_solv.gro')
hexagon_coat_solv.gro 를 불러온다.
- 물을 채운 구조를 읽어 universe에 저장한다.
- #!vmd 줄은 원작성자가 구조를 눈으로 볼 때 쓰던 명령이 주석으로 남은 것으로 무시해도 된다.
# 5. 잔기를 분자 이름으로 바꾸기
# 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)
잔기를 처음부터 하나씩 분자 이름으로 바꿔 matches에 저장한다.
맨 위의 주석의 요지는 "완전히 일반적인 방법은 아니니, 분자마다 잔기명이 고유하면 안전하다"를 의미한다.
변수는 다음과 같다.
- matches : 결과를 담을 리스트
- in_RNA : 지금까지 누적한 RNA 비드 수 (아니면 False)
universe.residues 를 파일 순서대로 돌면서 잔기마다 네 갈래 중 하나로 처리한다.
- 잔기명이 곧 분자 이름인 경우.
MC3H, DSPC, CHOL, W, WN 이 여기에 해당한다. - 이온인 경우
이온은 잔기명이 모두 ION 이라 구분이 안된다. 그래서 잔기의 첫 번째 비드 이름(CL)을 분자 이름으로 쓴다. - RNA인 경우
RNA는 잔기 42개가 이어진 한 분자이다. 잔기마다 기록하면 안 되므로 비드 수를 더해 가다가, 357개가 되는 순간 RNA 하나만 기록한다. 길이만 보는 단순한 방식이지만, 이 시스템에서는 RNA 종류가 하나라 충분하다. - 모르는 잔기인 경우
경고를 출력하고 반복을 즉시 멈춘다. 그러면 matches 가 중간에서 끊겨 검산이 실패하고, 파일도 만들어지지 않는다. 결과는 잔기별 이름이 나열된 긴 리스트이다.
# 6. 연속 구간으로 압축하기
successive_elements = []
previous_match = None
counter = None
for match in matches:
if match == previous_match:
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))
같은 분자가 연속되면 .top 에 적을 (이름, 개수) 형태로 압축하는 코드이다.
끝에서 CHOL은 아직 저장되지 않는다. 구간은 다음 분자가 나올 때 저장 되므로 마지막 구간 뒤에는 다음 분자가 없기 때문에 저장되지 않는다. 그래서 반복문 밖의 마지막 줄이 필요하다.
if previous_match is not None 은 첫 분자일 때 (None, None) 이 저장되는 것을 막는다.
#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}')
분자 1개의 (비드 수 * 개수) 를 모두 더해 .gro의 실제 원자 수와 비교한다.
all_fine은 다음 셀에서 파일을 쓸지 결정한다. False 이면 .top 목록 파일이 만들어지지 않는다.
이 검산은 총 비드 수만 본다.
분자 순서는 .gro 순서대로 읽기 때문에 맞춰지고, 비드 이름 불일치는 이후 grompp가 경고로 알려준다.
# 8. 분자 목록 파일 쓰기
# Writing the mol_counting file which can be copied or imported in the top.
if all_fine:
with open('Hexagon_coat_solv_mol_counting.top', 'w') as f:
#with open('Hexagon_coat_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:
'No ouput has been generated as the amount of atoms does not match.'
검산을 통과했을 때만 분자 목록을 파일로 쓴다.
Hexagon_coat_solv_mol_counting.top 을 작성한다.
# 9. 메인 토폴로지 파일 쓰기
topol = open('hexagon_coat_solv.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_coat_solv_mol_counting.top"\n')
topol.close()
.itp 둘을 #include 하고, 위 목록을 불러오는 메인 파일 작성 단계이다.
hexagon_coat_solv.top 을 만드는 셀로, topol.write 로 고정된 내용을 차례로 적고 topol.close() 로 저장한다.
- martini_v3.0.0.itp : 힘장 기본 규칙([defaults]), 비드 타입, 비드 간 비결합 상호 작용
- martini_v3.0_ffbonded.itp : 결합/각도 파라미터를 b_C1_C1 같은 이름으로 정의한 매크로
- RNA. itp : RNA 분자 정의
- martini_v3.0.0_phospholipids_v1.itp : DSPC 등 인지질 정의
- ILs.itp : MC3H 등 이온화 지질 정의
- martini_v3.0_sterols_v1.0.itp : CHOL 정의
- martini_v3.0.0_solvents_v1.itp : W, WN 등 용매 정의
- martini_v3.0.0_ions_v1.itp : NA, CL 등 이온 정의
순서가 중요하다. 비드 타입과 매크로를 분자 정의들이 가져다 쓰므로, 힘장 파일 두 개가 먼저 나와야 한다.
[ molecules ] 의 순서와 개수가 .gro 와 일치해야 한다.