gromosff

pygromos/files/forcefield/_generic_force_field.py

@@ -18,9 +18,12 @@
 
 
 class _generic_force_field:
-    def __init__(self, name: str = "generic", path_to_files: List(str) = None, auto_import: bool = True):
+    def __init__(
+        self, name: str = "generic", path_to_files: List(str) = None, auto_import: bool = True, verbose: bool = False
+    ):
         self.name = name
         self.path_to_files = path_to_files
+        self.verbose = verbose
         if auto_import:
             self.auto_import_ff()
 

pygromos/files/forcefield/amber/amberff.py

@@ -0,0 +1,74 @@
+import os
+import shutil
+from typing import List
+from pygromos.files.coord.cnf import Cnf
+
+from pygromos.files.forcefield._generic_force_field import _generic_force_field
+from pygromos.files.topology.top import Top
+from pygromos.files.forcefield.amber.amber2gromos import amber2gromos
+
+
+class AmberFF(_generic_force_field):
+    def __init__(
+        self, name: str = "amber", path_to_files: List(str) = None, auto_import: bool = True, verbose: bool = False
+    ):
+        super().__init__(name, path_to_files=path_to_files, auto_import=auto_import, verbose=verbose)
+
+    def auto_import_ff(self):
+        # check path
+        if self.path is not None:
+            if isinstance(self.path, List) and len(self.path) > 0 and isinstance(self.path[0], str):
+                self.amber_basedir = self.path[0]
+        elif shutil.which("tleap") is not None:
+            has_amber = True  # ambertools in path
+            self.amber_basedir = os.path.abspath(os.path.dirname(shutil.which("tleap")) + "/../")
+        else:
+            has_amber = False
+            raise ImportError(
+                "Could not import GAFF FF as ambertools was missing! " "Please install the package for this feature!"
+            )
+
+        if self.verbose:
+            print("Found amber: " + str(has_amber))
+
+        self.amber_bindir = self.amber_basedir + "/bin"
+        self.leaprc_files = [
+            self.amber_basedir + "/dat/leap/cmd/leaprc.gaff",
+            self.amber_basedir + "/dat/leap/cmd/leaprc.water.tip3p",
+        ]
+        self.frcmod_files = [self.amber_basedir + "/dat/leap/parm/frcmod.chcl3"]
+
+        for leaprc in self.leaprc_files:
+            if not os.path.isfile(leaprc):
+                raise ImportError("could not find ff file " + leaprc)
+
+        for frcmod in self.frcmod_files:
+            if not os.path.isfile(frcmod):
+                raise ImportError("could not find ff file " + frcmod)
+
+    def create_top(self, in_top: Top = None, mol: str = None) -> Top:
+        if self.amber is None:
+            self.create_mol2()
+            self.amber = amber2gromos(
+                in_mol2_file=self.in_mol2_file,
+                mol=self.mol,
+                forcefield=self.Forcefield,
+                gromosPP=self.gromosPP,
+                work_folder=self.work_folder,
+            )
+        self.top = Top(self.amber.get_gromos_topology())
+
+    def create_cnf(self, in_top: Top = None, mol: str = None) -> Cnf:
+        if self.amber is None:
+            self.create_mol2()
+            self.amber = amber2gromos(
+                in_mol2_file=self.in_mol2_file,
+                mol=self.mol,
+                forcefield=self.Forcefield,
+                gromosPP=self.gromosPP,
+                work_folder=self.work_folder,
+            )
+        self.cnf = Cnf(self.amber.get_gromos_coordinate_file())
+
+    def create_mol2(self, mol: str = None):
+        pass

pygromos/files/forcefield/gromos/gromosff.py

