[LNP] Step1_construct_core (2)

- 이전 글

2026.10.14 - [무한도전!!/MD 시뮬레이션 돌리자!!] - [LNP] Step1_construct_core (1)

[LNP] Step1_construct_core (1)

 

 

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

## 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

 

 

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

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

 

 

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

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

 

 

 

 

# Step 1: Create initial channel

- Write input file for PCG

O = open('input.str', 'w')

O.write("[Lipids List]\n")	# 섹션 1: 어떤 지질을 얼마나
O.write("Domain 0\n")	# 막 전체를 '도메인 0번' 하나로 취급함
O.write("MC3H 0.67  0.67 0.60\n")	# 형식: 이름 바깥층비율 안쪽층비율 지질당면적(nm^2)
O.write("CHOL 0.33  0.33 0.60\n")
O.write("End\n")
O.write("\n")
O.write("[Shape Data]\n")	# 섹션 2: 어떤 모양으로
O.write("ShapeType Cylinder\n")
O.write("Box 10 10 10\n")
O.write("Density 2\n")
O.write("Thickness 2\n")
O.write("Radius 3\n")
O.write("End\n")
O.close()

 

이 셀은 PCG가 읽을 input.str를 쓰는 셀이다.

내용은 "MC3H와 콜레스테롤을 2:1로 섞어서, 반지름 3 nm * 길이 10 nm짜리 원통형 이중층을 만들어라" 이다.

실제 생성은 PCG를 실행할 때 이루어진다.

 

 

---

- Build the cylinder with PCG

p = subprocess.Popen(f"{PCG_path} -str input.str -Bondlength 0.2 -LLIB Martini3.LIB -defout system -function analytical_shape"
                    , stdin=subprocess.PIPE, shell=True, universal_newlines=True)
p.wait()

 

Python에서 PCG 프로그램을 실행하는 셀이다.

앞에서 작성한 input.str을 PCG에 넘기면, 이중층 원통의 좌표파일 system.gro 와 지질 개수 파일 system.top 이 만들어진다.

인자	의미
shell=True	명령 문자열을 터미널에 친 것처럼 셀이 해석해서 실행
stdin=subprocess.PIPE	PCG에 입력을 보낼 통로를 열어둠 (여기선 아무것도 안 보내서 사실 필요 없음)
universal_newlines=True	입출력을 바이트가 아닌 텍스트로 처리
p.wait()	PCG가 끝날 때까지 기다림. 종료 코드를 돌려주는데 0이면 정상 종료

 

 

PCG 옵션은 다음과 같다. (TS2CG 1.1 소스에서 확인)

옵션	값	의미	기본값
-function	analytical_shape	메시 파일 대신 input.str의 [Shape Data] (원통 수식)로 막을 만들라는 뜻	backmap (PLM으로 만든 메시 폴더 -dts 가 필요)
-str	input.str	앞에서 작성한 설계도 (지질 조성 + 모양)	없음 (필수)
-LLIB	Martini3.LIB	지질 모양 템플릿 라이브러리	없음
-Bondlength	0.2	템플릿 좌표에 곱할 길이 (nm)	0.1
-defout	system	출력파일 이름 앞부분 -> system.gro, system.top	output

 

지정하지 않은 -seed 는 기본값인 9474가 쓰인다. 따라서 몇 번을 실행해도 같은 구조가 나온다.

 

 

---

- Extract the inner leaflet for continuation

u = md.Universe('system.gro')

L = LeafletFinder(u, 'name NP N1 ROH', cutoff=15)

inner_res = ' '.join(f'{i}' for i in L.groups()[1].resids)
inner = u.select_atoms(f'resid {inner_res}')
inner.atoms.write('Inner.pdb')

 

PCG가 만든 이중층 튜브에서 바깥층을 버리고 안쪽 층만 남긴다. 

머리가 원통 중심을 향하는 단층 채널, 즉 역육각상(H_II)의 기본 단위를 얻는 단계다.

 

 

1. PCG 결과를 MDAnalysis로 불러온다.

u = md.Universe('system.gro')

 

- Universe는 MDAnalysis에서 "분자 시스템 전체"를 담는 객체다. 원자 이름, 잔기 이름, 좌표가 모두 들어 있다.

- .gro 파일은 nm 단위지만, MDAnalysis는 내부적으로 Å로 자동 변환한다.

 

 

