基于 GROMACS 和 Python 的分子动力学模拟

1. Conda 创建环境

将以下内容保存为 environment.yml 文件，然后使用 Conda 创建环境：

name: md_simulation
channels:
  - conda-forge
  - defaults
dependencies:
  - python=3.9
  - gromacs=2022=nompi*     # GROMACS 2022，无 MPI 支持（适合单机）
  - matplotlib              # 用于数据可视化
  - nglview                # 用于3D可视化
  - mdtraj                 # 用于轨迹分析
  - notebook               # Jupyter Notebook 支持
  - numpy                  # 科学计算
  - requests               # 网络请求
  - pandas                 # 数据处理
  - seaborn                # 高级可视化
  - compilers              # 编译器（用于 GROMACS 插件）
  - openmpi                # MPI 支持（可选）
  - cmake                  # 用于构建
  - rdkit                  # 配体处理
  - openbabel              # 化学文件格式转换
  - ipywidgets             # Jupyter 小部件支持
  - pip                    # 支持 pip 包
  - pip:
      - gmxapi             # Python API for GROMACS

创建环境

运行以下命令以创建 Conda 环境：

conda env create -f environment.yml
conda activate md_simulation

---

1. 初始化工作目录**

import os
import subprocess
from rdkit import Chem
from rdkit.Chem import AllChem
from openbabel import pybel

# 设置工作目录
WORKDIR = "simulation"
os.makedirs(WORKDIR, exist_ok=True)

---

2. 从 SMILES 转换为 PDB 文件

def smiles_to_pdb(smiles, output_pdb):
    """
    将 SMILES 转换为 PDB 文件
    参数:
        smiles: str，SMILES 字符串
        output_pdb: str，输出 PDB 文件路径
    """
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)  # 添加氢原子
    AllChem.EmbedMolecule(mol, AllChem.ETKDG())  # 使用 ETKDG 算法生成 3D 坐标
    AllChem.UFFOptimizeMolecule(mol)  # 优化分子几何结构
    
    # 使用 OpenBabel 生成标准 PDB 文件
    tmp_pdb = os.path.join(WORKDIR, "tmp.pdb")
    Chem.MolToPDBFile(mol, tmp_pdb)
    ob_conversion = pybel.ob.OBConversion()
    ob_conversion.SetInAndOutFormats("pdb", "pdb")
    ob_mol = pybel.readfile("pdb", tmp_pdb).__next__()
    ob_mol.write("pdb", output_pdb, overwrite=True)
    print(f"已生成配体 PDB 文件：{output_pdb}")

# 测试 SMILES 转换
smiles = "CC(C)(O)CC(=O)NCCn1ccc2ncnc(Nc3ccc(Oc4cccc(c4)C(F)(F)F)c(Cl)c3)c12"
ligand_pdb = os.path.join(WORKDIR, "ligand.pdb")
smiles_to_pdb(smiles, ligand_pdb)

---

3. 生成拓扑文件

def generate_topology(pdb_file, output_dir):
    """
    使用 GROMACS 的 pdb2gmx 工具生成拓扑文件
    参数:
        pdb_file: str，输入 PDB 文件路径
        output_dir: str，输出目录
    """
    os.makedirs(output_dir, exist_ok=True)
    os.chdir(output_dir)
    subprocess.run([
        "gmx", "pdb2gmx",
        "-f", pdb_file,
        "-o", "processed.gro",
        "-p", "topol.top",
        "-i", "posre.itp",
        "-ff", "amber99sb",  # 力场
        "-water", "spce"     # 水模型
    ], check=True)
    print(f"拓扑文件已生成，目录：{output_dir}")
    os.chdir("..")

# 测试生成拓扑
generate_topology(ligand_pdb, os.path.join(WORKDIR, "topology"))

---

4. 合并蛋白质和配体

def combine_protein_ligand(protein_pdb, ligand_gro, output_gro):
    """
    合并蛋白质和配体到同一模拟盒中
    参数:
        protein_pdb: str，蛋白质的 PDB 文件
        ligand_gro: str，配体的 GRO 文件
        output_gro: str，合并后的 GRO 文件
    """
    subprocess.run([
        "gmx", "editconf",
        "-f", protein_pdb,
        "-o", "protein_box.gro",
        "-c", "-d", "1.0", "-bt", "cubic"
    ], check=True)
    subprocess.run([
        "gmx", "solvate",
        "-cp", "protein_box.gro",
        "-cs", "spc216.gro",
        "-o", output_gro,
        "-p", "topol.top"
    ], check=True)
    print(f"已合并蛋白质和配体，生成：{output_gro}")

---

5. 能量最小化

def energy_minimization(system_gro, topol, output_gro):
    """
    对系统进行能量最小化
    参数:
        system_gro: str，系统的 GRO 文件
        topol: str，拓扑文件
        output_gro: str，最小化后的 GRO 文件
    """
    subprocess.run([
        "gmx", "grompp",
        "-f", "minim.mdp",
        "-c", system_gro,
        "-p", topol,
        "-o", "em.tpr"
    ], check=True)
    subprocess.run([
        "gmx", "mdrun",
        "-v", "-deffnm", "em"
    ], check=True)
    print(f"能量最小化完成，输出文件：{output_gro}")

---

6. 分子动力学模拟

def run_md(system_gro, topol, output_trajectory):
    """
    运行分子动力学模拟
    参数:
        system_gro: str，系统的 GRO 文件
        topol: str，拓扑文件
        output_trajectory: str，轨迹文件
    """
    # 准备 NVT 模拟
    subprocess.run([
        "gmx", "grompp",
        "-f", "nvt.mdp",
        "-c", system_gro,
        "-p", topol,
        "-o", "nvt.tpr"
    ], check=True)
    subprocess.run([
        "gmx", "mdrun",
        "-v", "-deffnm", "nvt"
    ], check=True)

    # 准备 NPT 模拟
    subprocess.run([
        "gmx", "grompp",
        "-f", "npt.mdp",
        "-c", "nvt.gro",
        "-p", topol,
        "-o", "npt.tpr"
    ], check=True)
    subprocess.run([
        "gmx", "mdrun",
        "-v", "-deffnm", "npt"
    ], check=True)

    # 生产模拟
    subprocess.run([
        "gmx", "grompp",
        "-f", "md.mdp",
        "-c", "npt.gro",
        "-p", topol,
        "-o", "md.tpr"
    ], check=True)
    subprocess.run([
        "gmx", "mdrun",
        "-v", "-deffnm", "md"
    ], check=True)
    print(f"分子动力学模拟完成，轨迹文件：{output_trajectory}")

---

7. 运行工作流

protein_pdb = os.path.join(WORKDIR, "protein.pdb")  # 替换为实际蛋白质文件路径
ligand_gro = os.path.join(WORKDIR, "topology/processed.gro")
system_gro = os.path.join(WORKDIR, "system.gro")
topol = os.path.join(WORKDIR, "topology/topol.top")
output_trajectory = os.path.join(WORKDIR, "trajectory.xtc")

combine_protein_ligand(protein_pdb, ligand_gro, system_gro)
energy_minimization(system_gro, topol, system_gro)
run_md(system_gro, topol, output_trajectory)

---

注意事项

1. MDP 文件
   需要提供 minim.mdp、nvt.mdp 和 npt.mdp 参数文件。这些文件可以从 GROMACS 示例目录中获取或根据需要自行编写。

2. 测试环境
   在运行代码之前，使用以下命令测试 GROMACS 是否正确安装：

   gmx --version