diff --git a/modules/mol/mm/doc/buildingblock.rst b/modules/mol/mm/doc/buildingblock.rst
index c9dbcfd7d936986f5ed18ffa1a3905e852fe6918..1e9f7bf02b6f08eb7a3ab4f82cf6d2eda24eb74e 100644
--- a/modules/mol/mm/doc/buildingblock.rst
+++ b/modules/mol/mm/doc/buildingblock.rst
@@ -1,7 +1,7 @@
Blocks
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The most basic type of residue description is the BuildingBlock. It contains
information on atom names and their corresponding types, charges and
diff --git a/modules/mol/mm/doc/forcefield.rst b/modules/mol/mm/doc/forcefield.rst
index 9c5754152381c0712e176100ffe0d42a8cbb4f3b..9ddbe77ff15a7543604a5993c0232df551ca3d22 100644
--- a/modules/mol/mm/doc/forcefield.rst
+++ b/modules/mol/mm/doc/forcefield.rst
@@ -1,7 +1,7 @@
Forcefields
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The forcefields are a dump for interactions with their parameters, but also
for atom specific information or residue definitions in the form of a
@@ -36,26 +36,26 @@ Loading the standard forcefields provided by OpenStructure
Reading forcefields
--------------------------------------------------------------------------------
-The :class:`FFReader` object is rather experimental. It has nevertheless been
-thoroughly tested for loading the CHARMM and AMBER forcefields in the
-Gromacs format. The reader is capable of resolving the preprocessor statements
-as they are used in Gromacs.
-
.. class:: FFReader(base_dir)
- :param base_dir: Base path of the reader.
- All loaded files must be defined relative to this base
- path.
-
- :type base_dir: :class:`str`
-
The :class:`FFReader` builds up a :class:`Forcefield`, that gets updated with
every call to the read functions. If the read files contain preprocessor
statements as they are used in Gromacs, they will be applied to all
subsequent lines read in. Parsed preprocessor statements are:
#include, #define, #ifdef, #ifndef, #else and #endif
+ Note that this class is rather experimental. It has nevertheless been
+ thoroughly tested for loading the CHARMM and AMBER forcefields in the
+ Gromacs format. The reader is capable of resolving the preprocessor statements
+ as they are used in Gromacs.
+
+ :param base_dir: Base path of the reader.
+ All loaded files must be defined relative to this base
+ path.
+
+ :type base_dir: :class:`str`
+
.. method:: ReadGromacsForcefield()
Searches and reads the forcefield.itp and atomtypes.atp files
@@ -100,8 +100,6 @@ as they are used in Gromacs.
:returns: The reader internal :class:`Forcefield`
-
-
.. code-block:: python
path = "path_to_gromacs/share/top/charmm27.ff"
@@ -111,9 +109,9 @@ as they are used in Gromacs.
reader.ReadGromacsForcefield()
#we also want to read several residue databases
- reader.Read("aminoacids")
- reader.Read("rna")
- reader.Read("dna")
+ reader.ReadResidueDatabase("aminoacids")
+ reader.ReadResidueDatabase("rna")
+ reader.ReadResidueDatabase("dna")
#ions and water are also nice to have, they're stored in itp files
reader.ReadITP("tip3p")
@@ -138,9 +136,40 @@ as they are used in Gromacs.
ff = new_reader.GetForcefield()
ff.Save("charmm_forcefield.dat")
+Generating forcefields with Antechamber
+--------------------------------------------------------------------------------
+
+The antechamber submodule of mol.mm defines functions to use Antechamber (from
+AmberTools) to automatically generate force field parameters and load the
+results into :class:`~ost.mol.mm.Forcefield` objects.
+
+**Example usage**:
+
+.. code-block:: python
+
+ from ost.mol import mm
+
+ # create parameters for TYR using PDB's component dictionary
+ mm.antechamber.RunAntechamber('RVP', 'components.cif', base_out_dir='ligands')
+
+ # create force field
+ ff = mm.Forcefield()
+ ff = mm.antechamber.AddFromPath(ff, 'ligands/RVP')
+ # since Antechamber cannot deal with ions, you can do it manually
+ ff = mm.antechamber.AddIon(ff, 'CL', 'CL', 35.45, -1.0, 0.4401, 0.4184)
+ # save it
+ ff.Save('ligands/ff.dat')
+
+**Functions**:
+
+.. automodule:: ost.mol.mm.antechamber
+ :members:
+
The Forcefield Class
--------------------------------------------------------------------------------
+.. currentmodule:: ost.mol.mm
+
.. class:: Forcefield
@@ -154,7 +183,7 @@ The Forcefield Class
- .. method:: Load(filename)
+ .. staticmethod:: Load(filename)
reads in binary forcefield file
diff --git a/modules/mol/mm/doc/integrators.rst b/modules/mol/mm/doc/integrators.rst
index ed94d86955d0046806e872480146bf05a38393d3..5c4d3524174b2439d110c14f262af8a5b8d196d9 100644
--- a/modules/mol/mm/doc/integrators.rst
+++ b/modules/mol/mm/doc/integrators.rst
@@ -1,7 +1,7 @@
Integrators
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
.. class:: Integrator
diff --git a/modules/mol/mm/doc/interaction.rst b/modules/mol/mm/doc/interaction.rst
index 95d193e1a9e64b9e312590033bc3214d1b198fc9..1b0a705d9a47bab54d4d1775fb45facc04f956fc 100644
--- a/modules/mol/mm/doc/interaction.rst
+++ b/modules/mol/mm/doc/interaction.rst
@@ -1,7 +1,7 @@
Interactions
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The :class:`Interaction` object is intended to build a basic container that can
be used
diff --git a/modules/mol/mm/doc/molmm.rst b/modules/mol/mm/doc/molmm.rst
index f243b7dfb4869e7123ccb0757479745fc3a85b8b..2de5d956abd2ff0ac0179094be97a2e147ecc02e 100644
--- a/modules/mol/mm/doc/molmm.rst
+++ b/modules/mol/mm/doc/molmm.rst
@@ -1,7 +1,7 @@
The mm Module
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
Introduction
--------------------------------------------------------------------------------
diff --git a/modules/mol/mm/doc/observers.rst b/modules/mol/mm/doc/observers.rst
index 1aecc69fed3ff22253aef5f7a8ddfa8e31df1bbd..265236442ccfc7d9dc7df67458c4075dfc6124ea 100644
--- a/modules/mol/mm/doc/observers.rst
+++ b/modules/mol/mm/doc/observers.rst
@@ -1,7 +1,7 @@
Observers
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
Observers can be registered to a :class:`Simulation` and get called at a
defined interval.
diff --git a/modules/mol/mm/doc/settings.rst b/modules/mol/mm/doc/settings.rst
index 648be8bcc11b04995057433bf2644eedc679a94d..d8b3a13b4a713791172284bb79ed797374f7ef70 100644
--- a/modules/mol/mm/doc/settings.rst
+++ b/modules/mol/mm/doc/settings.rst
@@ -1,7 +1,7 @@
The Molecular Mechanics Settings
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The :class:`Settings` define all parameters to control the buildup of a
:class:`Topology` in the :class:`TopologyCreator` and the final setup
diff --git a/modules/mol/mm/doc/simulation.rst b/modules/mol/mm/doc/simulation.rst
index 217c1cc7df013f3b0c427586e62144c1d681e6b9..9d20aac61fa83325437ac4193e793c5e623e8659 100644
--- a/modules/mol/mm/doc/simulation.rst
+++ b/modules/mol/mm/doc/simulation.rst
@@ -1,7 +1,7 @@
Simulation
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The simulation finally connects a :class:`Topology` with an
:class:`EntityHandle`. While applying minimization or
diff --git a/modules/mol/mm/doc/topology.rst b/modules/mol/mm/doc/topology.rst
index 421be2fe9cf5d95ea3cd89352994f5e67bdcbbe6..b616479960ba581b1af8c39dc7210eb784242091 100644
--- a/modules/mol/mm/doc/topology.rst
+++ b/modules/mol/mm/doc/topology.rst
@@ -1,7 +1,7 @@
Topology
================================================================================
-.. currentmodule:: ost.mol
+.. currentmodule:: ost.mol.mm
The :class:`Topology` object is an abstract representation of a protein
structure or any combination of particles that interact with each other in
diff --git a/modules/mol/mm/pymod/CMakeLists.txt b/modules/mol/mm/pymod/CMakeLists.txt
index 3156429e1c24cf25b6c7685fd6507254a02e0f8f..1f47d275beb512dc0911c0390a7b737984415998 100644
--- a/modules/mol/mm/pymod/CMakeLists.txt
+++ b/modules/mol/mm/pymod/CMakeLists.txt
@@ -14,7 +14,8 @@ set(OST_MOL_MM_PYMOD_SOURCES
)
set(OST_MOL_MM_PYMOD_MODULES
- "__init__.py"
+ __init__.py
+ antechamber.py
)
if (NOT ENABLE_STATIC)
diff --git a/modules/mol/mm/pymod/__init__.py b/modules/mol/mm/pymod/__init__.py
index 5c1138cd0354113ea8d1744912a0cc0868202d2b..a4a7c650e1fe7c9bbb0643bdce4a8d4be6e8e6bd 100644
--- a/modules/mol/mm/pymod/__init__.py
+++ b/modules/mol/mm/pymod/__init__.py
@@ -1,6 +1,26 @@
+#------------------------------------------------------------------------------
+# This file is part of the OpenStructure project <www.openstructure.org>
+#
+# Copyright (C) 2008-2016 by the OpenStructure authors
+#
+# This library is free software; you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation; either version 3.0 of the License, or (at your option)
+# any later version.
+# This library is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this library; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+#------------------------------------------------------------------------------
+
import os.path
from _ost_mol_mm import *
-import ost
+import antechamber
+import ost
def LoadAMBERForcefield():
return Forcefield.Load(os.path.join(ost.GetSharedDataPath(),'forcefields','AMBER03.dat'))
diff --git a/modules/mol/mm/pymod/antechamber.py b/modules/mol/mm/pymod/antechamber.py
new file mode 100644
index 0000000000000000000000000000000000000000..4776b2a96a986ecb82a0fc399a6fa557f3330bdd
--- /dev/null
+++ b/modules/mol/mm/pymod/antechamber.py
@@ -0,0 +1,511 @@
+#------------------------------------------------------------------------------
+# This file is part of the OpenStructure project <www.openstructure.org>
+#
+# Copyright (C) 2008-2016 by the OpenStructure authors
+#
+# This library is free software; you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation; either version 3.0 of the License, or (at your option)
+# any later version.
+# This library is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this library; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+#------------------------------------------------------------------------------
+
+# Functions to use Antechamber (from AmberTools) to automatically generate force
+# field parameters. Allows the execution of Antechamber and the parsing of files
+# generated by it.
+
+import _ost_mol_mm as mm
+import ost
+from ost import settings, mol, geom
+import os, subprocess, math
+
+###############################################################################
+# helper functions
+def _GetInteraction(functype, atoms, params):
+ """Get an mm.Interaction with the given func-type and params for the given
+ atoms (name and types extracted from there)."""
+ interaction = mm.Interaction(functype)
+ interaction.SetNames([a.name for a in atoms])
+ interaction.SetTypes([a.GetStringProp('type') for a in atoms])
+ interaction.SetParam(params)
+ return interaction
+
+def _MakeComponentBuildingBlock(eh, ff_dict):
+ """Take EntityHandle eh (from ParseModifiedPDB) and ff_dict (from
+ ParseAmberForceField) and return BuildingBlock."""
+ # converters: length: A -> nm, angles: deg -> rad, energy: kcal -> kJ
+ dist_scale = 1./10.0
+ angle_scale = math.pi/180.
+ # bond strength in OpenMM needs a factor of 2 compared to Amber
+ bond_k_scale = 418.4*2.
+ angle_k_scale = 4.184
+
+ # get atom typing (dictionary listing all atoms of a certain type)
+ atype_dict = {}
+ for a in eh.atoms:
+ atype = a.GetStringProp('type')
+ if not atype in atype_dict:
+ atype_dict[atype] = [a.handle]
+ else:
+ atype_dict[atype].append(a.handle)
+
+ # set masses in entity handle (charges and types set in ParseModifiedPDB)
+ for atype, mass in ff_dict['MASS']:
+ for a in atype_dict[atype]: a.SetMass(mass)
+
+ # start by adding atoms
+ bb = mm.BuildingBlock()
+ for a in eh.atoms:
+ bb.AddAtom(a.name, a.GetStringProp('type'), a.GetCharge(), a.GetMass())
+
+ # add bonds: first all bonds from the entity, then force-field params
+ bl = eh.GetBondList()
+ for b in bl:
+ a1 = b.GetFirst()
+ a2 = b.GetSecond()
+ bond = mm.Interaction(mm.FuncType.HarmonicBond)
+ bond.SetNames([a1.name, a2.name])
+ bond.SetTypes([a1.GetStringProp('type'), a2.GetStringProp('type')])
+ bb.AddBond(bond)
+ added_bonds = []
+ for atype1, atype2, d0, k in ff_dict['BOND']:
+ for a1 in atype_dict[atype1]:
+ for a2 in atype_dict[atype2]:
+ # check connectivity and uniqueness of bond
+ if not mol.BondExists(a1, a2): continue
+ if [a1, a2] in added_bonds or [a2, a1] in added_bonds: continue
+ added_bonds.append([a1,a2])
+ params = [d0*dist_scale, k*bond_k_scale]
+ bond = _GetInteraction(mm.FuncType.HarmonicBond, [a1, a2], params)
+ bb.AddBond(bond, replace_existing=True)
+
+ # add angles
+ added_angles = []
+ for atype1, atype2, atype3, a0, k in ff_dict['ANGL']:
+ # a2 is the central atom
+ for a2 in atype_dict[atype2]:
+ for a1 in atype_dict[atype1]:
+ if not mol.BondExists(a1, a2): continue
+ for a3 in atype_dict[atype3]:
+ # check connectivity and uniqueness of angle
+ if not mol.BondExists(a2, a3): continue
+ if a1 == a3: continue
+ if [a1, a2, a3] in added_angles or [a3, a2, a1] in added_angles:
+ continue
+ added_angles.append([a1, a2, a3])
+ angle = _GetInteraction(mm.FuncType.HarmonicAngle, [a1, a2, a3],
+ [a0*angle_scale, k*angle_k_scale*2])
+ bb.AddAngle(angle)
+
+ # add dihedrals
+ for atype1, atype2, atype3, atype4, idiv, period, phase, k in ff_dict['DIHE']:
+ # there can be multiple ones for the same set of types!
+ added_dihedrals = []
+ for a1 in atype_dict[atype1]:
+ for a2 in atype_dict[atype2]:
+ if not mol.BondExists(a1, a2): continue
+ for a3 in atype_dict[atype3]:
+ if not mol.BondExists(a2, a3): continue
+ if a1 == a3: continue
+ for a4 in atype_dict[atype4]:
+ # check connectivity and uniqueness of dihedral (can be mirrored)
+ if not mol.BondExists(a3, a4): continue
+ if a2 == a4: continue
+ if [a1, a2, a3, a4] in added_dihedrals or \
+ [a4, a3, a2, a1] in added_dihedrals: continue
+ added_dihedrals.append([a1, a2, a3, a4])
+ dihe = _GetInteraction(mm.FuncType.PeriodicDihedral, [a1, a2, a3, a4],
+ [period, phase*angle_scale, k*angle_k_scale])
+ bb.AddDihedral(dihe)
+
+ # add impropers
+ added_impropers = []
+ for atype1, atype2, atype3, atype4, period, phase, k in ff_dict['IMPR']:
+ # third atom is the central atom in amber force-field
+ for ac in atype_dict[atype3]:
+ for a1 in atype_dict[atype1]:
+ if not mol.BondExists(ac, a1): continue
+ for a2 in atype_dict[atype2]:
+ if not mol.BondExists(ac, a2): continue
+ if a1 == a2: continue
+ for a4 in atype_dict[atype4]:
+ # check connectivity and uniqueness of impr. (same central)
+ if not mol.BondExists(ac, a4): continue
+ if a2 == a4 or a1 == a4: continue
+ if [ac, a1, a2, a4] in added_impropers or \
+ [ac, a1, a4, a2] in added_impropers or \
+ [ac, a2, a1, a4] in added_impropers or \
+ [ac, a2, a4, a1] in added_impropers or \
+ [ac, a4, a1, a2] in added_impropers or \
+ [ac, a4, a2, a1] in added_impropers: continue
+ added_impropers.append([ac, a1, a2, a4])
+ impr = _GetInteraction(mm.FuncType.PeriodicImproper, [a1, a2, ac, a4],
+ [period, phase*angle_scale, k*angle_k_scale])
+ bb.AddImproper(impr)
+ return bb
+
+def _ParseModifiedPDB(filename):
+ """Read mpdb file produced by antechamber and return tuple of:
+ - EntityHandle with connectivity, atom types (property 'type') and charges
+ - Residue name as extracted from the mpdb file
+ A RuntimeError is raised if the file can contains multiple residues.
+ """
+ eh = mol.CreateEntity()
+ rname = ''
+ edi = eh.EditXCS(mol.BUFFERED_EDIT)
+ chain = edi.InsertChain('A')
+ bond_list = []
+ # get all atoms and bonds from file
+ with open(filename, 'r') as in_file:
+ for line in in_file:
+ # atom or connectivity
+ # -> fixed column format assumed for both
+ if line.startswith('ATOM'):
+ aname = line[12:17].strip()
+ # extract res. name and ensure uniqueness
+ if not rname:
+ rname = line[17:20].strip()
+ r = edi.AppendResidue(chain, rname)
+ elif rname != line[17:20].strip():
+ raise RuntimeError("More than one residue in file " + filename +\
+ ". Cannot parse!")
+ # extract and store type and charge
+ charge = float(line[54:66])
+ atype = line[78:80].strip()
+ a = edi.InsertAtom(r, aname, geom.Vec3())
+ a.SetStringProp('type', atype)
+ a.SetCharge(charge)
+ elif line.startswith('CONECT'):
+ ai1 = int(line[6:11])
+ # max. 4 bond partners...
+ for i in range(4):
+ try:
+ j = 11 + 5*i
+ ai2 = int(line[j:j+5])
+ # only unique bonds added to list
+ s = set([ai1, ai2])
+ if not s in bond_list: bond_list.append(s)
+ except:
+ # exception thrown for empty strings or non-integers
+ # -> skip
+ continue
+ # set all bonds in entity
+ for indices in bond_list:
+ indices = list(indices)
+ a1 = r.atoms[indices[0]-1]
+ a2 = r.atoms[indices[1]-1]
+ edi.Connect(a1, a2)
+ # finalize
+ edi.UpdateICS()
+ return eh, rname
+
+def _ParseAmberForceField(filename):
+ """Read frcmod file produced by parmchk2 and return dictionary with all
+ entries for masses, bonds, angles, dihedrals, impropers and non-bonded (LJ)
+ interactions. Stored as key/list-of-value pairs:
+ - 'MASS': [atype, mass]
+ - 'BOND': [atype1, atype2, d0, k]
+ - 'ANGL': [atype1, atype2, atype3, a0, k]
+ - 'DIHE': [atype1, atype2, atype3, atype4, idiv, period, phase, k/idiv]
+ - 'IMPR': [atype1, atype2, atype3, atype4, period, phase, k]
+ - 'NONB': [Rvdw, epsilon]
+ """
+ keywords = ['MASS', 'BOND', 'ANGL', 'DIHE', 'IMPR', 'NONB']
+ with open(filename, 'r') as in_file:
+ ff_dict = {}
+ for line in in_file:
+ # look for keywords
+ keyword = line[:4]
+ if not keyword in keywords: continue
+ # loop until empty line found
+ ff_dict[keyword] = []
+ line = in_file.next()
+ while len(line.strip()) > 0:
+ # fixed column format -> extract entries dep. on current keyword
+ if keyword == 'MASS':
+ atype = line[0:2].strip()
+ s = line[2:].split()
+ mass = float(s[0])
+ ff_dict[keyword].append([atype, mass])
+ elif keyword == 'BOND':
+ atype1 = line[:2].strip()
+ atype2 = line[3:5].strip()
+ s = line[5:].split()
+ k = float(s[0])
+ d0 = float(s[1])
+ ff_dict[keyword].append([atype1, atype2, d0, k])
+ elif keyword == 'ANGL':
+ atype1 = line[:2].strip()
+ atype2 = line[3:5].strip()
+ atype3 = line[6:8].strip()
+ s = line[8:].split()
+ k = float(s[0])
+ a0 = float(s[1])
+ ff_dict[keyword].append([atype1, atype2, atype3, a0, k])
+ elif keyword == 'DIHE':
+ atype1 = line[:2].strip()
+ atype2 = line[3:5].strip()
+ atype3 = line[6:8].strip()
+ atype4 = line[9:11].strip()
+ s = line[11:].split()
+ idiv = float(s[0])
+ k = float(s[1])
+ phase = float(s[2])
+ # negative periods = there is more than one term for that dihedral
+ # -> no need to do anything special about that here...
+ period = abs(float(s[3]))
+ ff_dict[keyword].append([atype1, atype2, atype3, atype4, idiv,
+ period, phase, k/float(idiv)])
+ elif keyword == 'IMPR':
+ atype1 = line[:2].strip()
+ atype2 = line[3:5].strip()
+ atype3 = line[6:8].strip()
+ atype4 = line[9:11].strip()
+ s = line[11:].split()
+ k = float(s[0])
+ phase = float(s[1])
+ period = float(s[2])
+ ff_dict[keyword].append([atype1, atype2, atype3, atype4, period,
+ phase, k])
+ elif keyword == 'NONB':
+ line = line.strip()
+ atype = line[0:2].strip()
+ s = line[2:].split()
+ Rvdw = float(s[0])
+ epsilon = float(s[1])
+ ff_dict[keyword].append([atype, Rvdw, epsilon])
+ # next...
+ line = in_file.next()
+ return ff_dict
+###############################################################################
+
+def RunAntechamber(res_name, filename, format='ccif', amberhome=None,
+ base_out_dir=None):
+ """Run Antechamber to guess force field parameters for a given residue name.
+
+ This requires an installation of AmberTools (tested with AmberTools15) with
+ binaries ``antechamber`` and ``parmchk2``.
+
+ This has the same restrictions as Antechamber itself and we assume the input
+ to be uncharged. Note that Antechamber cannot deal with metal ions and other
+ non-organic elements.
+
+ The results are stored in a separate folder named `res_name` within
+ `base_out_dir` (if given, otherwise the current working directory). The main
+ output files are ``frcmod`` and ``out.mpdb``. The former contains force field
+ parameters and masses. The latter maps atom names to atom types and defines
+ the partial charges. The same output could be obtained as follows:
+
+ .. code-block:: console
+
+ $ antechamber -i <FILENAME> -fi <FORMAT> -bk '<RES_NAME>' -o out.mol2 -fo mol2 -c bcc -pf yes
+ $ parmchk2 -i out.mol2 -f mol2 -o frcmod -a Y
+ $ antechamber -i out.mol2 -fi mol2 -o out.mpdb -fo mpdb -pf yes
+
+ The force field parameters can be manually modified if needed. It can for
+ instance happen that some parameters cannot be identified. Those lines will
+ be marked with a comment "ATTN, need revision".
+
+ :param res_name: Residue name for which we desire force field parameters.
+ :type res_name: :class:`str`
+ :param filename: Path to a file which contains the necessary information for
+ `res_name`. It must include all hydrogens.
+ :type filename: :class:`str`
+ :param format: Format of file given with `filename`. Common formats are 'ccif'
+ for PDB's component dictionary or 'pdb' for a PDB file
+ containing the desired residue with all hydrogens.
+ :type format: :class:`str`
+ :param amberhome: Base path of your AmberTools installation. If not None,
+ we look for ``antechamber`` and ``parmchk2`` within
+ ``AMBERHOME/bin`` additionally to the system's ``PATH``.
+ :type amberhome: :class:`str`
+ :param base_out_dir: Path to a base path, where the output will be stored.
+ If None, the current working directory is used.
+ :type base_out_dir: :class:`str`
+ """
+ # find antechamber binaries
+ if amberhome is None:
+ search_paths = []
+ else:
+ search_paths = [os.path.join(amberhome, 'bin')]
+ try:
+ antechamber = settings.Locate('antechamber', search_paths=search_paths)
+ parmchk2 = settings.Locate('parmchk2', search_paths=search_paths)
+ except settings.FileNotFound as ex:
+ ost.LogError("Failed to find Antechamber binaries. Make sure you have "
+ "AmberTools installed!")
+ raise ex
+
+ # prepare path
+ cwd = os.getcwd()
+ if base_out_dir is None:
+ base_out_dir = cwd
+ out_dir = os.path.abspath(os.path.join(base_out_dir, res_name))
+ if not os.path.exists(out_dir):
+ # note: this creates intermediate folders too
+ try:
+ os.makedirs(out_dir)
+ except Exception as ex:
+ ost.LogError("Failed to create output directory " + out_dir + "!")
+ raise ex
+
+ # execute it
+ os.chdir(out_dir)
+ try:
+ cmds = [antechamber + " -i " + filename + " -fi " + format + " -bk " \
+ + res_name + " -o out.mol2 -fo mol2 -c bcc -pf yes",
+ parmchk2 + " -i out.mol2 -f mol2 -o frcmod -a Y",
+ antechamber + " -i out.mol2 -fi mol2 -o out.mpdb -fo mpdb -pf yes"]
+ all_sout = "Generating force field parameters for " + res_name + "\n"
+ all_serr = ""
+ for cmd in cmds:
+ all_sout += "-"*70 + "\n" + "Stdout of: " + cmd + "\n" + "-"*70 + "\n"
+ all_serr += "-"*70 + "\n" + "Stderr of: " + cmd + "\n" + "-"*70 + "\n"
+ job = subprocess.Popen(cmd.split(" "), stdout=subprocess.PIPE,
+ stderr=subprocess.PIPE)
+ sout, serr = job.communicate()
+ all_sout += sout
+ all_serr += serr
+ if job.returncode != 0:
+ ost.LogError("Unsuccessful execution of " + cmd + ". Return code: "\
+ + str(job.returncode))
+ # write command line outputs
+ with open("console.stdout", "w") as txt_file:
+ txt_file.write(all_sout)
+ with open("console.stderr", "w") as txt_file:
+ txt_file.write(all_serr)
+ except Exception as ex:
+ ost.LogError("Failed to excecute antechamber binaries!")
+ raise ex
+
+ # get back to original path
+ os.chdir(cwd)
+
+ # check result
+ frcmod_filename = os.path.join(out_dir, 'frcmod')
+ mpdb_filename = os.path.join(out_dir, 'out.mpdb')
+ if not os.path.exists(frcmod_filename):
+ raise RuntimeError("Failed to generate frcmod file with Antechamber!")
+ if not os.path.exists(mpdb_filename):
+ raise RuntimeError("Failed to generate out.mpdb file with Antechamber!")
+
+def AddFromFiles(force_field, frcmod_filename, mpdb_filename):
+ """Add data from a frcmod and an mpdb file to a force field.
+
+ This will add a new :class:`~ost.mol.mm.BuildingBlock` to `force_field` for
+ the residue defined in those files (residue name is extracted from the mpdb
+ file which can only contain a single residue). Charges for each atom are
+ extracted from the mpdb file. According to the frcmod file, an
+ :class:`~ost.mol.mm.Interaction` is added for each bond, angle, dihedral and
+ improper. Atom types with masses and non-bonded interactions are added to
+ `force_field` as needed.
+
+ :param force_field: A force field object to which the new parameters are
+ added.
+ :type force_field: :class:`~ost.mol.mm.Forcefield`
+ :param frcmod_filename: Path to ``frcmod`` file as generated by ``parmchk2``.
+ :type frcmod_filename: :class:`str`
+ :param mpdb_filename: Path to mpdb file as generated by ``antechamber``.
+ :type mpdb_filename: :class:`str`
+ :return: The updated force field (same as `force_field`).
+ :rtype: :class:`~ost.mol.mm.Forcefield`
+ """
+ # check files
+ if not os.path.exists(frcmod_filename):
+ raise RuntimeError("Could not find frcmod file: " + frcmod_filename)
+ if not os.path.exists(mpdb_filename):
+ raise RuntimeError("Could not find mpdb file: " + mpdb_filename)
+ # read in files
+ try:
+ eh, res_name = _ParseModifiedPDB(mpdb_filename)
+ except Exception as ex:
+ ost.LogError("Failed to parse mpdb file: " + mpdb_filename)
+ raise ex
+ try:
+ ff_dict = _ParseAmberForceField(frcmod_filename)
+ except Exception as ex:
+ ost.LogError("Failed to parse frcmod file: " + frcmod_filename)
+ raise ex
+ ost.LogInfo("Adding force field for " + res_name)
+ # add atoms to force field
+ for aname, mass in ff_dict['MASS']:
+ force_field.AddMass(aname, mass)
+ # add LJs
+ lj_sigma_scale = 2./10./2**(1./6.) # Rvdw to sigma in nm
+ lj_epsilon_scale = 4.184 # kcal to kJ
+ for aname, Rvdw, epsilon in ff_dict['NONB']:
+ # fix 0,0 (from OpenMM's processAmberForceField.py)
+ if Rvdw == 0 or epsilon == 0:
+ Rvdw, epsilon = 1, 0
+ lj = mm.Interaction(mm.FuncType.LJ)
+ lj.SetTypes([aname])
+ lj.SetParam([Rvdw*lj_sigma_scale, epsilon*lj_epsilon_scale])
+ force_field.AddLJ(lj)
+ # add building block
+ bb = _MakeComponentBuildingBlock(eh, ff_dict)
+ force_field.AddBuildingBlock(res_name, bb)
+
+ return force_field
+
+def AddFromPath(force_field, out_dir):
+ """Add data from a directory created with :meth:`Run` to a force field.
+ See :meth:`AddFromFiles` for details.
+
+ :param force_field: A force field object to which the new parameters are
+ added.
+ :type force_field: :class:`~ost.mol.mm.Forcefield`
+ :param out_dir: Output directory as created with :meth:`Run`. Must contain
+ files ``frcmod`` and ``out.mpdb``.
+ :type out_dir: :class:`str`
+ :return: The updated force field (same as `force_field`).
+ :rtype: :class:`~ost.mol.mm.Forcefield`
+ """
+ frcmod_filename = os.path.join(out_dir, 'frcmod')
+ mpdb_filename = os.path.join(out_dir, 'out.mpdb')
+ return AddFromFiles(force_field, frcmod_filename, mpdb_filename)
+
+def AddIon(force_field, res_name, atom_name, atom_mass, atom_charge, lj_sigma,
+ lj_epsilon):
+ """Add a single atom as an ion to a force field.
+
+ Since Antechamber cannot deal with ions, you can add simple ones easily with
+ this function. This adds a :class:`~ost.mol.mm.BuildingBlock` to `force_field`
+ for the given residue name containing a single atom. The atom will have a type
+ with the same name as the atom name and the given mass, charge and non-bonded
+ (LJ) interaction parameters.
+
+ :param force_field: A force field object to which the ion is added.
+ :type force_field: :class:`~ost.mol.mm.Forcefield`
+ :param res_name: Residue name for the ion to be added.
+ :type res_name: :class:`str`
+ :param atom_name: Atom name which is also used as atom type name.
+ :type atom_name: :class:`str`
+ :param atom_mass: Mass of the atom.
+ :type atom_mass: :class:`float`
+ :param atom_charge: Charge of the atom.
+ :type atom_charge: :class:`float`
+ :param lj_sigma: The sigma parameter for the non-bonded LJ interaction.
+ :type lj_sigma: :class:`float` in nm
+ :param lj_epsilon: The sigma parameter for the non-bonded LJ interaction.
+ :type lj_epsilon: :class:`float` in kJ/mol
+ """
+ # add mass (atom_type = atom_name)
+ force_field.AddMass(atom_name, atom_mass)
+ # add building block
+ bb = mm.BuildingBlock()
+ bb.AddAtom(atom_name, atom_name, atom_charge, atom_mass)
+ force_field.AddBuildingBlock(res_name, bb)
+ # add dummy LJ
+ lj = mm.Interaction(mm.FuncType.LJ)
+ lj.SetTypes([atom_name])
+ lj.SetParam([lj_sigma, lj_epsilon])
+ force_field.AddLJ(lj)
+
+__all__ = ('RunAntechamber', 'AddFromFiles', 'AddFromPath', 'AddIon',)
diff --git a/modules/mol/mm/tests/CMakeLists.txt b/modules/mol/mm/tests/CMakeLists.txt
index 99f2c4575dbc77849567d5d9e08a0497dfc9118d..b3eb4d801cb1e650fa6b89b718a5786f84e19e65 100644
--- a/modules/mol/mm/tests/CMakeLists.txt
+++ b/modules/mol/mm/tests/CMakeLists.txt
@@ -6,6 +6,7 @@ set(OST_MOL_MM_UNIT_TESTS
test_forcefield.cc
test_simulation.cc
tests.cc
+ test_antechamber.py
)
ost_unittest(MODULE mol_mm SOURCES "${OST_MOL_MM_UNIT_TESTS}")
diff --git a/modules/mol/mm/tests/TYR/ff_AA.dat b/modules/mol/mm/tests/TYR/ff_AA.dat
new file mode 100644
index 0000000000000000000000000000000000000000..34dc63afee16df9ae3e00bd21c39c47ccc9f0bd7
Binary files /dev/null and b/modules/mol/mm/tests/TYR/ff_AA.dat differ
diff --git a/modules/mol/mm/tests/TYR/frcmod b/modules/mol/mm/tests/TYR/frcmod
new file mode 100644
index 0000000000000000000000000000000000000000..956433fb0d20f94e566d639a0ce392ded96aae29
--- /dev/null
+++ b/modules/mol/mm/tests/TYR/frcmod
@@ -0,0 +1,104 @@
+Remark line goes here
+MASS
+n3 14.010 0.530
+c3 12.010 0.878
+c 12.010 0.616
+o 16.000 0.434
+ca 12.010 0.360
+oh 16.000 0.465
+hn 1.008 0.161
+h1 1.008 0.135
+hc 1.008 0.135
+ha 1.008 0.135
+ho 1.008 0.135
+
+BOND
+c3-n3 320.60 1.470
+hn-n3 394.10 1.018
+c -c3 328.30 1.508
+c3-c3 303.10 1.535
+c3-h1 335.90 1.093
+c -o 648.00 1.214
+c -oh 466.40 1.306
+c3-ca 323.50 1.513
+c3-hc 337.30 1.092
+ca-ca 478.40 1.387
+ca-ha 344.30 1.087
+ca-oh 386.10 1.362
+ho-oh 369.60 0.974
+
+ANGLE
+c -c3-n3 66.590 111.140
+c3-c3-n3 66.180 110.380
+h1-c3-n3 49.390 109.920
+c3-n3-hn 47.130 109.920
+c3-c -o 68.030 123.110
+c3-c -oh 69.840 112.200
+c3-c3-ca 63.250 112.090
+c3-c3-hc 46.370 110.050
+c -c3-c3 63.790 110.530
+c -c3-h1 47.630 107.660
+c -oh-ho 51.190 107.370
+o -c -oh 77.380 122.880
+c3-c3-h1 46.360 110.070
+c3-ca-ca 63.840 120.630
+ca-c3-hc 46.960 110.150
+ca-ca-ca 67.180 119.970
+ca-ca-ha 48.460 120.010
+ca-ca-oh 69.850 119.940
+ca-oh-ho 48.850 109.470
+hn-n3-hn 41.300 107.130
+hc-c3-hc 39.430 108.350
+
+DIHE
+o -c -c3-n3 6 0.000 180.000 2.000
+oh-c -c3-n3 6 0.000 180.000 2.000
+ca-c3-c3-n3 9 1.400 0.000 3.000
+hc-c3-c3-n3 9 1.400 0.000 3.000
+c3-c -oh-ho 2 4.600 180.000 2.000
+c3-c3-ca-ca 6 0.000 0.000 2.000
+c -c3-n3-hn 6 1.800 0.000 3.000
+c -c3-c3-ca 9 1.400 0.000 3.000
+c -c3-c3-hc 9 1.400 0.000 3.000
+o -c -c3-c3 6 0.000 180.000 2.000
+o -c -c3-h1 1 0.800 0.000 -1.000
+o -c -c3-h1 1 0.080 180.000 3.000
+o -c -oh-ho 1 2.300 180.000 -2.000
+o -c -oh-ho 1 1.900 0.000 1.000
+c3-c3-n3-hn 6 1.800 0.000 3.000
+oh-c -c3-c3 6 0.000 180.000 2.000
+c3-ca-ca-ca 4 14.500 180.000 2.000
+c3-ca-ca-ha 4 14.500 180.000 2.000
+ca-c3-c3-h1 9 1.400 0.000 3.000
+ca-ca-ca-ca 4 14.500 180.000 2.000
+ca-ca-ca-ha 4 14.500 180.000 2.000
+hc-c3-ca-ca 6 0.000 0.000 2.000
+ca-ca-ca-oh 4 14.500 180.000 2.000
+ca-ca-oh-ho 2 1.800 180.000 2.000
+ha-ca-ca-oh 4 14.500 180.000 2.000
+oh-c -c3-h1 6 0.000 180.000 2.000
+h1-c3-n3-hn 6 1.800 0.000 3.000
+h1-c3-c3-hc 9 1.400 0.000 3.000
+ha-ca-ca-ha 4 14.500 180.000 2.000
+
+IMPROPER
+c3-o -c -oh 1.1 180.0 2.0
+c3-ca-ca-ca 1.1 180.0 2.0
+ca-ca-ca-ha 1.1 180.0 2.0 Using general improper torsional angle X- X-ca-ha, penalty score= 6.0)
+ca-ca-ca-oh 1.1 180.0 2.0 Using the default value
+
+NONBON
+ n3 1.8240 0.1700
+ c3 1.9080 0.1094
+ c 1.9080 0.0860
+ o 1.6612 0.2100
+ ca 1.9080 0.0860
+ oh 1.7210 0.2104
+ hn 0.6000 0.0157
+ h1 1.3870 0.0157
+ hc 1.4870 0.0157
+ ha 1.4590 0.0150
+ ho 0.0000 0.0000
+
+
+
diff --git a/modules/mol/mm/tests/TYR/out.mpdb b/modules/mol/mm/tests/TYR/out.mpdb
new file mode 100644
index 0000000000000000000000000000000000000000..ef3dca1510bc20cd082af1607caa0750f1b7e8a6
--- /dev/null
+++ b/modules/mol/mm/tests/TYR/out.mpdb
@@ -0,0 +1,48 @@
+ATOM 1 N TYR 1 1.320 0.952 1.428 -0.899800 1.55 n3
+ATOM 2 CA TYR 1 -0.018 0.429 1.734 0.132500 1.70 c3
+ATOM 3 C TYR 1 -0.103 0.094 3.201 0.637100 1.70 c
+ATOM 4 O TYR 1 0.886 -0.254 3.799 -0.548000 1.52 o
+ATOM 5 CB TYR 1 -0.274 -0.831 0.907 -0.071100 1.70 c3
+ATOM 6 CG TYR 1 -0.189 -0.496 -0.559 -0.128300 1.70 ca
+ATOM 7 CD1 TYR 1 1.022 -0.589 -1.219 -0.090000 1.70 ca
+ATOM 8 CD2 TYR 1 -1.324 -0.102 -1.244 -0.090000 1.70 ca
+ATOM 9 CE1 TYR 1 1.103 -0.282 -2.563 -0.182000 1.70 ca
+ATOM 10 CE2 TYR 1 -1.247 0.210 -2.587 -0.182000 1.70 ca
+ATOM 11 CZ TYR 1 -0.032 0.118 -3.252 0.129100 1.70 ca
+ATOM 12 OH TYR 1 0.044 0.420 -4.574 -0.494100 1.52 oh
+ATOM 13 OXT TYR 1 -1.279 0.184 3.842 -0.602100 1.52 oh
+ATOM 14 H TYR 1 1.977 0.225 1.669 0.365800 1.20 hn
+ATOM 15 H2 TYR 1 1.365 1.063 0.426 0.365800 1.20 hn
+ATOM 16 HA TYR 1 -0.767 1.183 1.489 0.097700 1.20 h1
+ATOM 17 HB2 TYR 1 0.473 -1.585 1.152 0.064700 1.20 hc
+ATOM 18 HB3 TYR 1 -1.268 -1.219 1.134 0.064700 1.20 hc
+ATOM 19 HD1 TYR 1 1.905 -0.902 -0.683 0.136000 1.20 ha
+ATOM 20 HD2 TYR 1 -2.269 -0.031 -0.727 0.136000 1.20 ha
+ATOM 21 HE1 TYR 1 2.049 -0.354 -3.078 0.145500 1.20 ha
+ATOM 22 HE2 TYR 1 -2.132 0.523 -3.121 0.145500 1.20 ha
+ATOM 23 HH TYR 1 -0.123 -0.399 -5.059 0.421000 1.20 ho
+ATOM 24 HXT TYR 1 -1.333 -0.030 4.784 0.446000 1.20 ho
+CONECT 1 2 14 15
+CONECT 2 1 3 5 16
+CONECT 3 2 4 13
+CONECT 4 3
+CONECT 5 2 6 17 18
+CONECT 6 5 7 8
+CONECT 7 6 9 19
+CONECT 8 6 10 20
+CONECT 9 7 11 21
+CONECT 10 8 11 22
+CONECT 11 9 10 12
+CONECT 12 11 23
+CONECT 13 3 24
+CONECT 14 1
+CONECT 15 1
+CONECT 16 2
+CONECT 17 5
+CONECT 18 5
+CONECT 19 7
+CONECT 20 8
+CONECT 21 9
+CONECT 22 10
+CONECT 23 12
+CONECT 24 13
diff --git a/modules/mol/mm/tests/test_antechamber.py b/modules/mol/mm/tests/test_antechamber.py
new file mode 100644
index 0000000000000000000000000000000000000000..27d0828fab8e53a7eb6b9e6d7c024c9f2099ffd8
--- /dev/null
+++ b/modules/mol/mm/tests/test_antechamber.py
@@ -0,0 +1,101 @@
+'''Unit tests for the Antechamber submodule.'''
+
+import unittest
+import os, sys
+
+from ost.mol import mm
+
+class TestAntechamber(unittest.TestCase):
+ ###########################################################################
+ # HELPERS
+ def _CompareInteractions(self, int_new, int_ref):
+ self.assertEqual(int_new.IsParametrized(), int_ref.IsParametrized())
+ if int_new.IsParametrized():
+ self.assertEqual(int_new.GetNames(), int_ref.GetNames())
+ self.assertEqual(int_new.GetTypes(), int_ref.GetTypes())
+ self.assertEqual(int_new.GetFuncType(), int_ref.GetFuncType())
+ params_new = int_new.GetParam()
+ params_ref = int_ref.GetParam()
+ self.assertEqual(len(params_new), len(params_ref))
+ for p_new, p_ref in zip(params_new, params_ref):
+ self.assertAlmostEqual(p_new, p_ref)
+
+ def _CompareInteractionLists(self, int_list_new, int_list_ref):
+ self.assertEqual(len(int_list_new), len(int_list_ref))
+ for int_new, int_ref in zip(int_list_new, int_list_ref):
+ self._CompareInteractions(int_new, int_ref)
+
+ def _CompareBuildingBlocks(self, ff_new, ff_ref, res_name):
+ # get BuildingBlocks
+ bb_new = ff_new.GetBuildingBlock(res_name)
+ bb_ref = ff_ref.GetBuildingBlock(res_name)
+ # get atoms and LJs
+ self.assertEqual(len(bb_new.GetAtoms()), len(bb_ref.GetAtoms()))
+ for aname, aname_ref in zip(bb_new.GetAtoms(), bb_ref.GetAtoms()):
+ self.assertEqual(aname, aname_ref)
+ atype = bb_new.GetType(aname)
+ self.assertEqual(atype, bb_ref.GetType(aname))
+ self.assertEqual(bb_new.GetMass(aname), bb_ref.GetMass(aname))
+ self.assertEqual(ff_new.GetMass(atype), ff_ref.GetMass(atype))
+ self.assertEqual(bb_new.GetCharge(aname), bb_ref.GetCharge(aname))
+ self._CompareInteractions(ff_new.GetLJ(atype), ff_ref.GetLJ(atype))
+ # compare all types of interactions
+ self._CompareInteractionLists(bb_new.GetBonds(), bb_ref.GetBonds())
+ self._CompareInteractionLists(bb_new.GetAngles(), bb_ref.GetAngles())
+ self._CompareInteractionLists(bb_new.GetDihedrals(),
+ bb_ref.GetDihedrals())
+ self._CompareInteractionLists(bb_new.GetImpropers(),
+ bb_ref.GetImpropers())
+ # those below are expected to be 0 btw
+ self._CompareInteractionLists(bb_new.GetCMaps(), bb_ref.GetCMaps())
+ self._CompareInteractionLists(bb_new.GetConstraints(),
+ bb_ref.GetConstraints())
+ self._CompareInteractionLists(bb_new.GetExclusions(),
+ bb_ref.GetExclusions())
+ ###########################################################################
+
+ def testParseFromFiles(self):
+ # get new ff and compare with ref.
+ ff_ref = mm.Forcefield.Load('TYR/ff_AA.dat')
+ # self check
+ self._CompareBuildingBlocks(ff_ref, ff_ref, 'TYR')
+ # check parsed files
+ ff = mm.Forcefield()
+ mm.antechamber.AddFromFiles(ff, 'TYR/frcmod', 'TYR/out.mpdb')
+ self._CompareBuildingBlocks(ff, ff_ref, 'TYR')
+
+ def testParseFromPath(self):
+ # check parsed path
+ ff = mm.Forcefield()
+ mm.antechamber.AddFromPath(ff, 'TYR')
+ ff_ref = mm.Forcefield.Load('TYR/ff_AA.dat')
+ self._CompareBuildingBlocks(ff, ff_ref, 'TYR')
+
+ def testAddIon(self):
+ # add pseudo iron and see
+ ff = mm.Forcefield()
+ mm.antechamber.AddIon(ff, 'FE2', 'FE', 55, 2, 1, 0)
+ bb = ff.GetBuildingBlock('FE2')
+ # check block
+ self.assertEqual(bb.GetAtoms(), ['FE'])
+ self.assertEqual(bb.GetTypes(), ['FE'])
+ self.assertEqual(bb.GetMasses(), [55])
+ self.assertEqual(bb.GetCharges(), [2])
+ # nothing else...
+ self.assertEqual(len(bb.GetAngles()), 0)
+ self.assertEqual(len(bb.GetCMaps()), 0)
+ self.assertEqual(len(bb.GetConstraints()), 0)
+ self.assertEqual(len(bb.GetImpropers()), 0)
+ self.assertEqual(len(bb.GetDihedrals()), 0)
+ self.assertEqual(len(bb.GetBonds()), 0)
+ self.assertEqual(len(bb.GetExclusions()), 0)
+ # check force field
+ lj = ff.GetLJ('FE')
+ self.assertTrue(lj.IsParametrized())
+ self.assertEqual(lj.GetTypes(), ['FE'])
+ self.assertEqual(lj.GetFuncType(), mm.FuncType.LJ)
+ self.assertEqual(lj.GetParam(), [1, 0])
+
+if __name__ == "__main__":
+ from ost import testutils
+ testutils.RunTests()