From bb94017cdc64cbccfe5e423d5de337def671f298 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Fri, 29 Jul 2022 12:05:15 +0200 Subject: [PATCH 1/9] fix multiblock input --- polyply/src/map_to_molecule.py | 39 ++++++++++++++++++--------- polyply/tests/test_map_to_molecule.py | 30 +++++++++++++++++++-- 2 files changed, 54 insertions(+), 15 deletions(-) diff --git a/polyply/src/map_to_molecule.py b/polyply/src/map_to_molecule.py index eee77466..608202f1 100644 --- a/polyply/src/map_to_molecule.py +++ b/polyply/src/map_to_molecule.py @@ -13,6 +13,7 @@ # limitations under the License. import networkx as nx +from vermouth.graph_utils import make_residue_graph from polyply.src.processor import Processor def tag_exclusions(node_to_block, force_field): @@ -135,14 +136,13 @@ def match_nodes_to_blocks(self, meta_molecule): if len(meta_molecule.nodes) == 1: regular_graph.add_nodes_from(meta_molecule.nodes) - # this breaks down when to proteins are directly linked - # because they would appear as one connected component - # and not two seperate components referring to two molecules - # but that is an edge-case we can worry about later + # this will falesly also connect two molecules with the + # same molecule name, if they are from_itp and consecutively + # we deal with that below for idx, jdx in nx.dfs_edges(meta_molecule): - # the two nodes are restart nodes if idx in restart_attr and jdx in restart_attr: - restart_graph.add_edge(idx, jdx) + if restart_attr[idx] == restart_attr[jdx]: + restart_graph.add_edge(idx, jdx) else: regular_graph.add_edge(idx, jdx) @@ -151,17 +151,30 @@ def match_nodes_to_blocks(self, meta_molecule): self.node_to_block[node] = meta_molecule.nodes[node]["resname"] # fragment nodes match parts of blocks, which describe molecules - # with more than one residue + # with more than one residue. Sometimes multiple molecules come + # after each other in that case the connected component needs to + # be an integer multiple of the block for fragment in nx.connected_components(restart_graph): - block_name = restart_attr[list(fragment)[0]] - if all([restart_attr[node] == block_name for node in fragment]): - self.fragments.append(fragment) - for node in fragment: - self.node_to_block[node] = block_name - self.node_to_fragment[node] = len(self.fragments) - 1 + frag_nodes = list(fragment) + block = self.force_field.blocks[restart_attr[frag_nodes[0]]] + block_res = make_residue_graph(block, attrs=('resid', 'resname')) + len_block = len(block_res) + len_frag = len(fragment) + # here we check if the fragment is an integer multiple + # of the multiresidue block + if len_frag%len_block == 0: + n_blocks = len_frag//len_block else: + # if it is not raise an error raise IOError + for fdx in range(0, n_blocks): + current_frag_nodes = frag_nodes[fdx*len_block: (fdx+1)*len_block] + self.fragments.append(current_frag_nodes) + for node in current_frag_nodes: + self.node_to_block[node] = restart_attr[frag_nodes[0]] + self.node_to_fragment[node] = len(self.fragments) - 1 + def add_blocks(self, meta_molecule): """ Add disconnected blocks to :class:`vermouth.molecule.Moleclue` diff --git a/polyply/tests/test_map_to_molecule.py b/polyply/tests/test_map_to_molecule.py index a827f264..a9c52d85 100644 --- a/polyply/tests/test_map_to_molecule.py +++ b/polyply/tests/test_map_to_molecule.py @@ -147,7 +147,34 @@ Interaction(atoms=(7, 8), parameters=['1', '0.37', '8000'], meta={})], [(0, 1), (1, 2), (2, 3), (5, 6), (6, 7), (7, 8)], 9, - {0: "MIX", 1: "MIX", 3: "MIX", 4:"MIX"}) + {0: "MIX", 1: "MIX", 3: "MIX", 4:"MIX"}), + # test two consecutive multiblock fragments + (""" + [ moleculetype ] + ; name nexcl. + MIX 1 + ; + [ atoms ] + 1 SN1a 1 R1 C1 1 0.000 45 + 2 SN1a 1 R1 C2 1 0.000 45 + 3 SC1 2 R2 C1 2 0.000 45 + 4 SC1 2 R2 C2 2 0.000 45 + [ bonds ] + ; back bone bonds + 1 2 1 0.37 7000 + 2 3 1 0.37 7000 + 3 4 1 0.37 8000 + """, + ["R1", "R2", "R1", "R2"], + [Interaction(atoms=(0, 1), parameters=['1', '0.37', '7000'], meta={}), + Interaction(atoms=(1, 2), parameters=['1', '0.37', '7000'], meta={}), + Interaction(atoms=(2, 3), parameters=['1', '0.37', '8000'], meta={}), + Interaction(atoms=(4, 5), parameters=['1', '0.37', '7000'], meta={}), + Interaction(atoms=(5, 6), parameters=['1', '0.37', '7000'], meta={}), + Interaction(atoms=(6, 7), parameters=['1', '0.37', '8000'], meta={})], + [(0, 1), (1, 2), (2, 3), (4, 5), (5, 6), (6, 7)], + 8, + {0: "MIX", 1: "MIX", 2: "MIX", 3:"MIX"}) )) def test_multiresidue_block(lines, monomers, bonds, edges, nnodes, from_itp): """ @@ -166,7 +193,6 @@ def test_multiresidue_block(lines, monomers, bonds, edges, nnodes, from_itp): # map to molecule new_meta_mol = polyply.src.map_to_molecule.MapToMolecule(ff).run_molecule(meta_mol) # check that the disconnected molecule is properly done - #print(new_meta_mol.nodes) for node in new_meta_mol.nodes: assert len(new_meta_mol.nodes[node]['graph'].nodes) != 0 From 67e967ceea808399082f8ccd0636f4bb1267fe61 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Fri, 29 Jul 2022 13:13:56 +0200 Subject: [PATCH 2/9] add error message details and test for when multiblock doesn't match resgraph sequence --- polyply/src/map_to_molecule.py | 7 +++- polyply/tests/test_map_to_molecule.py | 59 +++++++++++++++++++++++++++ 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/polyply/src/map_to_molecule.py b/polyply/src/map_to_molecule.py index 608202f1..74bd467c 100644 --- a/polyply/src/map_to_molecule.py +++ b/polyply/src/map_to_molecule.py @@ -166,7 +166,12 @@ def match_nodes_to_blocks(self, meta_molecule): n_blocks = len_frag//len_block else: # if it is not raise an error - raise IOError + molname = restart_attr[frag_nodes[0]] + msg = (f"When mapping the molecule {molname} onto the residue graph " + "nodes labeled with from_itp, a mismatch in the length between " + "the provided molecule and the residue graph is found. Make " + "sure that all residues are in the residue graph and input itp-file.") + raise IOError(msg) for fdx in range(0, n_blocks): current_frag_nodes = frag_nodes[fdx*len_block: (fdx+1)*len_block] diff --git a/polyply/tests/test_map_to_molecule.py b/polyply/tests/test_map_to_molecule.py index a9c52d85..0a3967d7 100644 --- a/polyply/tests/test_map_to_molecule.py +++ b/polyply/tests/test_map_to_molecule.py @@ -353,3 +353,62 @@ def test_riase_multiresidue_error(lines, monomers, bonds, edges, nnodes, from_it nx.set_node_attributes(meta_mol, from_itp, "from_itp") with pytest.raises(IOError): new_meta_mol = polyply.src.map_to_molecule.MapToMolecule(ff).run_molecule(meta_mol) + +@pytest.mark.parametrize('lines, monomers, from_itp', ( + # we are missing one residue in the multiblock graph + (""" + [ moleculetype ] + ; name nexcl. + MIX 1 + ; + [ atoms ] + 1 SN1a 1 R1 C1 1 0.000 45 + 2 SN1a 1 R1 C2 1 0.000 45 + 3 SC1 2 R2 C1 2 0.000 45 + 4 SC1 2 R2 C2 2 0.000 45 + [ bonds ] + ; back bone bonds + 1 2 1 0.37 7000 + 2 3 1 0.37 7000 + 3 4 1 0.37 8000 + """, + ["R1", "R2", "R1"], + {0: "MIX", 1: "MIX", 2: "MIX"}), + # we have a residue too many in the residue graph provided + (""" + [ moleculetype ] + ; name nexcl. + MIX 1 + ; + [ atoms ] + 1 SN1a 1 R1 C1 1 0.000 45 + 2 SN1a 1 R1 C2 1 0.000 45 + 3 SC1 2 R2 C1 2 0.000 45 + 4 SC1 2 R2 C2 2 0.000 45 + [ bonds ] + ; back bone bonds + 1 2 1 0.37 7000 + 2 3 1 0.37 7000 + 3 4 1 0.37 8000 + """, + ["R1", "R2", "R3", "R1", "R2"], + {0: "MIX", 1: "MIX", 2: "MIX", 3:"MIX"}) + )) +def test_error_missing_residues_multi(lines, monomers, from_itp): + """ + It can happen that there is a length mismatch between a multiblock as identified + from the residue graph sequence and provided in the blocks file. We check that + an error is raised. + """ + lines = textwrap.dedent(lines).splitlines() + ff = vermouth.forcefield.ForceField(name='test_ff') + polyply.src.polyply_parser.read_polyply(lines, ff) + # build the meta-molecule + meta_mol = MetaMolecule(name="test", force_field=ff) + meta_mol.add_monomer(0, monomers[0], []) + for node, monomer in enumerate(monomers[1:]): + meta_mol.add_monomer(node+1, monomer, [(node, node+1)]) + nx.set_node_attributes(meta_mol, from_itp, "from_itp") + # map to molecule + with pytest.raises(IOError): + polyply.src.map_to_molecule.MapToMolecule(ff).run_molecule(meta_mol) From 09df05a1d513e8ed378a251c943e9c9a8f4dc739 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Mon, 1 Aug 2022 14:41:21 +0200 Subject: [PATCH 3/9] fix multi-block for single residue --- polyply/src/map_to_molecule.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/polyply/src/map_to_molecule.py b/polyply/src/map_to_molecule.py index 74bd467c..878a5e45 100644 --- a/polyply/src/map_to_molecule.py +++ b/polyply/src/map_to_molecule.py @@ -129,6 +129,7 @@ def match_nodes_to_blocks(self, meta_molecule): regular_graph = nx.Graph() restart_graph = nx.Graph() restart_attr = nx.get_node_attributes(meta_molecule, "from_itp") + restart_graph.add_nodes_from(restart_attr.keys()) # in case we only have a single residue # we assume the user is sane and that residue is not from @@ -146,6 +147,7 @@ def match_nodes_to_blocks(self, meta_molecule): else: regular_graph.add_edge(idx, jdx) + print(restart_graph.nodes) # regular nodes have to match a block in the force-field by resname for node in regular_graph.nodes: self.node_to_block[node] = meta_molecule.nodes[node]["resname"] @@ -299,6 +301,7 @@ def run_molecule(self, meta_molecule): # a block which consists of multiple residues. This # gets entangled here self.match_nodes_to_blocks(meta_molecule) + print(self.node_to_block) # raise an error if a block is not known to the library _assert_blocks_in_FF(self.node_to_block.values(), self.force_field) From 49515e3025d0ddeb0848f81aab841f7bc7e07c30 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Mon, 1 Aug 2022 14:48:57 +0200 Subject: [PATCH 4/9] remove prints --- polyply/src/map_to_molecule.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/polyply/src/map_to_molecule.py b/polyply/src/map_to_molecule.py index 878a5e45..8972a43d 100644 --- a/polyply/src/map_to_molecule.py +++ b/polyply/src/map_to_molecule.py @@ -147,7 +147,6 @@ def match_nodes_to_blocks(self, meta_molecule): else: regular_graph.add_edge(idx, jdx) - print(restart_graph.nodes) # regular nodes have to match a block in the force-field by resname for node in regular_graph.nodes: self.node_to_block[node] = meta_molecule.nodes[node]["resname"] @@ -301,7 +300,6 @@ def run_molecule(self, meta_molecule): # a block which consists of multiple residues. This # gets entangled here self.match_nodes_to_blocks(meta_molecule) - print(self.node_to_block) # raise an error if a block is not known to the library _assert_blocks_in_FF(self.node_to_block.values(), self.force_field) From 8a253e88b379a17aa2fd3ff07147cb6cd4b7203e Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Mon, 1 Aug 2022 14:49:16 +0200 Subject: [PATCH 5/9] add test for two single residue mol multiblocks --- polyply/tests/test_map_to_molecule.py | 32 ++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/polyply/tests/test_map_to_molecule.py b/polyply/tests/test_map_to_molecule.py index 0a3967d7..aa62a500 100644 --- a/polyply/tests/test_map_to_molecule.py +++ b/polyply/tests/test_map_to_molecule.py @@ -174,7 +174,37 @@ Interaction(atoms=(6, 7), parameters=['1', '0.37', '8000'], meta={})], [(0, 1), (1, 2), (2, 3), (4, 5), (5, 6), (6, 7)], 8, - {0: "MIX", 1: "MIX", 2: "MIX", 3:"MIX"}) + {0: "MIX", 1: "MIX", 2: "MIX", 3:"MIX"}), + # test two consecutive multiblock fragments with single residues + # and different mol-name + (""" + [ moleculetype ] + ; name nexcl. + OTR 1 + ; + [ atoms ] + 1 SN1a 1 R2 C1 1 0.000 45 + 2 SN1a 1 R2 C2 1 0.000 45 + [ bonds ] + ; back bone bonds + 1 2 1 0.4 7000 + [ moleculetype ] + ; name nexcl. + MIX 1 + ; + [ atoms ] + 1 SN1a 1 R1 C1 1 0.000 45 + 2 SN1a 1 R1 C2 1 0.000 45 + [ bonds ] + ; back bone bonds + 1 2 1 0.37 7000 + """, + ["R1", "R2"], + [Interaction(atoms=(0, 1), parameters=['1', '0.37', '7000'], meta={}), + Interaction(atoms=(2, 3), parameters=['1', '0.4', '7000'], meta={}),], + [(0, 1), (2, 3),], + 4, + {0: "MIX", 1: "OTR"}) )) def test_multiresidue_block(lines, monomers, bonds, edges, nnodes, from_itp): """ From f285353c0120454bb2f62bc3dca9778d237012a9 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Tue, 2 Aug 2022 13:34:03 +0200 Subject: [PATCH 6/9] simulation runner for polyply topologies --- polyply/src/gmx_wrapper/workflow.py | 57 +++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 polyply/src/gmx_wrapper/workflow.py diff --git a/polyply/src/gmx_wrapper/workflow.py b/polyply/src/gmx_wrapper/workflow.py new file mode 100644 index 00000000..d13ccf46 --- /dev/null +++ b/polyply/src/gmx_wrapper/workflow.py @@ -0,0 +1,57 @@ +""" +Useful classes for integrating workflows with gromacs wrapper. +""" +import os +from pathlib import Path +import gromacs as gmx +from gromacs.run import MDrunner +import vermouth + +class TopologyGMXRunner(MDrunner): + """ + Run energy minization on a polyply topology. + To use this method instantiate the class and + then run the simulation with run() or run_check(). + The latter returns a success status of the simulation. + After the simulation the coordinates in the topolgy + object will be updated in place. + """ + def __init__(self, + topology, + toppath, + mdppath, + dirname=os.path.curdir, + **kwargs): + """ + Set some enviroment variables for the run. + """ + self.topology = topology + self.toppath = toppath + self.mdppath = mdppath + self.startpath = dirname / Path("start.gro") + self.outpath = dirname / Path(kwargs.get('deffnm', 'confout') + ".gro") + self.tprname = dirname / Path(kwargs.get('deffnm', 'topol') + ".tpr") + kwargs['s'] = self.tprname + super().__init__(dirname=dirname, **kwargs) + + def prehook(self, **kwargs): + """ + Write the coordinates of the current system and grompp + the input parameters. + """ + # write the system coordinates to file + system = self.topology.convert_to_vermouth_system() + vermouth.gmx.gro.write_gro(system, self.startpath, precision=7, + title="starting coordinates", + box=self.topology.box, defer_writing=False) + # grompp everything in the directory + gmx.grompp(f=self.mdppath, + c=self.startpath, + p=self.toppath, + o=self.tprname) + + def posthook(self, **kwargs): + """ + Read in the new coordinates from the completed run. + """ + self.topology.add_positions_from_file(self.outpath) From dbb6b9e3bb6883264de725f3b7053fdc7b19903a Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Tue, 2 Aug 2022 13:57:22 +0200 Subject: [PATCH 7/9] basic workflow for topologies --- bin/polyply | 1 + polyply/src/gen_coords.py | 3 +++ polyply/src/gmx_wrapper/workflow.py | 3 ++- 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/bin/polyply b/bin/polyply index a9a8dedd..0220ab08 100755 --- a/bin/polyply +++ b/bin/polyply @@ -69,6 +69,7 @@ def main(): # pylint: disable=too-many-locals,too-many-statements help='Input file (ITP|FF)', nargs="*") file_group.add_argument('-o', dest='outpath', type=Path, help='Output ITP (ITP)') + seq_group = file_group.add_mutually_exclusive_group(required=True) seq_group.add_argument('-seq', dest='seq', type=str, nargs='+', help='A linear sequence of residue names.') diff --git a/polyply/src/gen_coords.py b/polyply/src/gen_coords.py index ae33183a..22c6240c 100644 --- a/polyply/src/gen_coords.py +++ b/polyply/src/gen_coords.py @@ -28,6 +28,7 @@ from .annotate_ligands import AnnotateLigands, parse_residue_spec, _find_nodes from .build_file_parser import read_build_file from .check_residue_equivalence import check_residue_equivalence +from .gmx_wrapper.workflow import TopologyGMXRunner LOGGER = StyleAdapter(get_logger(__name__)) @@ -130,6 +131,8 @@ def gen_coords(toppath, Path to topology file outpath: :class:pathlib.Path Path to coordinate file + mdppath: :class:pathlib.Path + Path to mdp file default None name: str Name of the molecule build: :class:pathlib.Path diff --git a/polyply/src/gmx_wrapper/workflow.py b/polyply/src/gmx_wrapper/workflow.py index d13ccf46..edd0a2de 100644 --- a/polyply/src/gmx_wrapper/workflow.py +++ b/polyply/src/gmx_wrapper/workflow.py @@ -3,9 +3,10 @@ """ import os from pathlib import Path +import vermouth import gromacs as gmx from gromacs.run import MDrunner -import vermouth + class TopologyGMXRunner(MDrunner): """ From 18605665a719260addf43af771e34d809aace8ad Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Tue, 2 Aug 2022 15:23:07 +0200 Subject: [PATCH 8/9] workflow for replica generato --- bin/polyply | 26 ++++++++++-- polyply/__init__.py | 1 + polyply/src/gen_coords.py | 1 - polyply/src/gen_replicas.py | 65 +++++++++++++++++++++++++++++ polyply/src/gmx_wrapper/workflow.py | 36 ++++++++++------ 5 files changed, 112 insertions(+), 17 deletions(-) create mode 100644 polyply/src/gen_replicas.py diff --git a/bin/polyply b/bin/polyply index 0220ab08..8e8682d4 100755 --- a/bin/polyply +++ b/bin/polyply @@ -23,7 +23,7 @@ import argparse from pathlib import Path import numpy as np import polyply -from polyply import (gen_itp, gen_coords, gen_seq, DATA_PATH) +from polyply import (gen_itp, gen_coords, gen_seq, gen_replicas, DATA_PATH) from polyply.src.load_library import load_library from polyply.src.logging import LOGGER, LOGLEVELS @@ -49,7 +49,9 @@ def main(): # pylint: disable=too-many-locals,too-many-statements # List of Subparsers for the different tools parser_gen_itp = subparsers.add_parser('gen_params', aliases=['gen_itp']) - parser_gen_coords = subparsers.add_parser('gen_coords') + # Make this a parent parser so we can later wrap replica around it and + # get all the arguments we need + parser_gen_coords = argparse.ArgumentParser(add_help=False) parser_gen_seq = subparsers.add_parser('gen_seq') # ============================================================================= @@ -162,7 +164,25 @@ def main(): # pylint: disable=too-many-locals,too-many-statements backmap_group.add_argument('-back_fudge', dest='bfudge', type=float, default=0.4, help='factor by which to scale the coordinates when backmapping.') - parser_gen_coords.set_defaults(func=gen_coords) + # gen_coords specific arguments + subparser_gen_coords = subparsers.add_parser('gen_coords', parents=[parser_gen_coords]) + subparser_gen_coords.set_defaults(func=gen_coords) + + # gen_replica specific arguments + subparser_gen_replicas = subparsers.add_parser('gen_replicas', parents=[parser_gen_coords]) + replica_group = subparser_gen_replicas.add_argument_group('Options for replica generation') + replica_group.add_argument("-nrepl", dest='nrepl', type=int, + help="number of system replicas to generate") + replica_group.add_argument("-mdp", dest='mdppath', type=Path, default=Path("md.mdp"), + help="mdp file (.mdp)") + replica_group.add_argument("-log", dest='logpath', type=Path, default=Path("polyply.log"), + help="logfile name") + replica_group.add_argument("-wdir", dest='workdir', type=Path, default=os.path.curdir, + help="logfile name") + replica_group.add_argument("-timeout", dest='timeout', type=float, default=None, + help="time after which to timeout gen_coords") + + subparser_gen_replicas.set_defaults(func=gen_replicas) # ============================================================================ # Input Arguments for the sequence generation tool diff --git a/polyply/__init__.py b/polyply/__init__.py index 9ef7d384..e6db655d 100644 --- a/polyply/__init__.py +++ b/polyply/__init__.py @@ -49,3 +49,4 @@ from .src.gen_itp import gen_itp, gen_params from .src.gen_coords import gen_coords from .src.gen_seq import gen_seq +from .src.gen_replicas import gen_replicas diff --git a/polyply/src/gen_coords.py b/polyply/src/gen_coords.py index 22c6240c..513e234c 100644 --- a/polyply/src/gen_coords.py +++ b/polyply/src/gen_coords.py @@ -28,7 +28,6 @@ from .annotate_ligands import AnnotateLigands, parse_residue_spec, _find_nodes from .build_file_parser import read_build_file from .check_residue_equivalence import check_residue_equivalence -from .gmx_wrapper.workflow import TopologyGMXRunner LOGGER = StyleAdapter(get_logger(__name__)) diff --git a/polyply/src/gen_replicas.py b/polyply/src/gen_replicas.py new file mode 100644 index 00000000..74391eac --- /dev/null +++ b/polyply/src/gen_replicas.py @@ -0,0 +1,65 @@ +""" +Generate multiple system replicas using gen_coords. In essence +this is a wrapper for gen_coords which also makes use of the +gromacs wrapper libary. +""" +import os +from pathlib import Path +from vermouth.log_helpers import StyleAdapter, get_logger +from polyply import gen_coords +from polyply.src.gmx_wrapper.workflow import GMXRunner + +LOGGER = StyleAdapter(get_logger(__name__)) + +def gen_replicas(nrepl, + mdppath, + logpath, + workdir, + timeout, + **kwargs): + """ + Generate multiple coordinate replicas of a system. + this function is a wrapper around gen_coords, which + also has some options for running simulations. + + Paremeters + ---------- + nrepl: int + number of replicas to generate + mdppath: Path + file path to mdp-file + logpath: Path + file path to log-file + workdir: Path + workdir to use + timeout: float + time after which gen_coords is timed out + """ + mdppath = mdppath.resolve() + workdir = workdir.resolve() + logpath = logpath.resolve() + kwargs["toppath"] = kwargs["toppath"].resolve() + base_name = kwargs['outpath'].name + os.chdir(workdir) + for idx in range(0, nrepl): + os.mkdir(workdir / Path(str(idx))) + os.chdir(workdir / Path(str(idx))) + kwargs['outpath'] = workdir / Path(str(idx)) / Path(base_name) + gen_coords(**kwargs) + # Running energy minization + LOGGER.info("running energy minization") + sim_runner = GMXRunner(kwargs["outpath"], + kwargs["toppath"], + mdppath, + deffnm="min") + status = sim_runner.run() + # everything has worked out yay! + with open(logpath, "a") as logfile: + if status != 0: + logfile.write(f"replica {idx} has failed with unkown gmx issue\n") + if not sim_runner.success: + logfile.write(f"replica {idx} has failed with infinite forces\n") + + logfile.write(f"replica {idx} was successfully generated\n") + + os.chdir(workdir) diff --git a/polyply/src/gmx_wrapper/workflow.py b/polyply/src/gmx_wrapper/workflow.py index edd0a2de..e0f72d62 100644 --- a/polyply/src/gmx_wrapper/workflow.py +++ b/polyply/src/gmx_wrapper/workflow.py @@ -3,12 +3,12 @@ """ import os from pathlib import Path -import vermouth import gromacs as gmx from gromacs.run import MDrunner +import gromacs.environment +gromacs.environment.flags['capture_output'] = "file" - -class TopologyGMXRunner(MDrunner): +class GMXRunner(MDrunner): """ Run energy minization on a polyply topology. To use this method instantiate the class and @@ -18,7 +18,7 @@ class TopologyGMXRunner(MDrunner): object will be updated in place. """ def __init__(self, - topology, + coordpath, toppath, mdppath, dirname=os.path.curdir, @@ -26,12 +26,14 @@ def __init__(self, """ Set some enviroment variables for the run. """ - self.topology = topology + self.coordpath = coordpath self.toppath = toppath self.mdppath = mdppath + self.success = False self.startpath = dirname / Path("start.gro") self.outpath = dirname / Path(kwargs.get('deffnm', 'confout') + ".gro") self.tprname = dirname / Path(kwargs.get('deffnm', 'topol') + ".tpr") + self.deffnm = kwargs.get('deffnm', None) kwargs['s'] = self.tprname super().__init__(dirname=dirname, **kwargs) @@ -40,19 +42,27 @@ def prehook(self, **kwargs): Write the coordinates of the current system and grompp the input parameters. """ - # write the system coordinates to file - system = self.topology.convert_to_vermouth_system() - vermouth.gmx.gro.write_gro(system, self.startpath, precision=7, - title="starting coordinates", - box=self.topology.box, defer_writing=False) # grompp everything in the directory gmx.grompp(f=self.mdppath, - c=self.startpath, + c=self.coordpath, p=self.toppath, o=self.tprname) def posthook(self, **kwargs): """ - Read in the new coordinates from the completed run. + Make sure that the energy minization did not result into + infinite forces. """ - self.topology.add_positions_from_file(self.outpath) + if self.deffnm: + logpath = Path(self.deffnm + ".log") + else: + logpath = Path("md.log") + + with open(logpath, "r") as logfile: + for line in logfile.readlines(): + if "Norm" in line: + maxforce = line.strip().split()[-1] + if maxforce == "inf": + self.success = False + else: + self.success = True From 9db2724ad8a2cb4e5a930c20d9e1a1617eb328d5 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Tue, 2 Aug 2022 15:58:49 +0200 Subject: [PATCH 9/9] allow for mdp defined workflow --- bin/polyply | 4 +-- polyply/src/gen_replicas.py | 56 ++++++++++++++++++++++++------------- 2 files changed, 38 insertions(+), 22 deletions(-) diff --git a/bin/polyply b/bin/polyply index 8e8682d4..2a795530 100755 --- a/bin/polyply +++ b/bin/polyply @@ -173,8 +173,8 @@ def main(): # pylint: disable=too-many-locals,too-many-statements replica_group = subparser_gen_replicas.add_argument_group('Options for replica generation') replica_group.add_argument("-nrepl", dest='nrepl', type=int, help="number of system replicas to generate") - replica_group.add_argument("-mdp", dest='mdppath', type=Path, default=Path("md.mdp"), - help="mdp file (.mdp)") + replica_group.add_argument("-mdp", dest='mdppaths', type=Path, default=[Path("md.mdp")], + help="mdp file (.mdp)", nargs="*") replica_group.add_argument("-log", dest='logpath', type=Path, default=Path("polyply.log"), help="logfile name") replica_group.add_argument("-wdir", dest='workdir', type=Path, default=os.path.curdir, diff --git a/polyply/src/gen_replicas.py b/polyply/src/gen_replicas.py index 74391eac..0c142563 100644 --- a/polyply/src/gen_replicas.py +++ b/polyply/src/gen_replicas.py @@ -12,7 +12,7 @@ LOGGER = StyleAdapter(get_logger(__name__)) def gen_replicas(nrepl, - mdppath, + mdppaths, logpath, workdir, timeout, @@ -26,7 +26,7 @@ def gen_replicas(nrepl, ---------- nrepl: int number of replicas to generate - mdppath: Path + mdppaths: list[Path] file path to mdp-file logpath: Path file path to log-file @@ -35,31 +35,47 @@ def gen_replicas(nrepl, timeout: float time after which gen_coords is timed out """ - mdppath = mdppath.resolve() + mdps_resolved = [] + for mdppath in mdppaths: + mdps_resolved.append(mdppath.resolve()) + workdir = workdir.resolve() logpath = logpath.resolve() kwargs["toppath"] = kwargs["toppath"].resolve() base_name = kwargs['outpath'].name os.chdir(workdir) + for idx in range(0, nrepl): - os.mkdir(workdir / Path(str(idx))) - os.chdir(workdir / Path(str(idx))) - kwargs['outpath'] = workdir / Path(str(idx)) / Path(base_name) + replica_dir = workdir / Path(str(idx)) + os.mkdir(replica_dir) + os.chdir(replica_dir) + kwargs['outpath'] = replica_dir / Path(base_name) gen_coords(**kwargs) - # Running energy minization - LOGGER.info("running energy minization") - sim_runner = GMXRunner(kwargs["outpath"], - kwargs["toppath"], - mdppath, - deffnm="min") - status = sim_runner.run() - # everything has worked out yay! - with open(logpath, "a") as logfile: - if status != 0: - logfile.write(f"replica {idx} has failed with unkown gmx issue\n") - if not sim_runner.success: - logfile.write(f"replica {idx} has failed with infinite forces\n") - logfile.write(f"replica {idx} was successfully generated\n") + for jdx, mdppath in enumerate(mdps_resolved): + mdpname = str(mdppath.stem) + deffnm = str(jdx) + "_" + mdpname + # Get previous coordinates + if jdx == 0: + startpath = kwargs["outpath"] + else: + startpath = replica_dir / Path(prev_deffnm + ".gro") + + # Running GROMACS protocol + LOGGER.info(f"running simulation {mdpname}") + sim_runner = GMXRunner(startpath, + kwargs["toppath"], + mdppath, + deffnm=deffnm) + status = sim_runner.run() + prev_deffnm = deffnm + + with open(logpath, "a") as logfile: + if status != 0: + logfile.write(f"replica {idx} step {jdx} has failed with unkown gmx issue\n") + if not sim_runner.success: + logfile.write(f"replica {idx} step {jdx} has failed with infinite forces\n") + + logfile.write(f"replica {idx} step {jdx} successfully completed\n") os.chdir(workdir)