1.首先你需要安装openmmforcefields

 conda install -c conda-forge openmmforcefields

2.然后产生一个小分子模版生成器（a small molecule residue template generator），按照下例

 # 创建一个 openforcefield Molecule 对象（苯，SMILES模式）
from openforcefield.topology import Molecule
molecule = Molecule.from_smiles('c1ccccc1')
# 创建一个SMIRNOFF模版（这个模版会有着最新的Open Force Field参数信息）
from openmmforcefields.generators import SMIRNOFFTemplateGenerator
smirnoff = SMIRNOFFTemplateGenerator(molecules=molecule)
# 创建一个OpenMM力场对象（AMBER ff14SB，水TIP3P）
from simtk.openmm.app import ForceField
forcefield = ForceField('amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml')
# 注册SMIRNOFF模版生成器
forcefield.registerTemplateGenerator(smirnoff.generator)

或者你也可以使用GAFF力场

 # 创建一个 openforcefield Molecule 对象（苯，SMILES模式）
from openforcefield.topology import Molecule
molecule = Molecule.from_smiles('c1ccccc1')
# 创建一个GAFF模版生成器
from openmmforcefields.generators import GAFFTemplateGenerator
gaff = GAFFTemplateGenerator(molecules=molecule)
# 创建一个OpenMM力场对象（AMBER ff14SB，水TIP3P）
from simtk.openmm.app import ForceField
forcefield = ForceField('amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml')
# 注册GAFF模版生成器
forcefield.registerTemplateGenerator(gaff.generator)

# 下面的话会比较麻烦，这是和原文差别很大的一部分
# 我们在导入pdbfile之前，需要先对他进行fix，这是因为你的蛋白格式可能并不是标准格式
# 力场会导入正确的GAFF参数，antechamber将会产生小分子参数
from pdbfixer import PDBFixer
from simtk.openmm.app.pdbfile import PDBFile
import os
def fix_pdb(pdb_id):
    path = os.getcwd()
    if len(pdb_id) != 4:
        print("Creating PDBFixer...")
        fixer = PDBFixer(pdb_id)
        print("Finding missing residues...")
        fixer.findMissingResidues()

        chains = list(fixer.topology.chains())
        keys = fixer.missingResidues.keys()
        for key in list(keys):
            chain = chains[key[0]]
            if key[1] == 0 or key[1] == len(list(chain.residues())):
                print("ok")
                del fixer.missingResidues[key]

        print("Finding nonstandard residues...")
        fixer.findNonstandardResidues()
        print("Replacing nonstandard residues...")
        fixer.replaceNonstandardResidues()
        print("Removing heterogens...")
        fixer.removeHeterogens(keepWater=True)

        print("Finding missing atoms...")
        fixer.findMissingAtoms()
        print("Adding missing atoms...")
        fixer.addMissingAtoms()
        print("Adding missing hydrogens...")
        fixer.addMissingHydrogens(7)
        print("Writing PDB file...")

        PDBFile.writeFile(
            fixer.topology,
            fixer.positions,
            open(os.path.join(path, "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)),
                 "w"),
            keepIds=True)
        return "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)
fix_pdb('4w52_ben.pdb')
from simtk.openmm.app import PDBFile
pdbfile = PDBFile('4w52_ben_fixed_pH_7.pdb')
system = forcefield.createSystem(pdbfile.topology)
from pdbfixer import PDBFixer

from simtk.openmm.app.pdbfile import PDBFile
import os
def fix_pdb(pdb_id):
    path = os.getcwd()
    if len(pdb_id) != 4:
        print("Creating PDBFixer...")
        fixer = PDBFixer(pdb_id)
        print("Finding missing residues...")
        fixer.findMissingResidues()

        chains = list(fixer.topology.chains())
        keys = fixer.missingResidues.keys()
        for key in list(keys):
            chain = chains[key[0]]
            if key[1] == 0 or key[1] == len(list(chain.residues())):
                print("ok")
                del fixer.missingResidues[key]

        print("Finding nonstandard residues...")
        fixer.findNonstandardResidues()
        print("Replacing nonstandard residues...")
        fixer.replaceNonstandardResidues()
        print("Removing heterogens...")
        fixer.removeHeterogens(keepWater=True)

        print("Finding missing atoms...")
        fixer.findMissingAtoms()
        print("Adding missing atoms...")
        fixer.addMissingAtoms()
        print("Adding missing hydrogens...")
        fixer.addMissingHydrogens(7)
        print("Writing PDB file...")

        PDBFile.writeFile(
            fixer.topology,
            fixer.positions,
            open(os.path.join(path, "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)),
                 "w"),
            keepIds=True)
        return "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)
fix_pdb('4w52_ben.pdb')

from simtk.openmm.app import PDBFile
pdbfile = PDBFile('4w52_ben_fixed_pH_7.pdb')
system = forcefield.createSystem(pdbfile.topology)

 

3.使用SystemGenerator来管理你的力场

作为使用模版生成器的替代方案， openmmforcefields包提供了一个SystemGenerator 功能来简化生物聚合物以及小分子力场的管理。为了去使用这些，你可以十分容易的指定你想要使用的小分子力场。

 # 确定一些关键字参数喂给力场
from simtk import unit
from simtk.openmm import app
forcefield_kwargs = { 'constraints' : app.HBonds, 'rigidWater' : True,
'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu }
# Initialize a SystemGenerator using the Open Force Field Initiative 1.
2.0 force field (openff-1.2.0)
# 使用openff-1.2.0起始一个SystemGenerato
from openmmforcefields.generators import SystemGenerator
system_generator = SystemGenerator(forcefields=['amber/ff14SB.xml', 'am
ber/tip3p_standard.xml'], small_molecule_forcefield='openff-1.2.0', for
cefield_kwargs=forcefield_kwargs, cache='db.json')
# 从OpenMM拓扑对象创建一个OpenMM系统以及一系列的openforcefield分子对象
molecules = Molecule.from_file('ben.sdf', file_format='sdf')
# 这里的pdbfile.topology用的是上面的pdbfile
system = system_generator.create_system(pdbfile.topology, molecules=mol
ecules)