[2] Coord gen new

polyply/src/backmap.py

@@ -0,0 +1,61 @@
+# Copyright 2020 University of Groningen
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from .processor import Processor
+from .generate_templates import find_atoms
+"""
+Processor implementing a template based back
+mapping to lower resolution coordinates for
+a meta molecule.
+"""
+
+class Backmap(Processor):
+    """
+    This processor takes a a class:`polyply.src.MetaMolecule` and
+    places coordinates form a higher resolution createing positions
+    for the lower resolution molecule associated with the MetaMolecule.
+    """
+    @staticmethod
+    def _place_init_coords(meta_molecule):
+        """
+        For each residue in a class:`polyply.src.MetaMolecule` the
+        positions of the atoms associated with that residue stored in
+        attr:`polyply.src.MetaMolecule.molecule` are created from a
+        template residue located in attr:`polyply.src.MetaMolecule.templates`.
+
+        Parameters
+        ----------
+        meta_molecule: :class:`polyply.src.MetaMolecule`
+        """
+        count = 0
+        for node in meta_molecule.nodes:
+            resname = meta_molecule.nodes[node]["resname"]
+            cg_coord = meta_molecule.nodes[node]["position"]
+            template = meta_molecule.templates[resname]
+            resid = node + 1
+            low_res_atoms = find_atoms(meta_molecule.molecule, "resid", resid)
+            for atom_low  in low_res_atoms:
+                atomname = meta_molecule.molecule.nodes[atom_low]["atomname"]
+                vector = template[atomname]
+                new_coords = cg_coord + vector
+                if meta_molecule.molecule.nodes[atom_low]["build"]:
+                    meta_molecule.molecule.nodes[atom_low]["position"] = new_coords
+
+    def run_molecule(self, meta_molecule):
+        """
+        Apply placing of coordinates to meta_molecule. For more
+        detail see `self._place_init_coords`.
+        """
+        self._place_init_coords(meta_molecule)
+        return meta_molecule

polyply/src/gen_coords.py

@@ -12,6 +12,45 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+"""
+High level API for the polyply coordinate generator
+"""
+
+import networkx as nx
+import vermouth.forcefield
+from .generate_templates import GenerateTemplates
+from .random_walk import RandomWalk
+from .backmap import Backmap
+from .topology import Topology
+
 def gen_coords(args):
-    print("template file")
-    return None
+    # Read in the topology
+    topology = Topology.from_gmx_topfile(name=args.name, path=args.toppath)
+    topology.preprocess()
+
+    # check if molecules are all connected
+    for molecule in topology.molecules:
+        if not nx.is_connected(molecule):
+            msg = ('\n Molecule {} consistes of two disconnected parts. '
+                   'Make sure all atoms/particles in a molecule are '
+                   'connected by bonds, constraints or virual-sites')
+            raise IOError(msg.format(molecule.name))
+
+    # read in coordinates if there are any
+    if args.coordpath:
+        topology.add_positions_from_file(args.coordpath)
+    else:
+        for molecule in topology.molecules:
+            for node in molecule.molecule.nodes:
+                molecule.molecule.nodes[node]["build"] = True
+
+    # Build polymer structure
+    GenerateTemplates().run_system(topology)
+    RandomWalk().run_system(topology)
+    Backmap().run_system(topology)
+    #energy_minimize().run_system(topology)
+
+    system = topology.convert_to_vermouth_system()
+    # Write output
+    vermouth.gmx.gro.write_gro(system, args.outpath, precision=7,
+                               title='polyply structure', box=(10, 10, 10))