2. 막의 층(leaflet)을 자동으로 구분하는 도구 사용

LeafletFinder(u, 'name NP N1 ROH', cutoff=15)

 

작동 원리: 머리 bead를 점으로 보고, 15 Å 안에 있는 점끼리 선으로 이어서 서로 연결된 덩어리가 하나의 층이 된다.

- 'name NP N1 ROH' : 층 구분의 기준으로 쓸 머리 bead.

- cutoff=15 : 머리 bead끼리 15 Å 이내면 '이웃'으로 연결

 

 

3. L.groups()[1]이 왜 안쪽층일까?

L.groups()는 찾은 덩어리들을 크기가 큰 순서로 돌려준다.

원통에서는 바깥층 둘레가 항상 더 길어서 지질도 많다. 그래서 [1]이 항상 안쪽층이 된다.

 

 

4. 머리 bead → 지질 분자 전체로 확장

inner_res = ' '.join(f'{i}' for i in L.groups()[1].resids)
inner = u.select_atoms(f'resid {inner_res}')

 

LeafletFinder 결과에는 머리 bead만 들어 있다. 꼬리까지 포함한 지질 전체가 필요하니 잔기 번호(resid)로 다시 선택한다.

 

1) 안쪽층 머리들의 resid를 꺼낸다. (419, 420, 421, ..., 627)

2) 공백으로 이어 붙여 문자열을 만든다. (149 420 421 ... 627)

3) select_atoms('resid 419 120 ...') 로 그 잔기에 속한 모든 bead를 선택한다.

 

이 방식은 지질마다 resid가 겹치지 않아야 제대로 동작한다.

 

 

5. 안쪽 층을 PDB 파일로 저장하기

inner.atoms.write('Inner.pdb')

 

셀을 실행하면 Found no information for attr: 'altLocs' 같은 경고가 뜨는데, .gro에는 없는 PDB 전용 항목을 기본값으로 채웠다는 뜻이라 무시해도 된다.

 

 

 

 

---

- Solvate one channel and add neutralizing ions using packmol

u = md.Universe('Inner.pdb')
sel = u.select_atoms('name N1 NP ROH')
sel_lipid = u.select_atoms('name N1 NP')
R = get_radius(sel.positions[:,0], sel.positions[:,1], sel.positions[:,2]) - 4 #radius of the channel, trying to substract 2 angstorm to avoid clashing

 

 

1. 방금 저장한 안쪽 층을 다시 불러오기

u = md.Universe('Inner.pdb')

 

변수 u가 이중층에서 안쪽층 단층으로 바뀐다. 

'Element information is missing' 경고는 PDB에 원소 기호 열이 비어있다는 뜻인데 CG bead라서 원소가 없는게 정상이니 무시해도 된다.

 

 

2. 용도가 다른 sel

- sel : 모든 머리 (NP + ROH) / 반지름 측정, 지질 총수 noL 용도
- sel_lipid : 전하를 띤 머리만 (NP) / Cl- 개수 계산 용도

 

sel_lipid를 따로 두는 이유는 전하 대문이다.

MC3H의 NP bead는 +1이지만 콜레스테롤의 ROH는 중성이다.

 

 

3. get_radius로 채널 안쪽 반지름을 재고, 여유분 4 Å를 뺀다.

R = get_radius(...) - 4

 

get_radius : 가장 안쪽 머리 bead까지 거리

R : 물과 이온을 채울 원통 반지름

 

 

* 바꿀 때 주의할 점

- input.str의 반지름이나 두께를 바꾸면 두 층 사이 거리가 달라지므로 cutoff 15가 계속 적절한지 확인한다.

  '[len(g) for g in L.groups()]'를 출력해서 딱 2개 그룹이 나오는지 보면 된다.
- 새 지질을 추가하면 그 지질의 머리 bead 이름을 'name NP N1 ROH'에 넣어야 한다.

  빠지면 그 지질은 어느 층에도 안 잡혀서 Inner.pdb에서 사라진다.
- 전하를 띤 지질이 바뀌면 sel_lipid 선택도 그에 맞게 고쳐야 Cl⁻ 개수가 맞다.

 

 

---

- Insert RNA

 

1. RNA 불러와서 원점에 맞추기

rna = md.Universe(RNA)
rna.atoms.positions = rna.atoms.positions - rna.atoms.center_of_mass()

 

