move header writing in gmx files to a system level meta attribute

bin/martinize2

@@ -912,19 +912,6 @@ def entry():
         write_canon=args.write_canon,
     )
 
-    # grompp has a limit in the number of character it can read per line
-    # (due to the size limit of a buffer somewhere in its implementation).
-    # The command line can be longer than this limit and therefore
-    # prevent grompp from reading the topology.
-    gromacs_char_limit = 4000  # the limit is actually 4095, but I play safe
-    command = " ".join(sys.argv)
-    if len(command) > gromacs_char_limit:
-        command = command[:gromacs_char_limit] + " ..."
-    system.meta["header"].extend(("This file was generated using the following command:",
-                                  command,
-                                  VERSION,
-                                 ))
-
     LOGGER.info("Read input.", type="step")
     for molecule in system.molecules:
         LOGGER.debug("Read molecule {}.", molecule, type="step")
@@ -1156,7 +1143,9 @@ def entry():
                            args.top_path,
                            itp_paths=itp_paths,
                            C6C12=False,
-                           defines=defines)
+                           defines=defines,
+                           command=sys.argv,
+                           version=VERSION)
 
     # Write a PDB file.
     vermouth.pdb.write_pdb(system, str(args.outpath), omit_charges=True)

vermouth/dssp/dssp.py

@@ -554,7 +554,7 @@ def gmx_system_header(system):
         )
     )
 
-    if None not in ss_sequence:
+    if (None not in ss_sequence) and (len(ss_sequence) > 0):
         system.meta["header"].extend(("The following sequence of secondary structure ",
                                       "was used for the full system:",
                                       "".join(ss_sequence),

vermouth/gmx/topology.py

@@ -118,7 +118,9 @@
                        itp_paths={"nonbond_params": "extra_nbparams.itp",
                                   "atomtypes": "extra_atomtypes.itp"},
                        C6C12=False,
-                       defines=()):
+                       defines=(),
+                       command=[],
+                       version=''):
     """
     Writes a Gromacs .top file for the specified system. Gromacs topology
     files are defined by directives for example `[ atomtypes ]`. However,
@@ -143,6 +145,10 @@
         C6C12 form
     defines: tuple(str)
         define statments to include in the topology
+    command: list
+        command used to generate the system
+    version: str
+        version of vermouth used to generate system
     """
     if not system.molecules:
         raise ValueError("No molecule in the system. Nothing to write.")
@@ -178,6 +184,19 @@
         system.molecules, key=lambda x: x.meta["moltype"]
     )
 
+    # grompp has a limit in the number of character it can read per line
+    # (due to the size limit of a buffer somewhere in its implementation).
+    # The command line can be longer than this limit and therefore
+    # prevent grompp from reading the topology.
+    gromacs_char_limit = 4000  # the limit is actually 4095, but I play safe
+    command_used = " ".join(command)
+    if len(command_used) > gromacs_char_limit:
+        command_used = command_used[:gromacs_char_limit] + " ..."
+    system.meta["header"].extend(("This file was generated using the following command:",
+                                  command_used,
+                                  version,
+                                  ))
+
     header = system.meta.get("header", [])
 
     for moltype, molecules in molecule_groups:

vermouth/tests/test_dssp.py

@@ -707,19 +707,30 @@ def test_convert_dssp_to_martini(sequence, expected):
     found = dssp.convert_dssp_to_martini(sequence)
     assert expected == found
 
-
-def test_dssp_system_header(test_molecule):
+@pytest.mark.parametrize('resnames, ss_string, secstruc',
+     (      # protein resnames with secstruc
+             ({0: "ALA", 1: "ALA", 2: "ALA",
+               3: "GLY", 4: "GLY",
+               5: "MET",
+               6: "ARG", 7: "ARG", 8: "ARG"},
+              'HHHH',
+              {1: "H", 2: "H", 3: "H", 4: "H"}
+              ),
+             # not protein resnames, no secstruc annotated or gets written
+             ({0: "A", 1: "A", 2: "A",
+               3: "B", 4: "B",
+               5: "C",
+               6: "D", 7: "D", 8: "D"},
+              "",
+              {1: "", 2: "", 3: "", 4: ""}
+              )
+     ))
+def test_gmx_system_header(test_molecule, resnames, ss_string, secstruc):
 
     atypes = {0: "P1", 1: "SN4a", 2: "SN4a",
               3: "SP1", 4: "C1",
               5: "TP1",
               6: "P1", 7: "SN3a", 8: "SP4"}
-    # the molecule resnames
-    resnames = {0: "A", 1: "A", 2: "A",
-                3: "B", 4: "B",
-                5: "C",
-                6: "D", 7: "D", 8: "D"}
-    secstruc = {1: "H", 2: "H", 3: "H", 4: "H"}
 
     system = create_sys_all_attrs(test_molecule,
                                   moltype="molecule_0",
@@ -728,6 +739,10 @@ def test_dssp_system_header(test_molecule):
                                   attrs={"resname": resnames,
                                          "atype": atypes})
 
+    # annotate the actual 'secstruct' attribute because create_sys_all_attrs actually annotates cgsecstruct
+    dssp.AnnotateResidues(attribute="secstruct",
+                          sequence="HHHH").run_system(system)
+
     dssp.AnnotateMartiniSecondaryStructures().run_system(system)
 
-    assert system.meta.get('header')
+    assert ss_string in system.meta.get('header', [''])