External EN bonds itp

bin/martinize2

@@ -539,6 +539,13 @@ def entry():
             "format: <start_resid_1>:<end_resid_1>, <start_resid_2>:<end_resid_2>..."
         ),
     )
+    rb_group.add_argument(
+        "-ew",
+        dest="rb_write",
+        default=False,
+        action="store_true",
+        help="Write the elastic network to an external .itp file",
+    )
     go_group = parser.add_argument_group("Virtual site based GoMartini")
     go_group.add_argument(
         "-go",
@@ -1059,6 +1066,7 @@ def entry():
             selector=selector,
             domain_criterion=domain_criterion,
             res_min_dist=args.res_min_dist,
+            write_out=args.rb_write,
         )
         rubber_band_processor.run_system(system)
 

vermouth/gmx/itp.py

@@ -18,6 +18,7 @@
 
 import copy
 import itertools
+from vermouth.file_writer import deferred_open
 
 __all__ = ['write_molecule_itp', ]
 
@@ -54,6 +55,27 @@ def _interaction_sorting_key(interaction):
 
     return (conditional, group)
 
+def write_bond_params(interactions, itp_path, correspondence, max_length):
+    """
+    Writes the [ bonds ] directive to file.
+
+    interactions: list
+        list of bonded interactions to write to file
+    itp_path: str
+        name of external file to write
+    """
+    with deferred_open(itp_path, "w") as itp_file:
+        itp_file.write('[ bonds ]\n')
+        for b_params in interactions:
+            atoms = ['{atom_idx:>{max_length[idx]}}'
+                     .format(atom_idx=correspondence[x],
+                             max_length=max_length)
+                     for x in b_params.atoms]
+            parameters = ' '.join(str(x) for x in b_params.parameters)
+            comment = ''
+
+            to_join = atoms + [parameters]
+            itp_file.write(' '.join(to_join) + comment + '\n')
 
 def write_molecule_itp(molecule, outfile, header=(), moltype=None,
                        post_section_lines=None, pre_section_lines=None):
@@ -194,6 +216,15 @@ def write_molecule_itp(molecule, outfile, header=(), moltype=None,
         # should be written under the [ dihedrals ] section of the ITP file.
         if name == 'impropers':
             name = 'dihedrals'
+        # el_bonds is a directive for elastic network bonds to be written out externally
+        if name == 'el_bonds':
+            outfile.write('#ifdef ELASTIC\n')
+            outfile.write('#include "en_bonds.itp"\n')
+            outfile.write('#endif\n\n')
+            write_bond_params(interactions, 'en_bonds.itp',
+                              correspondence, max_length)
+            continue
+
         outfile.write('[ {} ]\n'.format(name))
         seen_sections.add(name)
         for line in pre_section_lines.get(name, []):

vermouth/gmx/topology.py

@@ -158,6 +158,8 @@ def write_gmx_topology(system,
         _path = itp_paths['nonbond_params']
         write_nonbond_params(system, _path, C6C12)
         include_string += f'\n #include "{_path}"\n'
+
+
     # Write the ITP files for the molecule types, and prepare writing the
     # [ molecules ] section of the top file.
     # * We write one ITP file for each different moltype in the system, the

vermouth/processors/apply_rubber_band.py

@@ -227,7 +227,7 @@ def apply_rubber_band(molecule, selector,
                       lower_bound, upper_bound,
                       decay_factor, decay_power,
                       base_constant, minimum_force,
-                      bond_type, domain_criterion, res_min_dist):
+                      bond_type, domain_criterion, res_min_dist, eb_write):
     r"""
     Adds a rubber band elastic network to a molecule.
 
@@ -295,6 +295,9 @@ def apply_rubber_band(molecule, selector,
         Minimum separation between two atoms for a bond to be kept.
         Bonds are kept is the separation is greater or equal to the value
         given.
+    eb_write: bool
+        if True, EN parameters get returned to be written out to an external file.
+        if False, EN parameters are written into the molecule's topology.
     """
     selection = []
     coordinates = []
@@ -341,6 +344,11 @@ def apply_rubber_band(molecule, selector,
     # Multiply the force constant by 0 if the nodes cannot be linked.
     constants *= can_be_linked
     distance_matrix = distance_matrix.round(5)  # For compatibility with legacy
+    ebs = []
+    if eb_write:
+        interaction_type = 'el_bonds'
+    else:
+        interaction_type = 'bonds'
     for from_idx, to_idx in zip(*np.triu_indices_from(constants)):
         # note the indices in the matrix are not anymore the idx of
         # the full molecule but the subset of nodes in selection
@@ -350,12 +358,12 @@ def apply_rubber_band(molecule, selector,
         length = distance_matrix[from_idx, to_idx]
         if force_constant > minimum_force:
             molecule.add_interaction(
-                type_='bonds',
+                type_=interaction_type,
                 atoms=(from_key, to_key),
                 parameters=[bond_type, length, force_constant],
-                meta={'group': 'Rubber band'},
+                meta={'group': 'Rubber band'}
             )
-
+    return ebs
 
 def always_true(*args, **kwargs):  # pylint: disable=unused-argument
     """
@@ -493,7 +501,8 @@ def __init__(self, lower_bound, upper_bound, decay_factor, decay_power,
                  selector=selectors.select_backbone,
                  bond_type_variable='elastic_network_bond_type',
                  res_min_dist_variable='elastic_network_res_min_dist',
-                 domain_criterion=always_true):
+                 domain_criterion=always_true,
+                 write_out=False):
         super().__init__()
         self.lower_bound = lower_bound
         self.upper_bound = upper_bound
@@ -507,6 +516,7 @@ def __init__(self, lower_bound, upper_bound, decay_factor, decay_power,
         self.domain_criterion = domain_criterion
         self.res_min_dist = res_min_dist
         self.res_min_dist_variable = res_min_dist_variable
+        self.write_out = write_out
 
     def run_molecule(self, molecule):
         # Choose the bond type. From high to low, the priority order is:
@@ -535,5 +545,19 @@ def run_molecule(self, molecule):
                           minimum_force=self.minimum_force,
                           bond_type=bond_type,
                           domain_criterion=self.domain_criterion,
-                          res_min_dist=res_min_dist)
+                          res_min_dist=res_min_dist,
+                          eb_write=self.write_out)
+
         return molecule
+
+    def run_system(self, system):
+        """
+        Process `system`.
+
+        Parameters
+        ----------
+        system: vermouth.system.System
+            The system to process. Is modified in-place.
+        """
+        self.system = system
+        super().run_system(system)