# Molecular dynamics 課程時間:4小時 ## 🧪 一、什麼是分子動力學（MD）？ 分子動力學是一種模擬原子與分子在時間尺度下運動的計算方法，用於研究： - 蛋白質折疊、結合、穩定性 - 材料微觀結構 - 化學反應的動力學 ### 核心概念： - **力場（Force Field）/分子力學** :用來描述系統的勢能。 $E_{total}=E_{bond}+E_{angle}+E_{dihedral}+E_{nonbonded}$ #### 拆解項: | 項目 | 方程式 | 說明 | | -------- | -------- | -------- | | 鍵伸縮 | $E_b=k_b(r-r_o)^2$ | 描述共價鍵的彈性 | | 鍵角 | $E_{\theta}(\theta-\theta_0)^2$ | 原子間鍵角變化 | | 二面角 | $E_{\phi}=k_{\phi}[1+cos(n\phi-\delta)]$ | 扭轉能 | | 非鍵作用 | Lennard-Jones + Coulomb | 模擬范德瓦與靜電作用 | CHARMM36m forcefield: http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs #### 二面角 (Torsion angle) 扭轉角是由四個連續原子定義的角度。 它描述了當你固定前兩個鍵後，第四個原子相對於前三個原子形成平面的旋轉角度。 ![3-s2.0-B9780128210628000119-f01-05-9780128210628](https://hackmd.io/_uploads/BkZxaYWl-l.jpg) - 運動遵循牛頓運動方程式。 $m_i\frac{d^2r_i}{dt^2}=F_i=-\nabla_iU$ - 積分法 Verlet $r(t+Δt) = r(t) + r(t) - r(t-Δt) + a(t)Δt²$ - **時間步長（time step）**約為 1–2 fs（飛秒, $10^{-15}$s），模擬總時間從 ps 到 µs 不等。 - **統計系綜（Ensemble）**控制模擬條件（如 NVE、NVT、NPT）。 ## 🖥️ 二、GROMACS 軟體簡介 GROMACS 是一套開源、高效能的分子動力學模擬軟體，廣泛應用於生物分子模擬。 主要特色： - 支援 CPU 與 GPU 加速 - 開源、可自由使用 - 支援多種常見力場（AMBER、CHARMM、OPLS-AA 等） - 提供結構處理、模擬與後續分析工具 ## 🧰 三、GROMACS 模擬流程概覽 ```mermaid graph TD A[準備初始結構] --> B[產生拓撲檔] B --> C[定義模擬盒並溶劑化] C --> D[加鹽離子] D --> E[能量最小化] E --> F[平衡模擬（NVT / NPT）] F --> G[生產模擬（Production MD）] G --> H[結果分析] ``` ## 🔧 四、基本指令與常見檔案格式 常見輸入檔案： - `.pdb` / `.gro`：分子結構檔 - `.top`：拓撲主檔案 - `.mdp`：模擬參數 - `.itp`：拓撲附檔案 - `.tpr`：整合後的模擬二進位檔 常見指令範例： ```bash gmx pdb2gmx -f protein.pdb -o processed.gro -water spce gmx editconf -f processed.gro -o boxed.gro -c -d 1.0 -bt cubic gmx solvate -cp boxed.gro -cs spc216.gro -o solvated.gro -p topol.top gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -o ionized.gro -p topol.top -pname NA -nname CL -neutral ``` Gromacs提供`gmx dump`來檢查這些二進位檔案的內容: ``` gmx dump -s md.tpr | less gmx dump -s md.xtc | less ``` ## 🔬 五、平衡與生產模擬 ### 能量最小化： ```bash gmx grompp -f minim.mdp -c ionized.gro -p topol.top -o em.tpr gmx mdrun -deffnm em ``` ### 溫度與壓力平衡（NVT / NPT）： ```bash gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr gmx mdrun -deffnm npt ``` ### 生產模擬： ```bash gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr gmx mdrun -deffnm md ``` ### 消去週期性邊界(PBC) PBC 是為了解決模擬系統大小有限的問題。 模擬盒子外層的粒子，會有“影子盒子”來模擬無限系統。每當一個分子走出模擬盒子，它會從對面「回來」，彷彿是在無限大的環境中，但是這在用視覺化軟體觀察時會看到分子斷成兩半、順移或分子拖曳的奇怪現象，因此需要在模擬完成後消去PBC來產生合理的軌跡檔案。 PBC的去除沒有固定方法，需要根據分子軌跡檔中出現的現象來逐一處理，建議的流程參考如下: ``` # 當蛋白分子有部分原子在盒子外拖曳的現象，使用whole來維持蛋白質的完整性 gmx trjconv -s topol.tpr -f traj.xtc -o traj_whole.xtc -center -pbc whole # 如果模擬有兩個分子，可以組成cluster維持兩分子在附近的狀態 gmx trjconv -s topol.tpr -f traj_whole.xtc -o traj_cluster.xtc -pbc cluster # 前兩步驟處理完還是看到有一分子在盒子外順移，可以透過nojump把第一幀當作參考，把跳躍的分子拉回來(這一步只有在遇到分子順移時才需要) gmx trjconv -s equlibration.gro -f traj_cluster.xtc -o traj_nojump.xtc -pbc nojump # 只有一個分子可以用以下指令將分子放到盒子中心方便觀察。center指令會移動原子座標，這一步之後不可以再用nojump，因為參考位置已經不同了! gmx trjconv -s topol.tpr -f traj_nojump.xtc -o traj_center.xtc -center -pbc mol -ur compact # 分子會有旋轉及平移，我們可以用fit來消除，方便觀察 gmx trjconv -s topol.tpr -f traj_center.xtc -o traj_fit.xtc -fit rot+trans ``` ## 📊 六、軌跡與結果分析 常用分析工具： - RMSD：`gmx rms` - 分子半徑：`gmx gyrate` - 氫鍵分析：`gmx hbond` - 能量：`gmx energy` - 聚類分析：`gmx cluster` ## 作業 選擇一個PDB檔案，以CHARMM36m力場進行分子模擬。 載入Gromacs module: ``` module load gromacs/2025.03 ``` 下載charmm36m forcefield: http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs 放在要跑分子模擬的資料夾 ``` tar zxvf charmm36-jul2022.ff.tgz ``` 1. 將PDB檔案轉為gro檔案 選擇charmm36m forcefield 選擇TIP3P 選擇NH2+ 選擇COO- ``` gmx pdb2gmx -f protein.pdb -o protein.gro -ter -ignh 1 1 1 0 ``` 2. 設定模擬盒子 建立9x9x9 nm的模擬盒子 ``` gmx editconf -f protein.gro -o box.gro -c -box 9 9 9 -bt cubic ``` 3. 充填TIP3P水分子 ``` gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top ``` 4. 置換離子保持系統電中性及設定生理離子濃度150mM (0.15M NaCl) ``` gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr echo -e "SOL" | gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15 ``` 5. 能量最小化 ``` gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr gmx mdrun -v -nt 8 -deffnm em ``` 6. NVT平衡 ``` gmx grompp -f nvt_300K.mdp -c em.gro -r em.gro -p topol.top -o nvt_300K.tpr gmx mdrun -v -nt 8 -deffnm nvt_300K ``` 7. NPT平衡 ``` gmx grompp -f npt_300K.mdp -c nvt_300K.gro -r nvt_300K.gro -t nvt_300K.cpt -p topol.top -o npt_300K.tpr gmx mdrun -v -nt 8 -deffnm npt_300K ``` 8. 生產模擬 ``` gmx grompp -f md.mdp -c npt_300K.gro -t npt_300K.cpt -r npt_300K.gro -p topol.top -o md.tpr gmx mdrun -v -deffnm md -nt 30 ``` 9. 分析 去除PBC ``` echo -e "q"| make_ndx -f md.gro -o md.ndx echo -e "1\n1" |gmx trjconv -s md.tpr -f md.txc -o md_center.xtc -center -n md.ndx echo -e "1" |gmx trjconv -s md.tpr -f md_center.xtc -o md_whole.xtc -pbc whole -n md.ndx echo -e "1" |gmx trjconv -s md.tpr -f md_whole.xtc -o md_nojump.xtc -pbc nojump -n md.ndx echo -e "1" | gmx editconf -f md.gro -o md_noSolvent.gro -n md.ndx echo -e "4\n1" |gmx trjconv -s md.tpr -f md_nojump.xtc -o md_fit.xtc -fit rot+trans -n md.ndx ``` RMSD ``` echo -e "4\n4" | gmx rms -f md_fit.xtc -s md.tpr -o md_rmsd.xvg -tu ns -n md.ndx ``` RMSF ``` echo -e "1"|gmx rmsf -f md_fit.xtc -s md.tpr -o md_rmsf.xvg -n md.ndx -res ``` Rg ``` gmx gyrate -f md_fit.xtc -s md.tpr -o md_rg.xvg -n md.ndx -tu ns ``` SASA ``` gmx sasa -f md_fit.xtc -s md.tpr -o md_sasa.xvg -n md.ndx -tu ns ``` ## References 1. http://www.mdtutorials.com/gmx/