[LNP] Inverse-Hexagonal: Step1_construct_core.ipynb

④2026.09.29 - [무한도전!!/MD 시뮬레이션 돌리자!!] - [LNP] Step0. Inverse_Hexagonal_core 실행 준비하기

[LNP] Inverse-Hexagonal: Step0_실행 준비하기

 

 

---

 

이제 이어서 Step1 부터 공부해보자.

나는 비전공자라서 코드부터 이해해야 한다.

 

CodeBlock #1 - 라이브러리 불러오기

import numpy as np
import os
import subprocess 
import MDAnalysis as md
from MDAnalysis.analysis.leaflet import LeafletFinder
import re

 

numpy, os, subprocess, MDAnalysis, LeafletFinder가 쓰이는 듯 하다.

re의 경우 왜 있는지 모르겠는데 나의 경우에는 import re 를 삭제하고 실행시켰다.

 

① numpy : 수치 계산
       - 채널 중심 좌표 (centroid) 계산
② os : 환경변수, 셀 명령 실행
       - PATH에 프로그램 경로 추가
       - os.system(...)으로 packmol, gmx 실행
③ subprocess : 외부 프로그램 실행
       - PCG (막구조생성기)로 지질구조 만들기
       - gmx grompp, gmx mdrun으로 에너지 최소화
④ MDAnalysis : 분자 구조파일 읽고 원자를 선택/편집
       - md.Unvierse로 구조 불러오기
       - md.Merge로 RNA와 지질 채널 합치기
⑤ LeafletFinder : leaflet 생성
       - head 그룹 bead를 15Å (거리) 기준으로 묶어 안쪽/바깥쪽 층을 나눔

 

* PCG - 지질 구조 생성

* Packmol - 물과 이온 채우기

* GROMACS (gmx) - 박스 설정, 복제, 에너지 최소화

 

 

 

CodeBlock #2 - 프로그램 찾기

# 설치 위치 확인하기

import os
import shutil

gmx_bin = "/~경로~/software/gromacs-2026.3-install/bin"
packmol_bin = "/~경로~/software/packmol-20.14.3"
ts2cg_bin = "/~경로~/software/TS2CG1.1"
lnp_bin = "/~경로~/env/lnp/bin"

os.environ["PATH"] = os.pathsep.join([
    gmx_bin,
    packmol_bin,
    ts2cg_bin,
    lnp_bin,
    os.environ["PATH"],
])

PCG_path = shutil.which("PCG")
packmol = shutil.which("packmol")

print("gmx:", shutil.which("gmx"))
print("packmol:", packmol)
print("PCG:", PCG_path)

 

PATH를 설정해 설치 위치를 찾는 셀을 사이에 추가로 넣어주었다.

CodeBlock #3가 파일과 프로그램 위치를 변수로 저장한다면 이 셀은 셀이 프로그램을 찾을 수 있게 해준다.(환경설정)

* shutil.which() : PATH에서 실제로 프로그램을 찾아 경로를 돌려준다.

 

os.environ 부분의 코드는 프로그램을 찾아볼 폴더 목록(PATH) 맨 앞에 사용할 소프트웨어 폴더 4개를 끼워 넣는 것이다.

 

터미널에서 gmx라고만 쳐도 프로그램이 실행되는 건, 컴퓨터가 PATH에 적힌 폴더들을 순서대로 뒤져서 gmx라는 파일을 찾기 때문이다. 여기서 PATH는 그냥 폴더 경로들을 : 로 이어 붙인 긴 문자열을 의미한다.

ex) /usr/local/bin + /usr/bin + /bin
=> /usr/local/bin:/usr/bin:/bin

 

GROMACS나 packmol은 직접 설치한 프로그램이라 이 기본 폴더들에 없어서 위치를 알려줘야 한다.

 

 

@ 코드 조각별 의미