질량 중심을 원점으로 맞춘다. 

positions - center_of_mass()는 모든 bead에서 질량중심 좌표를 빼는 연산이다.

 

 

2. 지질 채널도 원점에 맞추고 → RNA + 지질 합치기 → Inner_rna.pdb

u.atoms.positions = u.atoms.positions - u.atoms.center_of_mass()
sys_rna = md.Merge(rna.select_atoms('all'), u.select_atoms('all'))
sys_rna.atoms.write('Inner_rna.pdb')

 

두 Unvierse의 원자들을 하나의 새 Universe로 합친다. 넣는 순서가 그대로 파일 속 순서가 된다.

 

이 순서가 중요한데, 나중에 write_topol 이 [ molecules ] 에 RNA 1 → MC3H 140 → CHOL 69 순서로 쓴다.

좌표 파일 순서와 다르면 GROMACS가 에러를 낸다.

 

 

3. 합친 파일 다시 불러오기

u = md.Universe('Inner_rna.pdb')

 

u 는 이제 RNA + 지질 채널을 가리킨다.

 

 

4. 머리 bead, RNA 백본 bead 선택

sel = u.select_atoms('name N1 NP ROH')
bb = u.select_atoms('name BB')

 

 

5. 물/이온 개수 계산

hydration = 3 #water beads per lipids, hence 12 water molecules per lipid. This is based on the paper with Nadine
noL = sel.resids.shape[0]
noCL = sel_lipid.resids.shape[0] - bb.resids.shape[0]
noW = (noL*hydration) - noCL

 

1) hydration = 3 : 지질 1개당 Martini 물 bead 3개다. Martini W bead 하나가 실제 물 분자 4개를 나타내서 지질당 물 분자 12개에 해당한다.

2) noCL : 전하 중성 맞추기 (GROMACS는 전체 전하가 0이어야 안정적이라서 남는 양전하만큼 이온을 넣는다)

3) noW : 이온만큼 물을 빼기 (CL 개수만큼 물 bead에서 빼는데, 채널 안 bead 총수를 일정하게 유지해 너무 빽빽해지지 않게 하려는 것이다. 

 

 

6. packmol 입력 파일 쓰기

write_packmol_input(R, noW=noW, noCL=noCL)

 

지금까지 구한 값을 넣어서 solvate_channel.inp를 만든다.

packmol이 이 파일을 읽어서 solvate_channel_rna.pdb를 만든다.

 

 

* 주의할 점

- 질량이 추측값이다.

PDB에는 질량 정보가 없어서 MDAnalysis가 bead 이름을 원소 기호로 보고 질량을 추측한다.

ex) BB → 붕소(B) 10.8, SC1 → 황(S) 32.1, NP → 질소(N) 14.0
ex) GLA, R1처럼 원소로 해석이 안 되는 이름 → 0

 

그래서 center_of_mass()가 실제 Martini 질량(bead당 36–72)과는 다른 값으로 계산된다

다만 여기선 중심 맞추기에만 쓰이고, 기하학적 중심과의 차이가 0.3 Å 정도라 결과에는 영향이 거의 없다.

더 깔끔하게 하려면 center_of_geometry()를 쓰면 된다. 실행 중 뜨는 Unknown masses are set to 0.0 경고가 이 내용이다.

 

- 다른 RNA로 바꿀 때 확인할 것

- 긴 축이 z 방향인지 (아니면 회전시켜야 함)
- 반지름이 채널(~18 Å)보다 충분히 작은지
- 음전하 bead 이름이 BB인지 (noCL 계산이 이걸 전제로 함)
- ITPs/RNA.itp도 같은 서열로 바꿨는지

 

 

---

- Run packmol - Only use packmol for solvate and add ions

os.system(f'{packmol} < solvate_channel.inp')

 

packmol은 solvate_channel.inp 대로 지질 채널 + RNA 1개를 놓고, 안쪽에 물과 CL을 겹치지 않게 채운다.

packmol 시드가 기본값 1234567 로 고정이라 매번 같은 결과가 나온다.

 

만약 다음과 같은 메시지가 뜨며 실패한다면,

ERROR: Could not find a solution...

 

hydration 을 줄이기 / tolerance를 낮추기 / 반지름에서 빼는 값(-4)을 줄이기 중 하나를 조정하면 된다.

 

 

---