@@ -0,0 +1,75 @@
+from typing import List
+from pygromos.data import topology_templates
+from pygromos.files.forcefield._generic_force_field import _generic_force_field
+from pygromos.files.topology.top import Top
+from pygromos.files.topology.ifp import Ifp
+from pygromos.files.topology.mtb import Mtb
+from pygromos.data.ff.Gromos2016H66 import ifp, mtb, mtb_orga
+from pygromos.data.ff.Gromos54A7 import ifp as ifp_54a7, mtb as mtb_54a7
+
+
+class GromosFF(_generic_force_field):
+    def __init__(
+        self, name: str = "gromos", path_to_files: List(str) = None, auto_import: bool = True, verbose: bool = False
+    ):
+        super().__init__(name, path_to_files, auto_import, verbose)
+
+    def auto_import_ff(self):
+        if self.name in ["2016H66", "gromos2016H66", "gromos"]:
+            self.ifp = Ifp(ifp)
+            self.mtb = Mtb(mtb)
+            self.mtb_orga = Mtb(mtb_orga)
+        else:
+            self.ifp = Ifp(ifp_54a7)
+            self.mtb = Mtb(mtb_54a7)
+            self.mtb_orga = None
+
+    def create_top(self, in_top: Top = None, mol: str = None) -> Top:
+        if in_top is None:
+            in_top = Top(in_value=topology_templates.topology_template_dir + "/blank_template+spc.top")
+
+        # check if mol is a residue or residue list
+        known_res_names = self.mtb.all_res_names
+        if self.mtb_orga is not None:
+            known_res_names += self.mtb_orga.all_res_names
+        if isinstance(mol, str):
+            if mol in known_res_names:
+                return self.convert_residue_to_top(in_top, mol)
+            else:
+                raise ValueError(
+                    f"{mol} is not a valid residue name"
+                )  # TODO: implement RDKit mol to residue name conversion
+        elif isinstance(mol, list):
+            if all([x in known_res_names for x in mol]):
+                for res in mol:
+                    in_top = self.convert_residue_to_top(in_top, res)
+                return in_top
+            else:
+                raise ValueError(
+                    f"{mol} is not a valid residue name"
+                )  # TODO: implement RDKit mol to residue name conversion
+        return in_top
+
+    def convert_residue_to_top(self, in_top: Top, residue: str) -> Top:
+        # find the mtb block for the resname. and store the block in data
+        data = None
+        if residue in self.mtb.mtb_solutes.keys():
+            data = self.mtb.mtb_solutes[residue]
+        elif residue in self.mtb.mtb_solvents.keys():
+            data = self.mtb.mtb_solvents[residue]
+        elif residue in self.mtb.mtb_ends.keys():
+            data = self.mtb.mtb_ends[residue]
+        elif self.mtb_orga is not None:
+            if residue in self.mtb.mtb_solutes.keys():
+                data = self.mtb.mtb_solutes[residue]
+            elif residue in self.mtb.mtb_solvents.keys():
+                data = self.mtb.mtb_solvents[residue]
+            elif residue in self.mtb.mtb_ends.keys():
+                data = self.mtb.mtb_ends[residue]
+        else:
+            raise ValueError(f"{residue} is not a valid residue name")
+
+        # create the topology by going through all fields in data
+        for atom in data.atoms:
+            in_top.add_new_atom()
+        return in_top

pygromos/files/forcefield/openff/openff.py

@@ -21,8 +21,10 @@
 
 
 class OpenFF(_generic_force_field):
-    def __init__(self, name: str = "generic", path_to_files: List(str) = None, auto_import: bool = True):
-        super().__init__(name, path_to_files, auto_import)
+    def __init__(
+        self, name: str = "openff", path_to_files: List(str) = None, auto_import: bool = True, verbose: bool = False
+    ):
+        super().__init__(name, path_to_files=path_to_files, auto_import=auto_import, verbose=verbose)
         self.atomic_number_dict = collections.defaultdict(str)
         self.gromosTop = None
 
@@ -80,11 +82,18 @@ def create_top(self, in_top: Top = None, mol: str = None) -> Top:
             self.gromosTop._orig_file_path = os.getcwd()
 
         # create molecule
+        self._init_mol_for_convert(mol=mol)
+
+        # convert molecule
+        self.convert()
+        return self.gromosTop
+
+    def _init_mol_for_convert(self, mol: str = None):
         if mol is None:
             raise ValueError("No molecule given!")
         elif isinstance(mol, Molecule):
             self.openFFmolecule = mol