코드	의미
os.environ	환경변수들이 담긴 딕셔너리 같은 객체
os.environ["PATH"]	현재 PATH 문자열
os.pathsep	경로 구분자 / (리눅스에서는 ":" 이고 윈도우에서는 ";" 이다)
"구분자".join([a, b, c])	리스트 요소들을 구분자로 이어 붙임 → "a:b:c"

 

* 원래 PATH를 마지막에 같이 넣는 이유 : 

os.environ["PATH"]를 빼고 4개만 넣으면 기존 경로가 지워진다.

그러면 ls, python 같은 기본 명령까지 못 찾게 되므로 "새 폴더 + 기존 PATH" 형태로 덧붙이는 것이다.

 

* 새 폴더를 앞에 두는 이유 : 

PATH는 앞에서부터 찾고 처음 발견한 걸 쓴다.

시스템에 다른 버전의 gmx가 있더라도, 앞에 둔 GROMACS 2026.3이 우선 사용되게 하려는 목적이다.

 

* 이 변경은 현재 노트북 커널(파이썬 프로세스)과 그 안에서 실행하는 명령에만 적용된다. 터미널이나 다른 프로그램의 PATH는 바뀌지 않는다. 커널을 재시작하면 원래대로 돌아가니, 이 셀을 다시 실행해야 한다.

 

 

 

 

CodeBlock #3 - 경로 변수 설정

main_dir = '.'
mdp_loc_BL = 'MDPs/'
ILs_itp = 'MDPs/ILs.itp'
insane = "insane-ILs.py"

PCG_path = '/projects/DAMM/TS2CG1.1-master/PCG' #This is only applicable for TS2CG1.1
packmol = '/projects/DAMM/packmol-20.14.3/packmol'

RNA = 'PDBs/RNA.pdb'
RNA_itp = 'ITPs/RNA.itp'

 

이 셀은 노트북 곳곳에서 쓸 파일과 프로그램 위치를 변수에 저장해두는 설정 셀이다.

계산은 하지 않고, 경로 문자열에 이름만 붙인다. 즉, 경로를 매번 쓰는 대신 변수 이름만 써서 사용하려는 목적이다.

 

해당 셀은 자신의 경로에 맞게 수정해줘야 한다.

PCG_path 줄의 주석을 보면 "TS2CG1.1"을 설치해야 한다고 나와있다.

 

코드블럭 2와 3의 경우 PCG_path와 packmol이 모두 정의되어 코드블럭 2에서 which로 찾은 값은 코드블럭 3의 직접 적은 경로가 덮어쓴다. 두 값이 같은 경로라서 결과에는 차이가 없다.

변수	값	가리키는 것
main_dir	.	현재 폴더( . = 지금 디렉토리)
mdp_loc_BL	MDPs/	GROMACS 시뮬레이션 설정파일(.mdp) 폴더. (min.mdp, rel*.mdp, prod.mdp 등이 있음)
ILs_itp	ITPs/ILs.itp	이온화 지질(MC3H 등)의 분자 정의 파일
insane	insane-ILs.py	막과 용매를 쌓아주는 파이썬 스크립트
PCG_path	PCG 절대경로	TS2CG의 PCG 프로그램 (원통형 지질막 생성)
packmol	packmol 절대경로	분자를 공간에 채워넣는 프로그램
RNA	PDBs/RNA.pdb	Martini 조립자(CG) RNA 구조. (BB, BB1, BB2, SC1, ... bead로 구성)
RNA_itp	ITPs/RNA.itp	RNA의 분자 정의 파일

 

 

 

CodeBlock #4 - 함수 정의

• get_radius ( ) : 지질 채널의 안쪽 반지름 계산

def get_radius (x,y,z):
    # Calculate the centroid of the cylinder (center axis)
    centroid_x = np.mean(x)
    centroid_y = np.mean(y)
    centroid_z = np.mean(z)

    # Calculate the distances from the centroid to each point
    distances = np.sqrt((x - centroid_x)**2 + (y - centroid_y)**2 + (z - centroid_z)**2)

    # Find the minimum distance, which represents the internal radius of the cylinder
    internal_radius = np.min(distances)

    return internal_radius

 