- Duplicate the channel 4 times

os.system('gmx editconf -f solvate_channel_rna.pdb -box 5.2 5.2 9.3 -o channel_solvated_combined.gro -bt triclinic -c -angles 90 90 60')
os.system('gmx genconf -f channel_solvated_combined.gro -o dublicate.gro -nbox 2 2 1 -dist 1')

 

1. gmx editconf: 채널 1개에 육각형 박스 씌우기

옵션	의미
-f / -o	입력 pdb, 출력 gro
-bt triclinic	직육면체가 아닌 기울어진 평행육면체 박스
-box 5.2 5.2 9.3	박스 세 변의 길이 (nm)
-angles 90 90 60	변 사이 각도. a와 b 사이가 60º
-c	분자를 박스 중앙으로 이동

 

60º 인 이유는 a-b 사이 각을 60º로 두면 xy 평면이 마름모가 되는데, 이 박스를 주기적으로 이어 붙이면 채널들이 육각형으로 배열된다. 그래서 이름이 역육각상(inverse hexagonal)이다. (각 채널 주위에 이웃 6개)

 

 

2. gmx gencof: 2 X 2로 복제

옵션	의미
-nbox 2 2 1	a 방향 2개 X b 방향 2개 X z 방향 1개 = 채널 4개
-dist 1	복제본 사이에 1 nm씩 추가 간격 (x,y,z 모두)

 

 

---

- Write topol.top

write_topol('system.top', noW, noCL)

 

실제로 생성된 파일은 아래와 같다.

#include "ITPs/martini_v3.0.0.itp"              ← Martini 3 기본 force field (bead 종류, 상호작용)
#include "ITPs/martini_v3.0_ffbonded.itp"
#include "ITPs/RNA.itp"                          ← RNA 분자 정의
#include "ITPs/martini_v3.0.0_phospholipids_v1.itp"
#include "ITPs/ILs.itp"                          ← MC3H 정의
#include "ITPs//martini_v3.0_sterols_v1.0.itp"   ← CHOL 정의
#include "ITPs/martini_v3.0.0_solvents_v1.itp"   ← W 정의
#include "ITPs/martini_v3.0.0_ions_v1.itp"       ← CL 정의
[ system ]
Hexagonal phase

[ molecules ]
RNA    1          ┐
     MC3H  140    │ extract_lines('system.top', 'lower monolayer')로 가져온 줄
     CHOL  69     │
W    529          │ 채널 1
CL    98          ┘
RNA    1          ← 채널 2 (같은 블록 반복)
...               ← 총 4번

 

- #include: "이 분자들이 어떻게 생겼는지"(bead, 결합, 전하) 정의한 파일을 불러온다.
- [ system ]: 시스템 이름이다. 아무 문자열이나 써도 괜찮다.
- [ molecules ]: 좌표 파일에 어떤 분자가 어떤 순서로 몇 개 있는지 적는 곳이다. GROMACS는 좌표 파일의 bead를 처음부터 차례로 읽으면서 이 목록과 하나씩 대조한다.

 

 

* 주의할 점

- RNA 1은 두 가닥(A, U) 합쳐서 한 분자다. RNA.itp가 357 bead를 하나의 분자로 정의하고 있어서 개수가 1이다.
- 4라는 숫자가 write_topol 안에 range(4)로 고정돼 있어서 -nbox를 바꾸면(예: 3 3 1 → 9개) 이 숫자도 직접 바꿔야 한다.
- noW와 noCL은 이전에 설정한 값을 그대로 쓰는데, packmol이 그 개수를 정확히 넣었다는 걸 전제로 한다. packmol이 실패한 채로 넘어가면 여기서는 에러가 안 나고 grompp에서 에러가 난다.

 

 

 

---

또 아래와 같은 경고문이 있다.

## Important note: Please do not proceed if you are not familiar with GROMACS.
### You can learn more about GROMACS here: http://www.mdtutorials.com/gmx/
#### Note that you may need to tune some .mdp parameters based on the GROMACS version being used.

 

만약 GROMACS에 익숙하지 않다면 아래 사이트에서 tutorial을 공부할 수 있다.

http://www.mdtutorials.com/gmx/

GROMACS Tutorials

 

 

# Step 2: Minimize System

지금까지 만든 구조에는 원자끼리 겹친 곳이 많아서 에너지 최소화가 필요하다.

