[LNP] Step7_Solvate
드디어 마지막 노트북이다!
조금만 더 힘을 내보자.
이 노트북은 앞 단계에서 만든 Inverse Hexagonal 코어 LNP 입자(지질 코팅 포함, Martini 3 CG)를 물과 이온이 있는 박스에 넣고, 평형화/production용 GROMACS 입력(.tpr)까지 만드는 과정을 담고있다.
- 전체 흐름
| 단계 | 하는 일 | 주요 파일 |
| 1. 용매화 | insane-ILs.py로 입자 주위를 WN 물로 채움. cubic 박스, 한 변은 입자 지름 + 7 nm | hexagon_coat_solv.gro |
| 2. 인덱스 생성 | genion 에서 쓸 WN 그룹 정의 | wat.ndx |
| 3. 토폴로지 (수동) | 별도 노트북에서 .top을 만들고, WN 줄을 main top의 [ molecules ] 로 옮김 | hexagon_coat_solv.top |
| 4. 중화 | genion -neutral 로 WN을 Na+/Cl-로 치환 | memion.gro |
| 5. NaCl 추가 | 150 mM을 목표로 이온 추가 | memion_2.gro |
| 6. 에너지 최소화 준비 | min.mdp 로 m.tpr 생성 | m.tpr |
| 7. 인덱스 재생성 | Particle / Solvent 그룹 분리 (온도 커플링 등에 사용) | index.ndx |
| 8. 평형화/production 준비 | 아래 표의 순서로 grompp 실행 | r1 ~ r5.tpr / p.tpr |
- 평형화 프로토콜 (310 K)
Timestep을 1 fs 에서 20 fs 까지 단계적으로 늘린다.
| 단계 | dt | 시간 | 비고 |
| min | - | 최대 1000 steps | steepest descent |
| rel1 | 1 fs | 10 ps | -DFLEXIBLE |
| rel2 | 2 fs | 100 ps | |
| rel 3 | 4 fs | 200 ps | |
| rel4 | 10 fs | 500 ps | |
| rel5 | 15 fs | 750 ps | |
| prod | 20 fs | 100 ns |
# 라이브러리 불러오기와 GROMACS 연결
import numpy as np
import os
import shutil
import subprocess
import MDAnalysis as md
# lmod 대신 직접 설치한 GROMACS 사용
GMX_BIN = '/NHNHOME/WORKSPACE/26msit002_B/PCH/software/gromacs-2026.3-install/bin'
os.environ['PATH'] = GMX_BIN + os.pathsep + os.environ['PATH']
assert shutil.which('gmx'), 'gmx를 찾을 수 없습니다'- os 는 insane 실행(os.system)과 환경변수 PATH 수정(os.environ)에 쓴다.
- shutil 은 gmx 실행 파일이 있는지 확인 (shutil.which) 하는 데 쓴다.
- subprocess 는 파이썬에서 gmx grompp, gmx genion 을 실행한다.
- MDAnalysis (md) 는 .gro 구조를 읽어 물 비드 수를 세고, 원자를 그룹으로 골라 .ndx 인덱스 파일을 쓴다.
- 외부 파일 경로 지정
insane = "insane-ILs.py"
MDP_loc = 'MDPs'
(수정본)
PYTHON = '/NHNHOME/WORKSPACE/26msit002_B/PCH/env/lnp/bin/python'
WORKDIR = '/NHNHOME/WORKSPACE/26msit002_B/PCH/LNP/Case_studies/Lipid_Nanoparticle_construction/Inverse_Hexagonal_core_LNP/Building_notebooks'
os.chdir(WORKDIR)
insane = "insane-ILs-py3.py"
MDP_loc = "MDPs"
경로를 변수로 만들어 두고, 이후 셀들에서 f-string 으로 끼워 넣는다.
- insane : 입자 주위에 물을 채우는 스크립트
- MDP_loc : GROMACS 시뮬레이션 설정 파일(.mdp) 폴더
Step 1: Solvate system
Step 1은 진공 상태의 LNP 입자를 시뮬레이션 박스 한가운데 놓고, 남은 빈 공간을 물로 채우는 과정이다.
입자가 실제 환경인 수용액에서 움직이도록 만드는 단계이다.
- 왜 Solvation을 시켜야 할까?
- LNP 구조는 물이 있어야 유지된다. 지질은 소수성 꼬리를 물에서 숨기고 친수성 머리를 물 쪽으로 향하면서 모인다. 물이 없으면 이 힘이 없어서, 표면 지질층의 구조와 거동이 현실과 달라진다.
- 입자가 자기 자신과 상호작용하지 않게 한다. MD는 박스가 사방으로 반복되는 주기 경계조건을 쓰므로, 박스가 작으면 입자가 옆 박스의 복제본과 닿는다. 그래서 물 층으로 충분히 떨어뜨린다.
- 이후 단계의 전제 조건이다. 이온은 물 비드를 치환해서 넣고(genion), 압력 조절로 실제 밀도를 맞추는 평형화도 용매가 있어야 의미가 있다.
그냥 W가 아니라 WN인 이유는 다음과 같다.
코어 안의 물과 바깥의 물을 구분하기 위해서이다. 입자 내부 RNA 채널 7개에는 이미 W 529개와 CL 98개씩이 들어 있다(Hexagon_coat_mol_counting.top). WN도 물리적으로는 같은 W 타입 비드지만 이름이 달라서, 이후 genion 이 바깥 물(WN)만 이온으로 치환할 수 있다. 덕분에 코어 조성은 그대로 유지된다.
# 1. insane으로 물 채우기
os.system(f'python2 {insane} -f hexagon_coat_m.gro -o hexagon_coat_solv.gro -d 7 -sol WN -excl 1 -pbc cubic')
(수정본)
ret = os.system(f'{PYTHON} {insane} -f hexagon_coat_m.gro -o hexagon_coat_solv.gro -d 7 -sol WN -excl 1 -pbc cubic')
assert ret == 0, 'insane 실행 실패'
os.system(...)은 문자열을 셀 명령으로 실행한다. 성공하면 0을 돌려주고, 실패해도 예외를 내지 않는다.
insane의 내부 동작은 다음과 같다.
먼저 hexagon_coat_m.gro 파일의 입력 구조를 읽는다. 그 다음, -pdc cubic -d 7 설정으로 인해 정육면체 박스 크기를 입자 지름 + 7 nm로 정한다. 그리고 입자를 박스 중앙으로 옮기면, 격자에서 입자와 겹치는 칸을 제외하고 남은 칸에 WN을 배치하고, 입자 + 물 구조를 hexagon_coat_solv.gro로 저장한다.
# 2. 물 개수 세기와 인덱스 파일 생성
u = md.Universe('hexagon_coat_solv.gro')
no_w = u.select_atoms('resname WN').resids.shape[0]
Waternumber = int( no_w / 4 )
### Create the index file
gro = 'hexagon_coat_solv.gro'
u = md.Universe(gro)
lipidsagg=u.select_atoms('all and not name NA CL W')
solventagg=u.select_atoms('name WN')
lipidsagg.write("wat.ndx", mode="w", name= 'Particle')
solventagg.write("wat.ndx", mode="a", name= 'WN')
u.atoms.write("wat.ndx", mode="a", name= 'System')- md.Universe(...) : .gro 파일을 읽어 원자 이름, 잔기 이름, 좌표를 담은 객체를 만든다.
- select_atoms('resname WN') : 잔기 이름이 WN인 원자만 고른다. 코어의 W는 포함되지 않는ㄴ다.
- .resids.shape[0] : 고른 원자 수다. WN은 비드 1개짜리 분자라서 곧 WN 분자 수이다.
- Waternumber : 실제 물 분자 수가 아니다.
이 코드는 현재 디렉토리 상태와도 맞지 않고 이후에는 python 2를 사용한다.
코드 상의 오류가 있어 수정된 코드로 다시 설명을 이어서 진행하겠다.
일단 먼저, 물을 채워 주는 프로그램인 insane-ILs.py를 지금 서버의 Python(Python 3)에서 실행할 수 있게 문법만 바꿔준다.
먼저 현재 lnp 환경을 활성화시켜주어야 한다. 그리고 아래의 단계들을 터미널에서 진행하면 된다.
# 1. lnp 환경 켜기
conda activate lnp
# 2. 원본을 보존하려고 사본 만들기
cp insane-ILs.py insane-ILs-py3.py
# 3. 사본을 Python 3 문법으로 변환
python -m lib2to3 -w -n insane-ILs-py3.py
/NHNHOME/WORKSPACE/26msit002_B/PCH/env/lnp/bin/python -m lib2to3 -w -n insane-ILs-py3.py
# 4. 잘 변환됐는지 확인 (옵션 없이 문법만 검사)
python -c "compile(open('insane-ILs-py3.py').read(), 'insane', 'exec'); print('OK')"
/NHNHOME/WORKSPACE/26msit002_B/PCH/env/lnp/bin/python -c "compile(open('insane-ILs-py3.py').read(), 'insane', 'exec'); print('OK')"
그리고 노트북의 라이브러리 불러오기 부분의 코드를 다음과 같이 수정해준다.
import numpy as np
import os
import shutil
import subprocess
import MDAnalysis as md
# lmod 대신 직접 설치한 GROMACS 사용
GMX_BIN = '/NHNHOME/WORKSPACE/26msit002_B/PCH/software/gromacs-2026.3-install/bin'
os.environ['PATH'] = GMX_BIN + os.pathsep + os.environ['PATH']
assert shutil.which('gmx'), 'gmx를 찾을 수 없습니다'PYTHON = '/NHNHOME/WORKSPACE/26msit002_B/PCH/env/lnp/bin/python'
WORKDIR = '/NHNHOME/WORKSPACE/26msit002_B/PCH/LNP/Case_studies/Lipid_Nanoparticle_construction/Inverse_Hexagonal_core_LNP/Building_notebooks'
os.chdir(WORKDIR)
insane = "insane-ILs-py3.py"
MDP_loc = "MDPs"
그 다음엔 README 파일에 있는 Step 6.1을 터미널에서 실행해준다.
Step 6.1 ) minimize in vacum the complex:
gmx editconf -f coated_test.pdb -d 7 -c -o coated_test.gro
(echo 'define = -DFLEXIBLE'; cat MDPs/min_vac.mdp) > MDPs/min_vac_flex.mdp
gmx grompp -f MDPs/min_vac_flex.mdp -c coated_test.gro -p hexagon_coat.top -o hexagon_coat_m.tpr
첫 번째 명령어를 실행한 결과는 다음과 같이 나왔다.
Command line:
gmx editconf -f coated_test.pdb -d 7 -c -o coated_test.gro
Note that major changes are planned in future for editconf, to improve usability and utility.
Read 54735 atoms
Volume: 0.001 nm^3, corresponds to roughly 0 electrons
No velocities found
system size : 23.342 21.835 16.065 (nm)
center : -0.028 0.007 0.036 (nm)
box vectors : 0.100 0.100 0.100 (nm)
box angles : 90.00 90.00 90.00 (degrees)
box volume : 0.00 (nm^3)
shift : 18.699 17.910 14.997 (nm)
new center : 18.671 17.918 15.032 (nm)
new box vectors : 37.342 35.835 30.065 (nm)
new box angles : 90.00 90.00 90.00 (degrees)
new box volume :40231.50 (nm^3)
WARNING: No boxtype specified - distance condition applied in each dimension.
If the molecule rotates the actual distance will be smaller. You might want
to use a cubic box instead, or why not try a dodecahedron today?
GROMACS reminds you: "Player Sleeps With the Fishes" (Ein Bekanntes Spiel Von ID Software)
위 사진과 같이 정상적으로 실행된 것을 볼 수 있다.
WARNING: No boxtype specified는 박스 형태를 지정하지 않았다는 안내인데 무시해도 된다.
이제 두 번째 명령어를 실행해준다. 실행 결과는 다음과 같다.
Command line:
gmx grompp -f MDPs/min_vac_flex.mdp -c coated_test.gro -p hexagon_coat.top -o hexagon_coat_m.tpr
Setting the LD random seed to -52695585
Generated 844 of the 356590 non-bonded parameter combinations
...
NOTE 1 [file Hexagon_coat_mol_counting.top, line 50]:
System has non-zero total charge: 1236.000000
Total charge should normally be an integer. See
https://manual.gromacs.org/current/user-guide/floating-point.html
for discussion on how close it should be to an integer.
Cleaning up constraints and constant bonded interactions with virtual sites
Cleaning up constraints and constant bonded interactions with virtual sites
Analysing residue names:
There are: 8803 Other residues
There are: 294 RNA residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Number of degrees of freedom in T-Coupling group rest is 135657.00
The integrator does not provide a ensemble temperature, there is no system ensemble temperature
NOTE 2 [file MDPs/min_vac_flex.mdp]:
You are using a plain Coulomb cut-off, which might produce artifacts.
You might want to consider using PME electrostatics.
This run will generate roughly 4 Mb of data
결과가 정상이고 파일도 min_vac_flex.mdp라서 FLEXIBLE 이 적용된 것을 확인하였다.
자유도도 135657로 괜찮으므로 이제 최소화를 실행한다.
gmx mdrun -deffnm hexagon_coat_m -ntmpi 1 -ntomp 8
실행 결과는 바로 나왔으며 다음과 같다.
Command line:
gmx mdrun -deffnm hexagon_coat_m -ntmpi 1 -ntomp 8
Reading file hexagon_coat_m.tpr, VERSION 2026.3 (single precision)
Update groups can not be used for this system because an incompatible virtual site type is used
1 GPU selected for this run.
Mapping of GPU IDs to the 1 GPU task in the 1 rank on this node:
PP:0
PP tasks will do non-perturbed short-ranged interactions on the GPU
PP task will update and constrain coordinates on the CPU
Using 1 MPI thread
Using 8 OpenMP threads
NOTE: The number of threads is not equal to the number of (logical) cpus
and the -pin option is set to auto: will not pin threads to cpus.
This can lead to significant performance degradation.
Consider using -pin on (and -pinoffset in case you run multiple jobs).
Steepest Descents:
Tolerance (Fmax) = 1.00000e+01
Number of steps = 5000
Energy minimization reached the maximum number of steps before the forces
reached the requested precision Fmax < 10.
writing lowest energy coordinates.
Steepest Descents did not converge to Fmax < 10 in 5001 steps.
Potential Energy = -1.3548814e+06
Maximum force = 2.8585504e+02 on atom 38217
Norm of force = 2.0012729e+01
GROMACS reminds you: "I can't go to a restaurant and order food because I keep looking at the fonts on the menu. Five minutes later I realize that it's also talking about food." (Donald Knuth)
에너지는 −1.35×10⁶ 로 내려갔으며 눌려 있던 지질이 모두 퍼졌고, 분자 사이의 심한 겹침도 사라졌다.
원자 수와 순서는 그대로이고 깨진 좌표도 없다. 비드의 최대 이동 거리는 1.44 nm로 구조가 흩어지지 않았다.
did not converge to Fmax < 10 이라는 안내가 떴지만 최대 힘 286은 진공 최소화의 목적(겹침 해소)에 충분하고, 물을 넣은 뒤 min.mdp로 다시 최소화하기에 괜찮다.
최소화된 최종 구조는 hexagon_coat_m.gro 이다. Pymol로 직접 확인해준다.
이제 Step 7 노트북의 맨 처음 코드부터 Step 1까지 실행해준다.
import numpy as np
import os
import shutil
import subprocess
import MDAnalysis as md
# lmod 대신 직접 설치한 GROMACS 사용
GMX_BIN = '/NHNHOME/WORKSPACE/26msit002_B/PCH/software/gromacs-2026.3-install/bin'
os.environ['PATH'] = GMX_BIN + os.pathsep + os.environ['PATH']
assert shutil.which('gmx'), 'gmx를 찾을 수 없습니다'PYTHON = '/NHNHOME/WORKSPACE/26msit002_B/PCH/env/lnp/bin/python'
WORKDIR = '/NHNHOME/WORKSPACE/26msit002_B/PCH/LNP/Case_studies/Lipid_Nanoparticle_construction/Inverse_Hexagonal_core_LNP/Building_notebooks'
os.chdir(WORKDIR)
insane = "insane-ILs-py3.py"
MDP_loc = "MDPs"ret = os.system(f'{PYTHON} {insane} -f hexagon_coat_m.gro -o hexagon_coat_solv.gro -d 7 -sol WN -excl 1 -pbc cubic')
assert ret == 0, 'insane 실행 실패'u = md.Universe('hexagon_coat_solv.gro')
no_w = u.select_atoms('resname WN').resids.shape[0]
Waternumber = no_w
### Create the index file
gro = 'hexagon_coat_solv.gro'
u = md.Universe(gro)
lipidsagg=u.select_atoms('all and not name NA CL W')
solventagg=u.select_atoms('name WN')
lipidsagg.write("wat.ndx", mode="w", name= 'Particle')
solventagg.write("wat.ndx", mode="a", name= 'WN')
u.atoms.write("wat.ndx", mode="a", name= 'System')
Step 2: Generate new topology in "Step7.5_create_top_hexagon"
Step 2에서는 7.5 notebook을 실행하라고 되어있지만 실제 파일 명은 7.1로 되어있다.
Step 7.1_create_top_hexagon.ipynb 를 돌리고 이어서 진행한다.
2026.11.09 - [무한도전!!/MD 시뮬레이션 돌리자!!] - [LNP] Step7.1_create_top_hexagon
[LNP] Step7.1_create_top_hexagon
Step 7.1 노트북을 실행한 뒤, 다시 Step 7_Solvate 노트북으로 돌아간다. 이 노트북의 목적은 hexagon_coat_solv.gro 를 읽어서 "어떤 분자가 어떤 순서로 몇 개 있는지"를 GROMACS 용 .top 파일로 적는 것이다.
www.cheonggwan.com
이제 두 파일을 수정해준다.
Hexagon_coat_mol_counting.top 에서 WN 275498 줄을 지우고, hexagon_coat_solv.top 의 맨 끝에 붙인다.
Step 3: Neutralize, minimize and equilibrate system
Step 3은 물만 채워진 시스템을 전기적으로 맞추고, 원자 간 충돌을 풀고, 천천히 안정화해서 본 시뮬레이션을 돌릴 수 있는 상태로 만드는 과정이다.
| 순서 | 단계 | 하는 일 | 왜 필요한가 |
| 1 | 중화 | WN 일부를 CL 1236개로 교체 | 총 전하 +1236을 0으로 맞춤. 주기 경계에서 정전기 계산이 왜곡되지 않게 함 |
| 2 | 150 mM NaCl 추가 | WN 일부를 NA, CL 각 2978개로 교체 | 생리적 염 농도를 재현 |
| 3 | 에너지 최소화 (m) | 비드 위치를 미세 조정 | 물과 이온을 넣으며 생긴 무리한 접촉을 풀어야 MD가 터지지 않음. |
| 4 | 인덱스 생성 (index.ndx) | 비드를 Particle과 Solvent 그룹으로 나눔 | 평형화에서 그룹별로 온도 조절 등을 따로 적용 |
| 5 | 평형화 (rel1 ~ rel5) | 짧은 MD를 단계별로 실행 | 온도와 압력(밀도)을 목표값에 서서히 맞춤 |
| 6 | 본 시뮬레이션 입력 (p.tpr) | 평형화된 구조로 입력 파일 생성 | 실제 데이터를 얻을 production MD 준비 |
이온은 모두 WN을 바꿔 넣는 방식이라, 입자 쪽 구조와 코어 안의 채널 물(W)은 건드리지 않는다.
# 1. 중화용 입력 파일 만들기 (grompp)
### Neutralize
subprocess.call(f'gmx grompp -f {MDP_loc}/min.mdp -c hexagon_coat_solv.gro -p hexagon_coat_solv.top -o memion.tpr -maxwarn 1'
, shell=True)
genion 은 어떤 비드에 어떤 전하가 있는지 알아야해서 .tpr 파일이 필요하다.
여기서는 시뮬레이션을 돌리지 않고, genion 에 넘길 파일만 만든다.
- min.mdp : 시뮬레이션 설정 파일
- hexagon_coat_solv.gro : 좌표 (물을 채운 구조)
- hexagon_coat_solv.top : 분자 목록과 힘장
- memion.tpr : 출력 파일
- maxwarn 1 : 경고 1개까지 허용
# 2. CL로 중화 (genion)
#### NB ! Had to include the WN line under moelcules and remove it from the included top file
p = subprocess.Popen("gmx genion -s memion.tpr -o memion.gro -p hexagon_coat_solv.top -pname NA -pq +1 -nname CL -nq -1 -neutral -n wat.ndx"
, stdin=subprocess.PIPE, shell=True, universal_newlines=True)
p.communicate('WN')
p.wait()| 옵션 | 의미 |
| -s memion.tpr | 방금 만든 입력 파일 |
| -o memion.gro | 이온이 들어간 구조 |
| -p hexagon_coat_solv.top | 토폴로지의 분자 수를 직접 수정 |
| -pname NA -pq +1 | 양이온의 이름과 전하 |
| -nname CL -nq -1 | 음이온의 이름과 전하 |
| -neutral | 총 전하가 0이 되는 만큼만 추가 |
| -n wat.ndx | 교체 대상 그룹이 들어 있는 인덱스 |
genion 은 실행 중에 어느 그룹을 이온으로 바꿀지 묻는다.
stdin=subprocess.PIPE 와 p.commuicate('WN') 이 이 질문에 WN 이라고 자동으로 답하고, p.wait() 는 끝날 때까지 기다린다.
총 전하가 +1236 이므로 임의의 WN 1236 개가 CL로 바뀐다.
* 이 셀은 2번 실행하면 안된다.
# 3. NaCl 개수 계산과 입력 파일 다시 만들기
### add 150mM NaCl
naclNUM = int((0.15 * Waternumber*4)/55.5)
subprocess.call(f'gmx grompp -f {MDP_loc}/min.mdp -c memion.gro -p hexagon_coat_solv.top -o memion_2.tpr -maxwarn 4', shell=True)
물은 약 55.5 M 이므로 (이온 수 / 물 분자 수 = 0.15 / 55.5) 라는 비율로 계산한다.
Waternumber는 이전 셀에서 만든 변수라서, 중간에 커널을 재시작했다면 no_w 셀을 다시 실행해야 한다. 안 그러면 NameError가 발생한다.
그 다음 grompp는 중화된 memion.gro 와 수정된 .top 으로 memion_2.tpr을 만든다.
구조와 분자 수가 바뀌었으니 genion 용 입력 파일도 새로 필요하다. 이제 총 전하가 0이라 전하 NOTE는 사라진다.
# 4. NaCL 교체 (genion)
u = md.Universe('memion.gro')
u.select_atoms('name WN').write('wat2.ndx', mode='w', name='WN')
p = subprocess.Popen(f"gmx genion -s memion_2.tpr -o memion_2.gro -p hexagon_coat_solv.top -pname NA -pq +1 -nname CL -nq -1 -np {naclNUM} -nn {naclNUM} -n wat2.ndx",
stdin=subprocess.PIPE, shell=True, universal_newlines=True)
p.communicate('WN')
p.wait()
중화 과정에서 쓴 wat.ndx 는 재사용할 수 없다.
CL 1236개가 들어가면서 WN의 원자 번호가 바뀌었기 때문에, 현재 구조로 번호를 다시 뽑는다.
genion 옵션은 앞과 같고 두 가지만 다르다.
- -np {nac1NUM} -nn {nac1NUM} : NA와 CL을 각각 2978개 추가한다. 같은 수를 넣으니 총 전하는 0으로 유지된다.
- -n wat2.ndx : 새 인덱스를 쓴다.
이 또한 마찬가지로 한 번만 실행해야 한다.
# 5. 최소화 입력 파일 (m.tpr)
subprocess.call(f'gmx grompp -f {MDP_loc}/min.mdp -c memion_2.gro -p hexagon_coat_solv.top -o m.tpr -maxwarn 4', shell=True)
memion_2.gro 와 최종 .top 으로 m.tpr을 만든다.
min.mdp 는 steepest descent로 최대 1000 step을 수행한다.
이 셀은 파일만 만들고 최소화 자체는 하지 않으므로, 터미널에서 직접 돌려야 한다.
아래 코드를 터미널에서 실행해준다.
gmx mdrun -deffnm m -ntmpi 1 -ntomp 16
실행 결과는 다음과 같다.
Steepest Descents did not converge to Fmax < 10 in 1001 steps.
Potential Energy = -8.5537190e+06
Maximum force = 4.2618130e+03 on atom 48253
Norm of force = 1.3183920e+01
최소화가 정상적으로 마무리되었다.
did not converge 가 떴지만 문제는 아니고, 수렴 기준 Fmax < 10 은 매우 엄격한 값이라, min.mdp 의 1000 step 제한에 먼저 도달한 것 뿐이다.
결과를 표로 정리하면 아래와 같다.
| 항목 | 결과 | 판정 |
| 퍼텐셜 에너지 | −8.55×10⁶ kJ/mol | 큰 음수값으로 정상 범위 |
| 최대 힘 (Fmax) | 4262, 48253번 비드 | 조금 높지만 진행 가능 |
| 평균 수준의 힘 (Norm) | 13.2 | 시스템 전체는 잘 풀림 |
| m.gro | 330,233개, NaN 없음, 박스 34.219 nm | 정상 |
| 조성 | WN 268,306 / W 3,703 / NA 2,978 / CL 4,900 | 목표와 일치 |
| .top 끝 부분 | WN 268306, NA 2978, CL 2978, CL 1236 | .gro 와 일치 |
최대 힘이 걸린 48253번 비드는 RNA 10번 A 염기의 SC4 이고, 가장 가까운 다른 잔기는 염기쌍 상대로 보이는 33번 U (0.292 nm) 이다. 이온이나 물과 충돌한 게 아니라 RNA 안쪽의 뻣뻣한 부분이다.
# 6. 인덱스 파일 (index.ndx)
### Create the index file
gro = 'm.gro'
u = md.Universe(gro)
lipidsagg=u.select_atoms('all and not name NA CL W WN')
solventagg=u.select_atoms('name W WN NA CL')
lipidsagg.write("index.ndx", mode="w", name= 'Particle')
solventagg.write("index.ndx", mode="a", name= 'Solvent')
u.atoms.write("index.ndx", mode="a", name= 'System')
rel 과 prod 설정의 tc-grps, comm-grps 가 Solvent 와 Particle 이라는 그룹 이름을 요구해서 만든다.
반드시 터미널로 mdrun을 돌려서 새 m.gro가 생긴 뒤에 실행해야 한다.
# 7. 평형화 입력 파일 (r1 ~ r5)
이제 터미널에서 아래 명령어를 입력하여 평형화를 돌린다.
prev=m
for s in 1 2 3 4 5; do
gmx grompp -f MDPs/rel$s.mdp -c $prev.gro -p hexagon_coat_solv.top -o r$s.tpr -n index.ndx -maxwarn 4 || break
gmx mdrun -deffnm r$s -ntmpi 1 -ntomp 16 || break
prev=r$s
done
r5가 끝나고 결과를 확인한 뒤에 터미널에서 아래 prod 코드를 돌려주면 된다.
gmx grompp -f MDPs/prod.mdp -c r5.gro -p hexagon_coat_solv.top -o p.tpr -n index.ndx -maxwarn 4
nohup gmx mdrun -deffnm p -ntmpi 1 -ntomp 16 > p.out 2>&1 &
2시간 35분 41초 만에 100 ns 계산이 정상 종료되었다.
속도는 924.9 ns/day 이다. trajactory에는 1 ns 간격으로 프레임 101 개가 저장되었다.
| 항목 | 100 ns 평균 | 판정 |
| 온도 | 310.0 K (용매 310.0, 입자 310.0) | 목표와 일치 |
| 압력 | 0.95 ± 6.6 bar | 1 bar에 맞음 |
| 밀도 | 1023.7 kg/m³ | r5 와 같고 변화 없음 |
| 박스 | 33.367 nm | 안정 |
| 에너지 드리프트 | 100 ns 동안 −1326 kJ/mol | 전체(−8.1×10⁶) 대비 무시할 수준 |
prod 결과, 입자는 터지거나 지질이 떨어져 나가지 않았고, 처음의 납작한 모양에서 거의 구형으로 바뀌었다.
(시작과 끝 두 프레임만 비교했을 때) 코어 안쪽은 RNA가 넓게 퍼지고, RNA 주변의 물과 CL은 줄었다. 외부 물도 RNA 근처까지 더 들어왔다. 그러나 전체 프레임의 Rg로 보면 RNA는 퍼지지 않았다.
이제 거의 고지에 도달했다!!!
이제 p.xtc 전체 프레임으로 물 교환, RNA 주변 물/이온 수, 입자 크기(Rg)의 시간 변화를 그려볼 것이다.
위 그래프는 클로드가 분석해서 그려준 그래프들이다.
그래프를 해석해보면 다음과 같다.
- 입자 형태는 지질 4414개가 101개 프레임 모두에서 한 덩어리였고, 초반 20 ns 안에 거의 구형으로 안정됐다.
- 물은 주로 코어 안에서 움직인다. RNA 주변을 벗어난 W는 약 720개인데, 벌크로나간 W의 순증가는 89개뿐이다. 나머지는 입자 안에 머물러 있다. 입자 안이지만 RNA에서 1 nm 넘게 떨어진 W는 1124개에서 1753개로 늘었다.
- 코어 CL도 마찬가지로 원래 코어에 있던 CL 중 RNA 주변 수는 478개에서 331개로 줄었지만, 벌크로 나간 것은 1개뿐이다. 바깥에서 추가한 CL은 RNA까지 거의 오지 못한다.
- 외부 물은 RNA 근처까지 들어왔다. WN은 초반 약 20 ns 동안 늘었다가 340~370개 수준에서 거의 일정하다.
이제 p.gro (100 ns 시점의 최종 구조)를 p_core.pdb로 바꿔서 직접 pymol로 확인해준다.
# 1) 입자 + 코어 물/이온 그룹을 담은 인덱스 생성 (index.ndx는 그대로 둠)
printf 'a 1-54735\nname 3 Core\nq\n' | gmx make_ndx -f p.gro -n index.ndx -o view.ndx
# 2) 주기 경계로 잘린 입자를 한 덩어리로 모으기
printf 'Particle\nSystem\n' | gmx trjconv -s p.tpr -f p.gro -n view.ndx -pbc cluster -o p_tmp.gro
# 3) 지질 + RNA만 (50,346개)
printf 'Particle\nParticle\n' | gmx trjconv -s p.tpr -f p_tmp.gro -n view.ndx -pbc mol -center -conect -o p_particle.pdb
# 4) 지질 + RNA + 코어 W + 코어 CL (54,735개)
printf 'Particle\nCore\n' | gmx trjconv -s p.tpr -f p_tmp.gro -n view.ndx -pbc mol -center -conect -o p_core.pdb
터미널에서 위 명령어들을 실행한 다음 p_core.pdb를 확인해보면 된다.
그 다음 직접 무비로 md 결과를 확인하고 싶으면 아래 코드를 터미널에 입력하면 된다.
printf 'Particle\nCore\n' | gmx trjconv -s p.tpr -f p.xtc -n view.ndx -pbc mol -center -o p_core_movie.xtc
printf 'Particle\nCore\n' | gmx trjconv -s p.tpr -f p.xtc -n view.ndx -pbc mol -center -conect -dump 0 -o p_core_start.pdb
구형의 LNP를 확인할 수 있다.
드디어 끝~~~~!!!!!!