입력 : 원통 모양 지질 채널의 머리 그룹 bead 좌표

 

1 - 모든 점의 평균으로 원통의 중심점을 구한다.

2 - 중심점에서 각 점가지의 거리를 계산한다.

3 - 그중 가장 가까운 거리를 채널의 안쪽 반지름으로 쓴다.

 

---

• write_packmol_input ( ) : packmol 입력파일 작성 (물과 이온 채우기용)

def write_packmol_input(channel_radius, noW=1, noCL=1, channel_length=100):
    f = open('solvate_channel.inp', 'w') 
    f.write(f'tolerance 3.0\nfiletype pdb\n\n')
        
    f.write(f'structure Inner_rna.pdb\n')
    f.write(f'  number 1\n')
    f.write(f'  inside cylinder 0. 0. 0. 0. 0. 1. {(int(channel_radius)*2)+5}. {int(channel_length+10)}.\n')
    f.write(f'end structure\n')

    f.write(f'structure PDBs/water.pdb\n')
    f.write(f'  number {noW}\n')
    f.write(f'  inside cylinder 0. 0. 0. 0. 0. 1. {int(channel_radius)}. {int(channel_length)}.\n')
    f.write(f'end structure\n')

    f.write(f'structure PDBs/CL.pdb\n')
    f.write(f'  number {noCL}\n')
    f.write(f'  inside cylinder 0. 0. 0. 0. 0. 1. {int(channel_radius)}. {int(channel_length)}.\n')
    f.write(f'end structure\n')
            
    f.write('output solvate_channel_rna.pdb\n')
    f.close()
    return

 

 

이 함수는 solvate_channel.inp 파일을 만든다.

 

- tolerance 3.0 : 서로 다른 분자의 bead를 최소 3Å 떨어뜨림

- Inner_rna.pdb (지질 채널+RNA) : 반지름 2R+5, 길이 110

- water.pdb (Martini W bead) / 개수 noW : 반지름 R, 길이 100

- CL.pdb (염화 이온) / 개수 noCL : 반지름 R, 길이 100

 

packmol 문법 "inside cylinder a b c  dx dy dz  r  l"의 의미는 다음과 같다.

- 점 (a, b, c)에서 시작함

- (dx, dy, dz) 방향으로 뻗음

- 반지름은 r, 길이는 l

 

* 이 코드에서는 모두 (0, 0, 0)에서 z축 방향(0, 0, 1)으로 뻗는 원통이다.

지질 채널(+RNA)은 1개짜리 강체로 취급된다. Packmol이 이걸 큰 원통 안에 맞춰 넣고, 물과 이온은 안쪽 작은 원통에만 채운다. 그래서 결과적으로 물이 채널 내부에만 들어가게 된다.

- 'channel_length=100' Å는 PCG 입력의 'Box 10 10 10 (nm)', 즉 원통 길이 10nm와 맞춘 값

- int(channel_radius)는 소수점을 버림 (반지름이 살짝 줄어듦)

- 출력 파일은 solvate_channel_rna.pdb

 

noW, noCL은 다음과 같이 계산됨.

- noW =  지질 수 * 3 - noCL

- noCL = 양전하를 띤 지질머리(N1, NP) 수 - RNA 백본(BB) 수

* ILs.itp에서 MC3H의 NP bead 전하가 +1이고, RNA 뉴클레오타이드는 -1이므로 이 차이만큼 Cl-를 넣어 전체를 중성으로 만든다.

---

• extract_lines ( ) - 텍스트 파일에서 특정 줄 이후만 잘라오기

def extract_lines(input_file, pattern):
    extracted_lines = []
    pattern_found = False

    with open(input_file, 'r') as f:
        for line in f:
            if pattern_found:
                extracted_lines.append(line)

            if pattern in line:
                pattern_found = True

    return extracted_lines

 

