高分子建模

SPONGE 1.4 教程。

高分子建模

更新时间
2023/04/27

本文以50个链节的聚乙二醇的分子构建（如下图）为例。

聚乙二醇分子的构建包含下面四种残基：

其中，EG、eg、eG、Eg为任意与你现有力场不重复的命名。这里使用的EG是我自己的命名，表示Ethylene glycol的缩写。E表示有头部的氢，e表示无头部的氢；G表示有尾部的羟基，g表示无尾部的羟基。

最终目标是构建一个分子

mol = Eg + eg * 48 + eG

也即如下图所示：

1. 构建单体EG力场

import Xponge
import Xponge.forcefield.amber.gaff as gaff
"""从乙二醇的smiles获得assignment,以赋予力场"""
assign = Xponge.Get_Assignment_From_Smiles("OCCO")
"""获得对应的原子类型和电荷"""
assign.determine_atom_type("gaff")
"""如果对精度要求高，使用tpacm4电荷
assign.calculate_charge("resp")
如果对精度要求不高，也可以使用tpacm4电荷"""
assign.calculate_charge("tpacm4")
"""将assignment转为ResidueType"""
EG =  assign.to_residuetype("EG")
"""保存文件"""
save_mol2(EG, "EG.mol2")

打开EG.mol2文件可以检查原子类型和电荷等情况

@<TRIPOS>MOLECULE
EG
 10 9 1 0 1
SMALL
USER_CHARGES
@<TRIPOS>ATOM
     1    O   -1.711    0.268   -0.387   oh     1       EG  -0.633105
     2    C   -0.485   -0.387   -0.314   c3     1       EG   0.105505
     3   C1    0.603    0.581    0.097   c3     1       EG   0.105505
     4   O1    1.843   -0.049    0.177   oh     1       EG  -0.633105
     5    H   -2.186    0.250    0.500   ho     1       EG   0.405197
     6   H1   -0.289   -0.876   -1.280   h1     1       EG   0.061202
     7   H2   -0.563   -1.156    0.505   h1     1       EG   0.061202
     8   H3    0.320    0.955    1.112   h1     1       EG   0.061202
     9   H4    0.701    1.428   -0.585   h1     1       EG   0.061202
    10   H5    1.767   -1.015    0.175   ho     1       EG   0.405197
@<TRIPOS>BOND
     1      1      2 1
     2      1      5 1
     3      2      3 1
     4      2      6 1
     5      2      7 1
     6      3      4 1
     7      3      8 1
     8      3      9 1
     9      4     10 1
@<TRIPOS>SUBSTRUCTURE
    1       EG      1 ****               0 ****  ****

从文件中我们可以看见，名为H的氢原子与第一个羟基氧相连，而后面的羟基名字为O1和H5，这两个信息将用于生成其他的残基。

2. 构建残基eg、Eg和eG的力场

"""接着上面构建单体的内容"""
"""构建残基eG"""
eG = EG.deepcopy("eG")
"""Omit_Atoms的第一个参数是要删掉的原子，第二个参数是删掉以后的残基的电荷。
这里删掉了一个氢离子，所以残基电荷应该是-1"""
eG.Omit_Atoms([eG.H], charge=-1)
save_mol2(eG, "EG1.mol2")
"""重复上面过程，构建残基Eg"""
Eg = EG.deepcopy("Eg")
Eg.Omit_Atoms([Eg.O1, Eg.H5], charge=+1)
save_mol2(Eg, "EG2.mol2")
"""重复上面过程，构建残基eg"""
eg = EG.deepcopy("eg")
eg.Omit_Atoms([eg.H, eg.O1, eg.H5], charge=0)
save_mol2(eg, "EG3.mol2")

3. 连接残基

"""接着上面的构建过程"""
"""eg残基的头部原子名是O
eg残基的主链上紧邻头部原子的原子名是C
eg残基的头部原子成键是1.2埃"""
eg.head = "O"
eg.head_next = "C"
eg.head_length = 1.2
"""eg残基的尾部原子名是C1
eg残基的主链上紧邻尾部原子的原子名是C
eg残基的尾部原子成键是1.2埃"""
eg.tail = "C1"
eg.tail_next = "C"
eg.tail_length = 1.2
"""下面设置键角和二面角的连接规则
eg的O、C与其他残基的尾部原子形成109.5度的键角
eg的C、C1其他残基的头部原子与形成109.5度的键角
eg的C1、C、O与其他残基的尾部原子与形成-180度的二面角
eg的O、C、C1与其他残基的头部原子与形成-180度的二面角"""
PI = 3.141592654
eg.head_link_conditions.append({"atoms": ["C", "O"], "parameter": 109.5 / 180 * PI})
eg.tail_link_conditions.append({"atoms": ["C", "C1"], "parameter": 109.5 / 180 * PI})
eg.head_link_conditions.append({"atoms": ["C1", "C", "O"], "parameter": -PI})
eg.tail_link_conditions.append({"atoms": ["O", "C", "C1"], "parameter": -PI})
"""获得分子"""
mol = eg * 50
"""将第一个残基变为Eg
将最后一个残基变为eG
补充上Eg和eG的原子"""
mol.residues[0].set_type("Eg")
mol.residues[-1].set_type("eG")
mol.add_missing_atoms()
"""使用gaff的parmchk2寻找缺失的力场参数"""
gaff.parmchk2_gaff(mol, "eg.frcmod")
"""保存文件"""
save_mol2(mol, "PEG.mol2")
save_sponge_input(mol, "PEG")

4. 观察模型和动力学模拟

使用vmd观察最终生成的PEG.mol2，可见

5. 后记

本例只包含聚乙二醇部分的建模，更多的处理如溶剂化等与普通过程相同。
例子里的各二面角等可根据自己情况调整。
例子中只使用了单体乙二醇进行力场参数的分配，也可以使用聚三乙二醇进行力场参数分配，不同点在于omit_atoms省略的原子不同

如果想构建无限长的分子，可在最后不set_type，而是对mol连接第0个和最后一个残基的头和尾，也即

"""获得分子"""
mol = eg * 50
"""将第一个残基变为Eg
将最后一个残基变为eG
补充上Eg和eG的原子
mol.residues[0].set_type("Eg")
mol.residues[-1].set_type("eG")
mol.add_missing_atoms()"""
mol.add_residue_link(mol.residues[0].C1, mol.residues[-1].O)