【模拟计算 就找算筹】
以 TPU 单链到五链体系为例｜分享讲稿与核心代码
案例体系：(MDI-BDO-MDI-PTMG₁₂)₃ TPU 单链，左端为 -NCO、右端为 -OH；完成单链参数化后，构建 5 链、80 × 80 × 80 ų 的初始无定形盒，并输出 Amber / LAMMPS 可用拓扑。
分享主线：
参数化化学结构 → 3D 单链 → 小寡聚物电荷 → 全链 GAFF2 参数化→ 多链装箱 → Amber 拓扑 → LAMMPS 数据文件
传统聚合物建模常常依赖手工绘制、逐步点击软件和人工改文件。一旦要改变聚合度、软段长度、端基或链数，前面的工作就要重复。Python 的价值在于把这些变化转换为参数：修改少量变量即可重建整个体系。
本案例中，Python 不替代量化计算或经典力场工具，而是承担三件最重要的工作：
可以把它概括为一句话：把“搭一个模型”变成“定义一套生成模型的规则”。
本案例可调的核心参数如下：
N_PTMG = 12 # 每个软段的 PTMG（丁撑氧）单元数DP = 3 # 聚合度CHG_PTMG = 3 # 电荷用小寡聚物的 PTMG 单元数CHG_DP = 2 # 电荷用小寡聚物的聚合度
只要修改 N_PTMG、DP、CHG_PTMG 或 CHG_DP，就能生成另一种分子量或链结构的 TPU，同时仍沿用后续参数化流程。
2.用参数化 SMILES 自动构建 TPU 单链
首先要解决的问题不是三维坐标，而是“化学结构是否正确”。TPU 由硬段 MDI–BDO–MDI 和软段 PTMG 组成，且链两端具有不同端基。代码以 SMILES 表达键连关系，再由 OpenBabel 补氢、生成三维构象并进行初步几何优化。
下面这段代码用循环生成聚合物主链：每一个重复单元包含两段 MDI 及对应的氨酯连接；n_ptmg 决定软段长度，dp 决定重复次数。
def build_smiles(n_ptmg, dp): """左端自由异氰酸酯，右端 -OH 的 TPU 链。""" ring = [0] closers = [] s = 'O=C=N' # 左端 -N=C=O def mdi(): ring[0] += 1; a = ring[0] ring[0] += 1; b = ring[0] pre = f'c{_lab(a)}ccc(Cc{_lab(b)}ccc(' clo = f'){"cc"+_lab(b)})cc{_lab(a)}' return pre, clo for r in range(dp): pre, clo = mdi(); s += pre; closers.append(clo) s += 'NC(=O)O' + 'CCCC' + 'OC(=O)N' # 氨酯-BDO-氨酯 pre, clo = mdi(); s += pre; closers.append(clo) s += 'NC(=O)O' + 'CCCCO' * n_ptmg # 氨酯 + PTMG if r < dp - 1: s += 'C(=O)N' return s + ''.join(reversed(closers))
生成结构时，脚本还会用 RDKit 复核分子式、分子量、原子数和形式电荷，避免“结构已经生成但化学式不对”的问题。
def make3d_and_write(smiles, prefix, tag=''): m = pybel.readstring('smi', smiles) m.addh() m.make3D(forcefield=FF_OPT, steps=OPT_STEPS) m.localopt(forcefield=FF_OPT, steps=OPT_STEPS) m.write('mol2', f'{prefix}.mol2', overwrite=True) m.write('pdb', f'{prefix}.pdb', overwrite=True)
这里得到的是一个化学上正确、几何上可用的初始构象，不是平衡态聚合物构象。真实的链构象、密度和缠结状态要在后续的最小化、升温和 NPT 平衡中形成。
AM1-BCC 电荷是 GAFF2 常用的电荷方案，但对整条聚合物直接计算会越来越昂贵，也更容易受构象影响。本案例的做法是：构造一个包含端基、内部连接点和软段环境的小寡聚物，先在其上计算 AM1-BCC 电荷，再把电荷按局部化学环境映射回完整链。
这比“相同元素直接复制相同电荷”可靠得多：同样是 O 或 N，位于氨酯、醚键或端基时的电子环境并不一样。
小寡聚物与整链同时生成，但尺寸更小：
smi_full = build_smiles(N_PTMG, DP)make3d_and_write(smi_full, OUT_FULL, '整链')smi_chg = build_smiles(CHG_PTMG, CHG_DP)make3d_and_write(smi_chg, OUT_CHG, '电荷寡聚物')
随后在小寡聚物上运行 AM1-BCC：
antechamber -i tpu_charge.mol2 -fi mol2 \ -o tpu_charge_bcc.mol2 -fo mol2 \ -c bcc -at gaff2 -nc 0 -pf y -dr no
电荷映射的关键是“原子周围若干键范围内的化学环境”。本案例以 3 个键为半径，优先匹配最具体的局部环境；若无法匹配，逐步降低半径，最后才使用元素级兜底。
RADIUS = 3for i in range(fm.GetNumAtoms()): for r in range(RADIUS, 0, -1): k = env_key(fm, i, r) if k is not None and k in ref[r]: vals = ref[r][k] q[i] = sum(vals) / len(vals) used_r[r] += 1 break
映射后还要把整链电荷校正到目标净电荷；本例为中性分子：
NET_CHARGE = 0.0tot = sum(q)corr = (tot - NET_CHARGE) / len(q)q = [x - corr for x in q]
4. GAFF2 参数化：从结构与电荷到可计算拓扑
一个可进行分子动力学模拟的体系，不只是坐标和电荷。还必须有原子类型、键、键角、二面角、非键相互作用参数等完整拓扑。此处 Python 负责前处理和调度，AmberTools 负责 GAFF2 类型识别、参数补全和 Amber 拓扑构建。
先把已映射的电荷写入整链，并赋 GAFF2 原子类型：
printf("hello world!");antechamber -i tpu_full.mol2 -fi mol2 \ -o tpu_gaff.mol2 -fo mol2 \ -at gaff2 -c rc -cf charges.dat -nc 0 -pf y -dr no
然后由 parmchk2 检索 GAFF2 中没有直接覆盖的参数，生成补充参数文件：
parmchk2 -i tpu_gaff.mol2 -f mol2 -o tpu.frcmod -s gaff2
最后用 tleap 合并结构与参数，并写出 Amber 模拟所需的 prmtop 和 inpcrd：
source leaprc.gaff2mol = loadmol2 tpu_gaff.mol2loadamberparams tpu.frcmodcheck molsaveamberparm mol tpu.prmtop tpu.inpcrdsavepdb mol tpu_amber.pdb
parmchk2 能补“缺失的通用参数”，不等于这些参数都已经过目标体系验证。对力学性质、玻璃化温度、气体吸附或界面相互作用等目标，仍需要用实验或更高层级计算进行验证。
5. 从单链到五链无定形体系，并导入 LAMMPS
单链参数化完成后，不应重新为每一条链重复计算电荷和参数。正确做法是将经过验证的单链保存为库，再复制、装箱、重建多链拓扑。这保证每条 TPU 链使用相同的力场定义，也使体系规模扩展更直接。
为了让 tleap 能可靠识别重复链，先把原子名唯一化、残基名统一为 TPU：
INP, OUT, RESNAME = 'tpu_gaff.mol2', 'tpu_uniq.mol2', 'TPU'el = ''.join(ch for ch in c[1] if ch.isalpha())[:1].upper()cnt[el] = cnt.get(el, 0) + 1name = f'{el}{cnt[el]}'
使用 Packmol 放入 5 条链：
tolerance 2.0filetype pdboutput box5.pdbstructure tpu_single.pdb number 5 inside box 1.1.1.79.79.79.end structure
读取装箱后的 PDB，设置周期性盒子，并重建五链 Amber 拓扑：
source leaprc.gaff2loadamberparams tpu.frcmodloadoff tpu.libsys = loadpdb box5.pdbcheck sysset sys box {80.080.080.0}saveamberparm sys tpu5.prmtop tpu5.inpcrd
通过 InterMol 转换至 LAMMPS：
python -m intermol.convert --amb_in tpu5.prmtop tpu5.inpcrd \ --lammps --oname tpu5
最终交付不应只是一个结构图或 PDB 文件，而应是一套可复建的模拟体系：
printf("hello world!");参数化建链脚本 + 电荷映射脚本 + GAFF2 参数文件（mol2 / frcmod） + Amber 拓扑（prmtop / inpcrd） + LAMMPS 数据与输入文件 + 每一步的检查日志
一句话总结：Python 的价值不只是自动生成一条 TPU 链，而是把聚合物的化学结构、局部电荷、力场和多链模拟体系，固化为一条可复用、可检查、可扩展的建模流程。
算筹科技，旨在为高校、科研院所、医院等相关企业提供模拟计算、超算资源、实验测试等服务。我们拥有丰富的计算经验和深厚的案例积累，专业人员一对一全程服务，提供多种定制化的计算方案及算力支持。