이대로 바로 시뮬레이션 하면 힘이 너무 커서 원자들이 튕겨 나가며 시스템이 터진다.

그래서 먼저 원자를 조금씩 움직여 겹침을 푸는 과정이 필요하다.

 

 

1. 설정 + 구조 + 토폴로지를 하나로 묶어 m.tpr 이라는 파일로 만든다.

p = subprocess.Popen(f"gmx grompp -f {mdp_loc_BL}/min.mdp -c dublicate.gro -p topol.top -o m.tpr -maxwarn 1"
                    , stdin=subprocess.PIPE, shell=True, universal_newlines=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)
p.wait()

 

grompp(GROMACS preprocessor)는 세 가지 입력을 검사하고 합쳐서 실행용 바이너리 .tpr을 만든다.

옵션	파일	역할
-f	MDPs/min.mdp	어떻게 계산할지 (알고리즘, 스텝 수, cutoff 등)
-c	dublicate.gro	어디에 있는지 (좌표, 박스)
-p	topol.top	무엇인지 (분자, 결합, 전하)
-o	m.tpr	결과: 이 셋을 합친 실행파일
-maxwarn 1	-	WARNING 1개까지는 무시하고 진행

 

이 단계에서 좌표와 토폴로지가 일치하는지 검사한다.

 

 

min.mdp 주요 설정은 아래와 같다.

설정	값	의미
integrator	steep	최급강하법. 힘이 가리키는 방향으로 조금씩 이동
nsteps	1000	최대 1000 스텝
emtol	- (기본 10)	최대 힘이 10 KJ/mol/nm 아래로 떨어지면 수렴으로 판단
coulombtype	reaction-field	Martini 표준 정전기 계산 방식
epsilon_r	15	Martini 표준 유전율 (전하 상호작용 약회)
rcoulomb, rvdw	1.1 nm	Martini 3 표준 cutoff
pdb	xyz	3방향 주기 경계

 

 

---

2. 에너지 최소화를 진행하여 m.gro 파일을 생성한다.

p = subprocess.Popen(f"gmx mdrun -v -deffnm m -ntomp 5 -ntmpi 2"
                    , stdin=subprocess.PIPE, shell=True, universal_newlines=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)
p.wait()

옵션	의미
-deffnm m	모든 입출력 파일 이름을 m으로 통일. m.tpr 을 읽고 m.gro, m.log, m.edr, m.trr 을 씀
-v	진행 상황을 스텝마다 출력
-ntmpi 2	계산 영역을 2개로 나눠 병렬 처리
-ntomp 5	각 영역당 스레드 5개 (총 10코어 사용)

 

 

출력 파일은 아래와 같다.

파일	내용
m.gro	최소화된 최종 구조 (다음 단계에서 사용)
m.log	에너지 변화와 최종 결과
m.edr	에너지 기록 (gmx energy로 그래프 가능)
m.trr	궤적

 

* 더 수렴시키고 싶으면 min.mdp의 nsteps 를 늘리면 된다.

 

 

---

3. 원자 그룹을 정의해서 index.ndx 파일을 만든다.

### Create the index file
gro = 'm.gro'
u = md.Universe(gro)	# 최소화된 구조 불러오기
lipidsagg=u.select_atoms('resname MC3H LI2H CHOL A U')	# 이 잔기 이름에 해당하는 원자 ("Lipids" 그룹)
solventagg=u.select_atoms('not resname MC3H LI2H CHOL A U')	# 나머지 전부 "Solvent" 그룹
lipidsagg.write("index.ndx", mode="w", name= 'Lipids')	# 첫 그룹은 새 파일로 쓰기
solventagg.write("index.ndx", mode="a", name= 'Solvent')	# 나머지 그룹은 같은 파일에 이어 붙이기
u.atoms.write("index.ndx", mode="a", name= 'System')	# 전체 원자는 "System" 그룹

 

- 인덱스 파일 (.ndx) 이란 : GROMACS에서 원자 그룹에 이름을 붙여 두는 파일이다. 예를 들어 "지질 그룹과 용매 그룹에 온도 조절을 따로 걸어라" 같은 설정을 할 때 쓴다. 내용은 그룹 이름과 원자 번호 목록이다.

 

 

 

 

 

하... 드디어 Step 1 끝냈다!!!!!!!!!!!!!!!!!!!!!!!