파일을 한 줄씩 읽다가 pattern이 들어간 줄을 만나면, 그다음 줄부터 끝까지 모아서 돌려준다.

(pattern이 있는 줄 자체는 포함되지 않는다)

PCG가 만든 system.top에 쓰이는데 즉, 안쪽 leaflet(lower monolayer)의 지질 개수만 가져온다.

 

---

• write_topol ( ) : GROMACS 토폴로지 파일 topol.top 작성

def write_topol (topol, noW, noCL):
    '''Writes out a topol.top file with the solvate channel times 4 including the correct itp files'''
    extracted_lines = extract_lines(topol, 'lower monolayer')
    
    itpDir   = "ITPs"
    IL_itp   = 'ITPs/ILs.itp'
    Sterols  = "ITPs//martini_v3.0_sterols_v1.0.itp"    
    ffbonded = "ITPs/martini_v3.0_ffbonded.itp"
    RNA      = 'ITPs/RNA.itp'

    with open('topol.top', 'w') as topout:
        topout.write(f'#include "{itpDir}/martini_v3.0.0.itp"\n#include "{ffbonded}"\n#include "{RNA}"\n\
#include "{itpDir}/martini_v3.0.0_phospholipids_v1.itp"\
\n#include "{IL_itp}"\n#include "{Sterols}"\n#include "{itpDir}/martini_v3.0.0_solvents_v1.itp"\n\
#include "{itpDir}/martini_v3.0.0_ions_v1.itp"\n')
        topout.write(f'[ system ]\n')
        topout.write('Hexagonal phase\n')
        topout.write(f'\n')
        topout.write(f'[ molecules ]\n')
        for i in range(4):
                topout.write(f'RNA    1\n')
                for line in extracted_lines:
                    topout.write(line)
                topout.write(f'W    {noW}\n')
                topout.write(f'CL    {noCL}\n')
    return

 

#include 부분에서는 Martini 3 force field와 RNA, 인지질, 이온화 지질(ILs), 스테롤, 용매, 이온의 .itp 파일을 불러온다.

경로가 함수 안에 직접 적혀 있어서 ILs_itp, RNA_itp 같은 변수들은 쓰이지 않는다.

 

[ molecules ] 부분에서는 같은 블록을 4번 반복한다.

* GROMACS는 좌표 파일의 분자 순서와 토폴로지의 순서가 정확히 일치해야 한다.

이후에 나오는 'md.Merge(rna, u)' 코드로 인해 한 채널 안의 순서는 RNA → 지질이다. 그리고 packmol이 그 뒤에 W → CL을 붙인다.  이후의 코드 'gmx gencof -nbox 2 2 1'이 이 채널하나를 통째로 2*2=4개를 복제한다.

따라서 "RNA, 지질, 물, 이온" 블록이 4번 연속으로 나오는 구조가 된다.

 

 

 

 

이후에는 다음과 같이 경고문이 적혀져 있다.

## Importat note: Please do not proceed if you are not familiar with Packmol.
### You can learn more about Packmol here: https://m3g.github.io/packmol/userguide.shtml#tutorials
#### Note that you may need to tune tolerances, adjust lipid ratios, and make other parameter modifications

 

Packmol에 익숙하지 않으면 다음 과정을 진행하지 말라는 것이다. 

Packmol에 대해서 공부하고 싶다면 다음 사이트에서 자료를 찾을 수 있다.

 

https://m3g.github.io/packmol/userguide.shtml#tutorials

Packmol - Initial configurations for Molecular Dynamics

 

 

나는 다음 유튜브 영상을 참고해서 공부했다.

https://www.youtube.com/watch?v=U94Nz1IAzUk&t=390s

 

 

그런데 이 튜토리얼로 공부하는 게 더 좋을 듯하다.

https://www.youtube.com/watch?v=VfYsCpkvLLA