-        elif isinstance(mol, Chem.Mol):
+        elif isinstance(mol, Chem.rdchem.Mol):
             self.openFFmolecule = Molecule.from_rdkit(self.mol)
         elif isinstance(mol, str):
             self.openFFmolecule = Molecule.from_smiles(mol)
@@ -102,10 +111,6 @@ def create_top(self, in_top: Top = None, mol: str = None) -> Top:
         self.exclusionList13 = dict()
         self.exclusionList14 = dict()
 
-        # convert molecule
-        self.convert()
-        return self.gromosTop
-
     def convert(self):
         # print OpenFF warning in Title
         titleString = ""

pygromos/files/forcefield/serenityff/serenityff.py

@@ -0,0 +1,171 @@
+import collections
+from typing import List
+from simtk import unit
+from rdkit import Chem
+
+from pygromos.data import topology_templates
+from pygromos.files.forcefield.openff.openff import OpenFF
+from pygromos.files.forcefield._generic_force_field import _generic_force_field
+from pygromos.files.forcefield.serenityff.serenityff_data import serenityff_C12, serenityff_C6
+from pygromos.files.topology.top import Top
+
+
+class SerenityFF(_generic_force_field):
+    def __init__(
+        self, name: str = "serenity", path_to_files: List(str) = None, auto_import: bool = True, verbose: bool = False
+    ):
+        super().__init__(name, path_to_files, auto_import, verbose)
+        self.off = OpenFF(name, path_to_files, auto_import, verbose)
+
+    def auto_import_ff(self):
+        self.top = Top(in_value=topology_templates.topology_template_dir + "/blank_template+spc.top")
+        self.develop = False
+        self.C12_input = {}
+        self.partial_charges = collections.defaultdict(float)
+        self.serenityFFelements = ["H", "C", "N", "O", "F", "S", "Br", "I"]
+        self.C6_pattern = collections.defaultdict(list)
+        self.C12_pattern = collections.defaultdict(list)
+
+    def create_top(
+        self,
+        in_top: Top = None,
+        mol: str = None,
+        C12_input={"H": 0.0, "C": 0.0},
+        partial_charges=collections.defaultdict(float),
+    ):
+        self.top = in_top
+        self.off.gromosTop = in_top
+        self.off._init_mol_for_convert(mol)
+        self.off.convertResname()
+        self.off.convertBonds()
+        self.off.convertAngles()
+        self.off.convertTosions()
+        self.off.convertImproper()
+        self.off.convert_other_stuff()
+        self.create_serenityff_nonBonded(C12_input=C12_input, partial_charges=partial_charges)
+        self.top = self.off.gromosTop
+        return self.top
+
+    def read_pattern(self, C6orC12: str = "C6"):
+        for element in self.serenityFFelements:
+            if C6orC12 == "C6":
+                folder = serenityff_C6
+            else:
+                folder = serenityff_C12
+            try:
+                infile = open(folder + element + ".dat", "r")
+                for line in infile:
+                    content = line.strip().split()
+                    if C6orC12 == "C6":
+                        self.C6_pattern[element].append(content)
+                    else:
+                        self.C12_pattern[element].append(content)
+                infile.close()
+            except IOError:
+                raise Exception("WIP")
+
+    def _pattern_matching_for_one_element(self, element: str = "H") -> dict():
+        # TODO: add C12 support
+
+        return_dict = collections.defaultdict(list)
+        for pattern in reversed(self.C6_pattern[element]):
+            # create pattern and init some variables
+            mol_pattern = Chem.MolFromSmarts(pattern[0])
+            idx = 0
+            idx_in_rdkmol = 0
+
+            # find the atom for which the pattern was made
+            for atom in mol_pattern.GetAtoms():
+                if atom.GetAtomMapNum() == 1:
+                    idx = atom.GetIdx()
+
+            # get all matches
+
+            matches = self.mol.GetSubstructMatches(mol_pattern, uniquify=False)
+            if len(matches) >= 1:
+                for match in matches:
+                    idx_in_rdkmol = match[idx]
+                    return_dict[idx_in_rdkmol] = [element + str(pattern[1]), pattern[2]]
+        return return_dict
+
+    def get_LJ_parameters(self) -> dict:
+        return_dict = collections.defaultdict(list)
+        if self.develop:
+            # get all elements in self.mol
+            contained_elements_set = set()
+            for atom in self.mol.GetAtoms():
+                element = atom.GetSymbol()
+                contained_elements_set.add(element)
+
+            # pattern matching for all elements contained in self.mol
+            for element in contained_elements_set:
+                return_dict.update(self._pattern_matching_for_one_element(element=element))
+        else:
+            raise NotImplementedError("WIP")
+        return return_dict
+
+    def create_serenityff_nonBonded(
+        self, C12_input={"H": 0.0, "C": 0.0}, partial_charges=collections.defaultdict(float)
+    ):
+        if self.develop:
+            self.read_pattern()
+        c6dict = self.get_LJ_parameters()
+        if self.develop:
+            self.off.createVdWexclusionList()
+            moleculeItr = 1
+            for molecule in self.off.molecule_force_list:
+                panm_dict = collections.defaultdict(int)
+                for key in molecule["vdW"]:
+                    force = molecule["vdW"][key]
+                    ATNM = int(key[0]) + 1
+                    MRES = moleculeItr
+                    element_symbol = self.mol.GetAtomWithIdx(int(key[0])).GetSymbol()
+                    panm_dict[element_symbol] += 1
+                    PANM = element_symbol + str(panm_dict[element_symbol])
+                    IAC = 0  # will not be used if we use automatic
+                    MASS = self.off.openFFmolecule.atoms[int(key[0])].mass.value_in_unit(unit.dalton)
+                    CG = 0
+                    if self.develop:
+                        CG = partial_charges[int(key[0])]
+                    if ATNM == self.mol.GetNumAtoms():
+                        CGC = 1
+                    else:
+                        CGC = 0
+                    if str(key[0]) in self.off.exclusionList13:
+                        openFFexList13 = list(self.off.exclusionList13[str(key[0])])
+                        INE = [int(x) + 1 for x in openFFexList13]
+                    else:
+                        INE = list()
+                    if str(key[0]) in self.off.exclusionList14:
+                        openFFexList14 = list(self.off.exclusionList14[str(key[0])])
+                        INE14 = [int(x) + 1 for x in openFFexList14]
+                    else:
+                        INE14 = list()
+                    epsilon = float(force.epsilon.value_in_unit(unit.kilojoule_per_mole))
+                    rmin = 2 * force.rmin_half.value_in_unit(unit.nanometer)
+                    C6 = float(c6dict[key[0]][1])  # ** 2
+                    CS6 = 0.5 * C6
+                    C12 = epsilon * (rmin**12)
+                    if self.develop:
+                        C12 = C12_input[(c6dict[key[0]][0])]
+                    CS12 = 0.5 * C12
+                    IACname = c6dict[key[0]][0]
+                    self.off.gromosTop.add_new_atom(
+                        ATNM=ATNM,
+                        MRES=MRES,
+                        PANM=PANM,
+                        IAC=IAC,
+                        MASS=MASS,
+                        CG=CG,
+                        CGC=CGC,
+                        INE=INE,
+                        INE14=INE14,
+                        C6=C6,
+                        C12=C12,
+                        CS6=CS6,
+                        CS12=CS12,
+                        IACname=IACname,
+                    )
+                moleculeItr += 1
+        else:
+            raise NotImplementedError("WIP")

pygromos/files/forcefield/serenityff/serenityff_data/.gitignore

@@ -0,0 +1,2 @@
+C6
+C12

pygromos/files/topology/ifp.py

@@ -9,7 +9,7 @@
 from typing import Union
 
 
-class ifp(_general_gromos_file._general_gromos_file):
+class Ifp(_general_gromos_file._general_gromos_file):
     _gromos_file_ending = "ifp"
 
     def __init__(self, in_value: Union[str, dict]):