From 4e11731d3a849ef68702c7ff156b4cb1e6683c9f Mon Sep 17 00:00:00 2001
From: Gerardo Tauriello <gerardo.tauriello@unibas.ch>
Date: Fri, 20 May 2016 14:44:22 +0200
Subject: [PATCH] SCHWED-924: incorporated Antechamber based ff parameter
generation
Mostly based on scripts written by Niklaus.
The functionality is in functions in a new submodule ost.mol.mm.antechamber.
It allows for the execution of Antechamber and the parsing of files generated
by it.
---
modules/mol/mm/doc/buildingblock.rst | 2 +-
modules/mol/mm/doc/forcefield.rst | 65 ++-
modules/mol/mm/doc/integrators.rst | 2 +-
modules/mol/mm/doc/interaction.rst | 2 +-
modules/mol/mm/doc/molmm.rst | 2 +-
modules/mol/mm/doc/observers.rst | 2 +-
modules/mol/mm/doc/settings.rst | 2 +-
modules/mol/mm/doc/simulation.rst | 2 +-
modules/mol/mm/doc/topology.rst | 2 +-
modules/mol/mm/pymod/CMakeLists.txt | 3 +-
modules/mol/mm/pymod/__init__.py | 22 +-
modules/mol/mm/pymod/antechamber.py | 511 +++++++++++++++++++++++
modules/mol/mm/tests/CMakeLists.txt | 1 +
modules/mol/mm/tests/TYR/ff_AA.dat | Bin 0 -> 52853 bytes
modules/mol/mm/tests/TYR/frcmod | 104 +++++
modules/mol/mm/tests/TYR/out.mpdb | 48 +++
modules/mol/mm/tests/test_antechamber.py | 101 +++++
17 files changed, 843 insertions(+), 28 deletions(-)
create mode 100644 modules/mol/mm/pymod/antechamber.py
create mode 100644 modules/mol/mm/tests/TYR/ff_AA.dat
create mode 100644 modules/mol/mm/tests/TYR/frcmod
create mode 100644 modules/mol/mm/tests/TYR/out.mpdb
create mode 100644 modules/mol/mm/tests/test_antechamber.py
diff --git a/modules/mol/mm/doc/buildingblock.rst b/modules/mol/mm/doc/buildingblock.rst
index c9dbcfd7d..1e9f7bf02 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 9c5754152..9ddbe77ff 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 ed94d8695..5c4d35241 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 95d193e1a..1b0a705d9 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 f243b7dfb..2de5d956a 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 1aecc69fe..265236442 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 648be8bcc..d8b3a13b4 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 217c1cc7d..9d20aac61 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 421be2fe9..b61647996 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 3156429e1..1f47d275b 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 5c1138cd0..a4a7c650e 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 000000000..4776b2a96
--- /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 99f2c4575..b3eb4d801 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
GIT binary patch
literal 52853
zcmZQ%U}&%hVF3mP1||jw$VgOWYqW>+OptlS#f7#EPB6LT<Zl`djxb(EKC;S;WMp1)
z3W`cYWVwtyWL|y-149Gc%;ZFnK1PVc81g|<Fb;~{$;KdcFqh;de$z;Rn~`S>;;}L?
zFfc>S35g7nU|?W?ndAr8XY3uBwjZP#CgKb?CD~YH9lIUK7MM0?7$f<Qgt9%dj6X~!
zf0rlEevqT#<{`q)_`#fsdqKM3a_%rza-wh|^FELq%*D<wh6u?wfjht&5t2rDBwdko
zn#Nh78taOr6DA2W%Q*_hN=`gw!~*s^OvK*<#>&rF`|Hm>kRM^r@s9{WNLndNfdc>*
zB5-{fdCFe9?T};M1MU_qa`1S_F#I5P!WN_tW}k->LSu4%V+)$3G0ZU;n39N4%}AW)
z%nf!9+*lWcv6zyGP|ZmE@>c{^Cn8i4k}y|zAVMHNp!YheMnuSE<jY1q0ZW2Xsu}|W
zgCzq413NSU!zfUwGD31L$b1+ENhibood;(*JAyQTWc{?;#q19MYw`dI;4mR04`d!l
zb)GRu8Hfc|4`TSW)?Bc?!?qACK!gD>qdbgYM&UC7rY0GtEg9|tXPA()V_>$Fot{b?
zJ|`eT0l`A#X=f*CeObHK4duickpa(12qQck*WI{pt76}Y%?P+RV2pei173FcPgwB7
zHZ6(86C{HcjvygWSi?ghKLaM@3`?b;%(jQsu5jW$Yz}~h6~gn$i7+u|cv0l+p8QVA
zPD-zZI3qIP?s0a4Im5%r$YIHS+bI{Zh6Bt>gfp<4VN8q}a1S7P#U17jXBWeGc0Ri$
z6^HTqg@74GB$|O_JDllCszcy9U5wyNSF&A#D4(38$n_1P7$eb)j6}E(Tnu4M4;Mq}
zdG~BjtnDE#YLE?qhlq!Z5xIuAB3$B1wo6<QE^#HzCHWb!T!5%o{XLi~xb654&BYfY
z8Tl}c2qO@ch`-13%Hy`0EYpcK0v;U^Auy(YL<r0fkX4M3h7w2?$^{id4Dei!$i;|y
z0jX4m6;XbyUO((J{J349hJefiVGtLDL53oXLa>mkD?|x}sJP;%f3v?@#D&jPL^+Mk
zR1bJl0;xLmp77QF%F=USZ=m@f)2#>>BC1FvvtR+>0Wws{@F-Obg}K7Rh$vHGZb*hP
z;E4#~P(&T{xS7*oRv;g~Fhx|r@N|Gkj!0&LVi**&{t+SeD}LF3krKjhD$HMqm_jPn
z5ix~G3GQK;zwB29vXE*lQjLRTEX=tcPDV>^eYKZhK15t>A)HIdTw^@u!XgJAqzGdX
z;SF{WC`b{F4reEgx8LmV7O`Rr1U!aj7{Vm680z78y!nIudXM8^IkZdw(+77cBJwj3
zjT%^J`XkB~kBE>rci-Dz(7lGwShx{L=EJOncczfrJaSH)4lEB9vH2fMih?(Q5hVjs
z>j$13JRC)yeYH>eaS)qRu^5Z36@ipcVWHvdVz_f3o5PY`tkFe@xkkj8i*O|(7~sL|
z4DWNIwB@b}f3xSVK7uzK$S~KKG;<L)Ah#!A0p^U<=rj7q$mWpL#e+9JAQ=ZU1F5V*
zG7a8PaCSAMs=@HagR?97E=EdLa3*qVk8BsiO+#sb!h*^<igJUIoe8&%?4ZbiM<{Yq
zfJdH(i(wuwo5SXVr-)3(aO03#o`^=HGfGRAY?I-^iIm>p$-=|Mh-`zA!V}I!L?l&A
zMudPfB3wLN$q#0DU_`;4OdW&a#vww$8Kvn?b}(aeG9vz5$&J7K47ih#GYQN#e~-+c
ztPcJP_#0F31dJ%wkj;e;XCZe$NHQ0>(uX+>*)W)xhsXQAFZOxKXYl1bMB5p;Jq@e%
z!QB~9s|P+f>tVDi<h}ijGi&e~3a{A^O{t7zSlA+s6e3r5#;JUt>=*uA;>-dW6b4Zs
z4AKU34`vGmVF98WglMb4f&$Sr_G4gRaJbFy=>YN-$aqEu1`q|pAe&%15Q98e{E}x3
zQVKG`KOzL?HgJG|#7MUckx~%eL9h@JgWOJlITXb<2e7yH%~i7p`3HyB5N&0I9S9b}
zZz%CZ59{D9BjjcwEO`AG7>M#6qR9@AIII~4-j_vg(88kCBO(NrP8R-<bi@;wh=4<4
zVJQNTyAANLM5=uZLH2=C4+zu2J|hD5!GawghzRc?(ql#*Oc}Ps50mqN#a!yD1=!Lh
zOdSPQ!9vu-h<v*c9!ySzISI*IFe{MKraR0s92pK~3KowcEX#lo4q&wmo+6##UZaa`
z#zguJUx|#!ir7kIM1X>lBAy}{7L53;LiDNeS>@pfnq>v$ur`Y@2Q0aax>h-WE1K+t
z)pjIVg<M?09E9){EJ2|5?qNZQqqKq<gBgV6Sq7S}W?)cA{)a6PVG)Z+;|LaV*#T3B
zl%wE<vLDRFgbELsAqXoFETS!g8A7o5fGNb=I!GqaI)EAC>;#K7P>RHoL|_ULK0<0J
zV5?)W4|~8;jI%4ak^$v3jx|L@R%%$RK{x^_=8&o~cthOT6<o_uY#~w?6Tc5VT&Uwi
zLKeb(=;1=N4-o+gZ{{H46=_5aX($L5|A@8&r~o6b6@;)3FAHm*-5I&_152ZXk_=vJ
zkdh3JgaY?#6qPK5r+UQd0%Xs^qmYC|O2CUqHo}?6A&C$qBhg^B8*U<Ekp!|I;e#+F
z<b165z#=CZX#@~C0AcdZ$Rh|;iASWIh?H=UJq8N{58_KVyne)PBeHK1Hc}@Zkp{X5
zq*bzt8!Wy>wh=KNfh~YZC{eN4gUI$sHsZ(<Gz?Yb8FZL6sAFf;POHfI0B$2<;EsBs
ziY@0M(kf!yjrf#<b#wqJZeiZRK01J$hCqP|YT1A=p4JW3EClC9kYylD6${}tHTH%M
zxUm7U420=nA+bJ0<U4FEq|r!3l0>QVVX4aB0~CkEwBL{t2uw4=YGex$i4!?j!t_vN
zA#zs+mP&DW3ud1O+=J8q7vSy7VDTV!`w-<ReBmZynI_TpA!;E+Fd|rp0QEqO)Z_3f
zkugZTg;a_`yj`dacnyzzfE6Bpjxg^L%7XAkgNQ~qw(cR;WhtcgMG>ueGAx684L;t3
zW7rI)5SwiXdl988ve#gGk;mWR_Td;`g6W}_eS}7s;Pzq3dWfumXjph46=q0GL<S+g
zg@`{hAX|yZ42ZfGb&d;`KtcT{(CT*_gYt;6GJIn+h{1lW141yn?zl~Mz-t#QY!Tsz
zjfF@ZNDJSQI~d^l2vj`~vXDR&!kmL-HX_<kr+kU^A{jPfi%X<-6NP?+*GNchd~D4`
zq*O{Z8<FZHxHC{DVyP0Uh*}A$aqW!1ng@9i<|JR&5E01KAE+>AfKNH*85hK-ld&J9
zEeh9u5cyVa@OB7zx`PiZB^$3-eFEOr0hfah$Rr!@-qE`Uv>yZ}7YttM0&;C}anaBI
zeaO3VoLyn6l8xU4T7b=j`4qHV9HfhZfgwLbeSsv{G4NUjrY|FJ76<OV85nZ#1ymV^
zvHuREZNGr0fQ)4J57W>j;muqulI{qdhHg{JP<6T^bYe;(`gs|K82dpGsRElMlCc;%
zkvC`P--T?>5N2RtP+(wS&_GH)(C`IsFG5?<4+~&KEs1CfAWGQL%^9E|Byn>FB1n<D
zOtjtKff=2MCJLfP@o+c7vqywjGmObK!xfhIJzNc`U<TQKL2?H^zrYG3_~azKRf1Tr
zOy@l_h?Xp(T|&nVFYq#{xEPjL5X&)wU2-a{>|$Cizzb#w7B@&nz-tMjj6mLn08dbe
zWxk_Z2>d|n@<+E2fKn*h+HsH&s6a<_G~mk;5Ca*A+7EeE47pnf;8_MSfrVTvM5=wa
zx3d`*TL=(=i^)PNx8O}lloAZj)&xw`5EU#Tb4Ry44BIUa@Qj7j(1CYhP?|iXZKWXE
zTw^lLMXLCa+zU(P!7kMYe%qHU`43(Rg|qyEDMuKF$P~y`GEAMbtC8h>d(O3g@fnK<
zNrYj=#qc^7(G5bjVXDg)`}1iVDKQrw<kU46(fLC3LxWu?cQ2yfMrCsm;ega4BF4Q)
z+v8yoj#!omZjOO!J@{Y(9o8>{oPi_pVoPO+v2G9KRE9`;#_RWdvA_LxF{#N6$zbH<
z05{mxh*kzu!^QA%2WR9i3(Un>X4_z07I;qu-o8QP=%Ktu9G*iE-a)XC1~8Bs8K9I1
zO8tmt9kHvZ5seFk9SGLw8gT{&1`5}R!=n>XjbU#l!y+-*g}{<rm_kH6VPheRAnew`
ziwAdNCxo!sf)FF#K4K??5N2azA=!sDC%{69;8HqRd4jbVftiA+YQUv4T4sY8L1+;x
zOdSPQ!5l<z2`e7Ekb@JZ49QzCE0Fte1d=4o6f7P?ST?!>p3p*ISTN$N+D2EvL-uuz
zu7GDiWGzIq0I|;xDf=Rp+r#oU!Qum^5byF0TCRY{Qg0&si{v+W;g6W4b9RLlR|PGb
zu<b;Gc@9yBAdE&-(+J<8G}vL9X=)$b_sHEvn9bDpA9B9|W($)4;2WT*6`t6;moS^1
zsc#=5JdwJcq}zwLenT2YK%9|>J*C3ZC!r-cMB9h(AJ*{}BJIPQ*%95$(Y0C(BYmwF
zmSuK`A{1#0IPypbJi8&EV@G7l#8y(1Wg)Sh58^FEWIcFQgfi?;bsrM25XFahSBK$l
zZep7$BDfL_)}st|4dO#upaCC}%QF_akV3}#9^B`>t~iE!eb3Q<n{1KS@c6@HCqHBH
zB*-zTFc+h&)0r2CbDd6Z*k14<p|F&Ue1<B99MU=+yUH7=3u+McLxy4KK@BuXcmyLP
z;Wom3k&(aZ<0Y_bV7~A`UNQqa&y<&efdRDIN)9PTK)qKG2E{+YIc!AJ5@~d4^gL4r
z1_lzB%pgh>#DEOm^GsoJj#SRUi(QlwhTL@~NagSlTO)#0KEc~N{_cBNh3xKxHGn4r
zaHMPGg%z;y^?>Ioe|Ol)jCisJvJpsG55^o_n=!gJWAGgBh?G&_^%0^=;_uE^@x^}X
zYF_aARvdF#h~fvqLMqGP`472NBXimq$yg-U!g7{}lM#)~#W&3jDhWW9CDyJOOcLwp
z>FBAEB%B(FNT`VN6*-|I3hUhqKiYe-;XUUSk&X~a7O`p_>yAc$_pHjV_Bu{i!BL4L
zD-bdl-f%;9IWgwaWO5tcSU>~<q7RP}2#6d8pV1w86R(JdE+V`TEX2?Vyg!2|@Q|aM
zD&0NmPrPDD3<%F5TZbhLBQ^Yq^c~*Che#P5J}`~E3=y6hDLR}BYk`7PtH25rM8D7j
zKJ1vCpp9?X7!hy?7J1e|&esWf+=|UQSmlDH=0T*@VLCC1$Zp{B7*qk1mEG`!AtLS&
zQ)*Z(gQrI-9nFRBSZL&a9m+rsndLDeBO-=xU@0FtKw%rhKmkprRftf;XO+k3(OjS~
zW?(>^@Q++{z|t!s2w|xY{Twhz(J^{77g9hXax`M<4*y&Q(MNM(YxyCpL&PBV&DXFT
z=MS%8KxqI=K7%Qwo^_C9f}@6kS%)ZXvD7t)HY5J91m$eV)I1rUL$1(Zu0jq=SiT_n
zyf|zjNwj5fm%x`uQ8w=bZ-gOx4Q2^)(GRze&;%jeL_}glun>MHWFMib2IfYDi3nQ|
zi9aJBR#>6*E?}CmwX_OajIcF*U|NW`4@<KZmwf{<ClV4AAdNHv5AQ?d8CzaUC1Xy6
ze-5ra&tZJL;BCgRAVl;{l8yPfO~Lz<;d1bi&*bDFbp_PR>Ja<alT$u^IR(}Ti(hBN
zL1f8^$L?%Jk@Q2XCdf;)_~En<bj2TBCt~V9Iq_78BZ^LES9mlgCmM3BL6wYxFBMKs
zoRKGl>e?tIog8aqP$d!bG|7qU+#jMCi!%3fI0Z5>1S@VlU?G)}CwP1pc(w&5f-n<{
z9K4@F;QqYaqjfmvev&aH5iycs!o#Eub`C6DJrIGOkr(xKCYmHXbTSfmdGer2BEmHz
zQQeROO%h2blQui5BvKe8-c?~kl|-J6`VssXEGfglz@W##z~Iinz`zBnCXg^F%yEp`
zBMKWt(T&t89i5G$!)z2NdO)!?h;|Am!&@teLeAMexn9za!kh8vcDo+j1Ia0{as{!i
z%h|>7V=|9j+Pky(X0(tFbb+6|k5oK6yBA*<vb$uy0<Rfx+w+i4euSUR<mU=@NVV5y
z;ta_!ft5<AXS$@Exoca-(~s8>r1M-5h9FLCA^(a#<a2wB$-koy?x-m8FAId98|)lK
zwo9TAVd6oyOQH}iA<+<|bb>fw#?RHT(3Z<CQS$)4C`PPMMHqr~yqTY?VGHkh+Zd%~
z_zXd$dc-MB&aQCpc(@vqdvPHmU_GKpyK@l9C5V9ah$7o19`Gpk@F2z|h;db<xsK68
z1CSaDqlX3zo_h!3RUh&&B%%`h^6jU6Dig7HxFH&;SdE1h+N1XlGK}6kNb0?V@DzpA
zMnf!mKt2f9JD1f#VJbV3izyK1!rMH^%~)8Pb#^hjF^|XL>^w%i<|0hPmKWgte3U_^
z$1A_v%hceVNk=rX5q&Vko>9c+1Y~oKNi!Fz?TN_#2;U>OLc@x_*h{uc5E&3~gELI<
z8tma>IL}<bA*~R9NWfzW(e_A2w73ziAe6?pmO7h*2ru6IbCHZgSb=0Pyir4e!HCq6
zf(Tjoa0*J7qiE`1`|W%1-%bm6XC9oDj4&7x6n-eDJN6&{WB<gQA8+iz;|yVPMjpb+
z$cX_a>fve_q|f8vznj<#G~sSWv_GAZdk64lm!GTQBYhSJ>4zj=1cziYqCJUpUMTq{
zBOT3;2xk=M8&iIBB-YSyj-tjbn()9yYOOm*QNhhfp@9_S9#n8MqFdpNXxBSa!(<eP
z!qPMOp@Gx^L5@Lq0pj6mO!iHTNUno35dlLLld;C2M-;i2JR*e#q9^5y2y&{Jj0ggx
zwBZp&g>Xi8Ga?45VKO2JoDr6Icu*mnQL+atnGk0(qK-y1vk;{y_A{UTTn!^uvO3(5
zy98dhfn(7gqJG6%ryx=tQf=koYRKow<ACFOZdhvqX(R*D=z?FB2OhovjhMnGd+2Z`
z4VF1&*uX8mP31%Ryc9&+1|A8Bwkpzm0a6<n7J*1Z14J%tL$rMnb|6^D>*x@zD@11j
zzTC}^>ej*ij?{-B**e5(A$Xu7&+cLd8}*m%A$2(5Z9(iK1hDeO8FU#K<gzSKZo`rT
zVD3gl6cP*WQ$#sI$UgY&JP!L{ISyeoJa!Q2ahNXLONP%6fy-o22w@LS0yzLr7$PEX
z^hA--6GcENjMBAqkY!7wCyEU36GafkF}!9$v<Q)^P2?dxSS90!xOFb%aVwE`YGSbv
z(aa@aA=1@hge-)QXCP`Yq*@f-cR(6`Kp9Db)sBd@kY$C|q~4K)-9m)ji1k9qV}`I*
z$Am0I1SGs7MXChgYv7QLhWiS>6bH0E5z<T`u&55nFIW-^B9Kvr&5#lbjEQvPB>$Yp
zctaH7K`b^Rx*3FQM7mRicpI@KAEa!7lo^nQZs1I$d;z-n24CX~i#>?2L>fs&>YgI{
zW+*8ZmTahIBO;w5_2-e3EU`8s;t`8)5w{H?$0Ix_kp{#(T#bpe5o>=0;R^VmEY<~?
z$j&f?wZ}m9DV}H~-a>d*#Zl4`X(3V?L4+w5AL1(B;JFdhO@LJNgo+!4JqQ+78xiRS
zIk8a3MnoP&8d^sl)`KS;q;U&yUkKt!0xOlV1uIr7v1UiItR#?h;46;^rd2}yFId@6
zyoChPDzO$K1u9Z{fjI{;oP|=FAte_W6RA`szVgRnBl1El#CQZ!kq)0vp{%)qNI1xo
zHSh$Alr@nO4pL=BQCdZ^l2BTu%u1x#g{M`NnHPkIksJa#X%CVkDXXrqTZt%vP{v+~
ztUZu@1WTpJVzBgx6oDunMz)=XR-%L^+)Yu$){)3Q!tYCjf3Oy`&Qb7-@Te7+NbL=j
z2t@djI#wc0hT!rV%t{X``4Wefh%iQY*n>)OiDV^8(-G!1vYJ@<lMi+)v8GrLDisCD
zR$>cH)L9*><Vz$gv3d=4@{8(LB7BLMutJ^5qqddE_8@FW<W}Na==hTlQkwu_I7%4{
zuT+p41r(J&NLC{I5^fJtVNbl3h`tRr|01mPKuSJHOhj%a!Akh59^@HAL;(-W^f>1Z
zVG4-vr;qNFWI*0Y2pgkEn&5<W>Ik&xu<h$Z8qt9_Wsz-z#VqntGkC+Cz>p5L>?3-g
zB;K@%Xt5GZn@BBI0wX&3b0DGs_dqIVkeK9D#zb3*l+%%zh|-7nnvr-b5v4Z5Yq)Dh
zdhQ9uYZok)BNe0Y+<++O5P1MsR}xksf-c_0XCd5h#K<(l8brSiF;s*z$r5WLVoU`a
z3u)maQjLU&Oyq$U(rv^#x{6c}z?q1P#j*ET2?Z-QuffF-)gf|CkCZd8R)-W>iOByb
zRwBX`Q4}CGZK!1>!aqnk5PJiTYE~k=hA3ySH&dx)B~k<uDsn(p!qS>!aDV{zt>J3H
zxVMJ4z816t9d8S(gW=b;CQoS>wMUkLA0?81_6!fst>Jev7{Gg`VdamrJIv(dFF9&p
z17QyKcY&W^lRwQFasn?x4nEG5pK<j$WH&Wj60s#XKLgKJ>q#;Dz`KQEZb8~+oQEL?
zpWn$al>hPuZ8tUCF&W8?EJx8K;e&oyBoW0-Mn3HJX_z-q&e;ul{1og~eg*~xaRvqk
z6$S<dn7ctk=^zXW$<dw7SV93L1PUN{C4^`LBG2UoW=q-WskDI?)Da9vBqQKSgVrbY
z!u<hbz}MQrv!y>Vw<yE4!%`5U+VyuaB*p}oMnq2qd5IR&9#*@;iTkhxB+P1fR3H|9
zAv?kdj}b7t5ix)mfJZb~=y#Pf*8Px((Jy}&qtVlTL1h{!-ar_X6R~xa5Np<uyJN5f
zGkV%D38($S0xcgF-q_|w{ap+nH*=C;D$FpXEDK-yiJbk3y5AY*L~I2CVkIkl4$0re
z2#>KCh9SxZ#9C9V3kAtLJQra+a{U1dA7@8Arot2<j6krE$}yxm1(bI{d1lb<Dn*>c
zhb3sq-v9~=A(Z+69tbW*L|ns+B~fJ{>Pnbds7KuqVJzGTB=ccO58nJm&Ixi(IFBZV
z>BD9w+<}Obcu)cZkFl`8LnJB$3sLAI6$kKyhIhXz%(ZZjBbf>pL(U5@FL*c^aj<{I
zwf&O<V~z0`3kw9qiY7$ef>;j?Zk~WL7Sg!ls*v~gGtR65r)M1dJP?@_ISi2e3$q=$
zEn=L?_sM?Y&n3=Smb}0M0W&2aEEs|tCt%)!Cl7>g5G+KzASWJp2*Rru>ThI#ciE8h
z0?f_Sw+`M3q{KQzT?P+SL|u$%h9b6`AeTk(6iUKz$FT6gZyomP8g8Epd}^8GmCuNT
zg77+`phYATIExzV)R7MLL@s|3cEDMVu#rJfDS{<u!xSQVN$`L~_8Bas{1Kx8*mDDw
zt&0eO#V@$L19<{@4+R0+V0y8536ZN2r7e;bFr7pnvI)~bfmN`m0L_g<R_EfFI>KWY
zA`OmiK5-&6MFk5-e6<>)0f{eDgBDanj+bq-2*b9#nYvavINav<bjVIvZHFaC<F*R9
zAcm<z_zIrI5$7ZkzV91m3}z4_TLw28xv)YA!e=);C@6=K?Sr|Bu;(C=hU4gHm@!oJ
z9A5ii1psov3ll>wuV7+W>uH#y30=SqGlU|`V1`guWWh{9YW!o}p$bz6%59JW631LO
zEDSIs4&htG&5ZuYsS;)fa?99=fdLY4lw@xBY#u@DjOl9~c;Nt9<q27V2`-fJSO<#$
zybTiM$V0S}krEB8aPa_LQ~+6X3TiuGDK4pPAABh$MfSm0Baa?5J|Zrd#kVfO7<X}o
zsPqsWDy)?rtf35=ZKd2U_-%<OWf?5b`6Fs^kaeJZfDMCYdO0Aoy)X}g@*)U>lxHMz
zO10Vrb~M=HQ)Cj@)j&v5aq)qYZFZn@|M8fVoUA6W&hD1g4m&)GGV%|UnAw3`hfh)R
z*X8YYH{09rIXpRq80!qV1ZxSIl#$o)>6R@4>+&-mGR(EJj!Cq~<M8A}JZ^$TQNCBt
wbi3Yn3D^vRsUapPlZ`=^bw<C!7v6b^iD}F1ip0#Y*#=7}dB!kBpoqbT0Zx6WMgRZ+
literal 0
HcmV?d00001
diff --git a/modules/mol/mm/tests/TYR/frcmod b/modules/mol/mm/tests/TYR/frcmod
new file mode 100644
index 000000000..956433fb0
--- /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 000000000..ef3dca151
--- /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 000000000..27d0828fa
--- /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()
--
GitLab