Merge pull request #643 from marrink-lab/Go-contact-map

internal go contact map

bin/martinize2

@@ -48,7 +48,7 @@ from vermouth.map_input import (
     generate_all_self_mappings,
     combine_mappings,
 )
-from vermouth.rcsu.contact_map import read_go_map
+from vermouth.rcsu.contact_map import GenerateContactMap, read_go_map
 from vermouth.rcsu.go_pipeline import GoPipeline
 from vermouth.gmx.topology import write_gmx_topology
 
@@ -546,11 +546,11 @@ def entry():
     go_group.add_argument(
         "-go",
         dest="go",
-        required=False,
+        nargs='?',
+        const=True,
         type=Path,
-        default=None,
-        help="Contact map to be used for the Martini Go model."
-             "Currently, only one format is supported. See docs.",
+        help="Use Martini Go model. Accepts either an input file from the server, "
+             "or just provide the flag to calculate as part of Martinize."
     )
     go_group.add_argument(
         "-go-eps",
@@ -581,6 +581,29 @@ def entry():
         help=("Minimum graph distance (similar sequence distance) below which"
               "contacts are removed. "),
     )
+    go_group.add_argument(
+        "-go-write-file",
+        dest="go_write_file",
+        nargs='?',
+        const='contact_map_martinize.out',
+        default=False,
+        type=Path,
+        help=("Write out contact map to file if calculating as part of Martinize2.")
+    )
+    go_group.add_argument(
+        "-go-backbone",
+        dest="go_backbone",
+        default="BB",
+        type=str,
+        help=("Name of protein backbone bead on which to place a virtual interaction go site")
+    )
+    go_group.add_argument(
+        "-go-atomname",
+        dest="go_atomname",
+        default="CA",
+        type=str,
+        help=("Name of the virtual interaction go site atom")
+    )
 
     water_group = parser.add_argument_group("Apply water bias.")
     water_group.add_argument(
@@ -967,6 +990,17 @@ def entry():
     elif args.cystein_bridge != "auto":
         vermouth.AddCysteinBridgesThreshold(args.cystein_bridge).run_system(system)
 
+    if args.go:
+        system = vermouth.MergeAllMolecules().run_system(system)
+        # need this here because have to get contact map at atomistic resolution
+        if isinstance(args.go, Path):
+            LOGGER.info("Reading Go model contact map.", type="step")
+            read_go_map(system=system, file_path=args.go)
+        else:
+            LOGGER.info("Generating Go model contact map.", type="step")
+            GenerateContactMap(write_file=args.go_write_file).run_system(system)
+
+
     # Run martinize on the system.
     system = martinize(
         system,
@@ -988,18 +1022,17 @@ def entry():
         vermouth.ApplyPosres(node_selector, args.posres_fc).run_system(system)
 
     # Generate the Go model if required
-    if args.go:
+    if system.go_params["go_map"]:
         go_name_prefix = args.molname
-        LOGGER.info("Reading Go model contact map.", type="step")
-        go_map = read_go_map(args.go)
         LOGGER.info("Generating the Go model.", type="step")
         GoPipeline.run_system(system,
                               moltype=go_name_prefix,
-                              contact_map=go_map,
                               cutoff_short=args.go_low,
                               cutoff_long=args.go_up,
                               go_eps=args.go_eps,
-                              res_dist=args.go_res_dist,)
+                              res_dist=args.go_res_dist,
+                              go_anchor_bead=args.go_backbone,
+                              go_atomname=args.go_atomname)
 
         defines = ("GO_VIRT",)
         itp_paths = {"atomtypes": "go_atomtypes.itp",

doc/source/tutorials/go_models.rst

@@ -14,28 +14,49 @@ the form of Lennard-Jones interactions, which are written as an extra file to be
 
 The Gō model is described in the help::
 
- Virtual site based GoMartini:
-   -go GO                Contact map to be used for the Martini Go model. Currently, only one format is supported. See docs. (default: None)
+ Virtual site based GōMartini:
+   -go [GO]              Use Martini Go model. Accepts either an input file from the server, or just provide the flag to
+                         calculate as part of Martinize. (default: None)
    -go-eps GO_EPS        The strength of the Go model structural bias in kJ/mol. (default: 9.414)
-   -go-moltype GOVS_MOLTYPE
-                         Set the name of the molecule when using Virtual Sites GoMartini. (default: molecule_0)
    -go-low GO_LOW        Minimum distance (nm) below which contacts are removed. (default: 0.3)
    -go-up GO_UP          Maximum distance (nm) above which contacts are removed. (default: 1.1)
    -go-res-dist GO_RES_DIST
                          Minimum graph distance (similar sequence distance) below which contacts are removed. (default: 3)
+   -go-write-file [GO_WRITE_FILE]
+                         Write out contact map to file if calculating as part of Martinize2. (default: None)
+   -go-backbone GO_BACKBONE
+                         Name of protein backbone bead on which to place a virtual interaction go site (default: BB)
+   -go-atomname GO_ATOMNAME
+                         Name of the virtual interaction go site atom (default: CA)
+
+To add a Gō model to your protein, the first step is to calculate the contact map of your protein.
+The contact map can be obtained in two ways. Firstly, by uploading your pdb structure
+to the `web server <http://pomalab.ippt.pan.pl/GoContactMap/>`_, and downloading the associated ``contact_map.out`` file.
+Alternatively, the contact map can be calculated directly without the need for
+any external processes. While the implementations of the contact algorithm are identical, the Martinize2 implementation
+may be relatively slow for larger systems. Typically for proteins with fewer than 1000 residues, the calculation of the
+contact map as part of Martinize2 will add up to a minute of extra calculation. Note that while the implementations of
+the main algorithm are identical, there may be small differences in the resulting contact map files due to assumptions
+the server makes about the format of input pdb files, which the implementation in Martinize2 overcomes. If you want
+to check the contact map that Martinize2 has calculated, you can write it out using the ``-go-write-file`` argument.
+While the contact map files may have small differences, it is likely that they will still result in the same non-bonded
+file outputs, a result of how symmetrical contacts are further identified in the definition of the Gō model.
+
+The Gō model is added to the input system using the ``-go`` argument of martinize2. If you have used a contact map
+from the server, give the path to the contact map file as the argument:
 
-To add a Gō model to your protein, the first step is to calculate the contact map of your protein by uploading it
-to the `web server <http://pomalab.ippt.pan.pl/GoContactMap/>`_.
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go contact_map.out``
 
-The contact map is then used in your martinize2 command:
+Otherwise the contact map is calculated as part of Martinize2 by just specifying the ``-go`` argument:
 
-``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go contact_map.out``
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go``
 
 
 Without any further additions, this will:
- 1) Generate virtual sites with the atomname ``CA`` directly on top of all the backbone beads in your protein.
-    Each ``CA`` atom has an underlying atomtype (see below), which has a default name specified using the
-    ``-go-moltype`` flag.
+ 1) Generate virtual sites with the atomname ``CA`` directly on top of all the backbone beads ``BB`` in your protein.
+    The atomname and the name of the backbone bead are controlled using the ``-go-atomname`` and ``-go-backbone`` flags
+    respectively. Each ``CA`` atom has an underlying atomtype (see below), which has a default name inherited from the
+    ``-name`` flag.
  2) Use the contact map to generate a set of non-bonded parameters between specific pairs of ``CA`` atoms in your molecule
     with strength 9.414 kJ/mol (changed through the ``-go-eps`` flag).
  3) Eliminate any contacts which are shorter than 0.3 nm and longer than 1.1 nm, or are closer than 3 residues in the
@@ -51,23 +72,23 @@ interactions between them.
 For example, ``go_atomtypes.itp`` looks like any other ``[ atomtypes ]`` directive::
 
  [ atomtypes ]
- molecule_0_1 0.0 0 A 0.00000000 0.00000000
- molecule_0_2 0.0 0 A 0.00000000 0.00000000
- molecule_0_3 0.0 0 A 0.00000000 0.00000000
- molecule_0_4 0.0 0 A 0.00000000 0.00000000
- molecule_0_5 0.0 0 A 0.00000000 0.00000000
+ molecule_1 0.0 0 A 0.00000000 0.00000000
+ molecule_2 0.0 0 A 0.00000000 0.00000000
+ molecule_3 0.0 0 A 0.00000000 0.00000000
+ molecule_4 0.0 0 A 0.00000000 0.00000000
+ molecule_5 0.0 0 A 0.00000000 0.00000000
  ...
 
 Similarly, ``go_nbparams.itp`` looks like any ``[ nonbond_params ]`` directive (obviously, the exact parameters here
 depend on your protein)::
 
  [ nonbond_params ]
- molecule_0_17 molecule_0_13 1 0.59354169 9.41400000 ;go bond 0.666228018941817
- molecule_0_18 molecule_0_14 1 0.53798937 9.41400000 ;go bond 0.6038726468003999
- molecule_0_19 molecule_0_15 1 0.51270658 9.41400000 ;go bond 0.5754936778307316
- molecule_0_22 molecule_0_15 1 0.73815666 9.41400000 ;go bond 0.8285528398039018
- molecule_0_22 molecule_0_18 1 0.54218134 9.41400000 ;go bond 0.6085779754055839
- molecule_0_23 molecule_0_19 1 0.53307395 9.41400000 ;go bond 0.5983552758317587
+ molecule_17 molecule_0_13 1 0.59354169 9.41400000 ;go bond 0.666228018941817
+ molecule_18 molecule_0_14 1 0.53798937 9.41400000 ;go bond 0.6038726468003999
+ molecule_19 molecule_0_15 1 0.51270658 9.41400000 ;go bond 0.5754936778307316
+ molecule_22 molecule_0_15 1 0.73815666 9.41400000 ;go bond 0.8285528398039018
+ molecule_22 molecule_0_18 1 0.54218134 9.41400000 ;go bond 0.6085779754055839
+ molecule_23 molecule_0_19 1 0.53307395 9.41400000 ;go bond 0.5983552758317587
  ...
 
 To activate your Gō model for use in Gromacs, the `martini_v3.0.0.itp` master itp needs the additional files included.
@@ -125,7 +146,36 @@ not truly as intermolecular forces in the system as a whole. For more detail on
 see the paper by `Korshunova et al. <https://pubs.acs.org/doi/10.1021/acs.jctc.4c00677>`_.
 
 
-Visualising Go networks
+Multiple Gō models in the same system
+-------------------------------------
+
+If you have several different proteins that you have martinized with their own Gō models, then several extra steps are
+required to ensure that they can be simulated together.
+
+Firstly, when the proteins are coarse grained with Martinize2, the `-name` flag must be used to ensure that all the virtual
+sites created for the purposes of the Gō model have unique names and atomtypes. For example::
+
+ martinize2 -f protein.pdb -x cg.pdb -o topol.top -go contact_map.out -name my_protein
+
+will ensure that the atoms created to apply Gō sites to are names `my_protein_{0..n}` for a protein of n residues.
+
+Having martinized the proteins in this way, all the `go_atomtypes.itp` and `go_nbparams.itp` files generated for each
+protein should be concatenated into a single file, which may then be included in the master force field file as
+described above.
+
+Note on multiple Gō models
+--------------------------
+
+While the Gō model is expressed as a set of interactions between Gō sites on a protein, interactions are not
+generally extended over all copies of a protein. That is, if a simulation is set up with several copies of a
+multimeric protein where the monomers are held together by a Gō model, Gromacs will not permit the multimers
+to "fall apart" for monomers to find other monomers that were initially in other complexes. This limitation is
+a result of the nature of the Gromacs itp format, with the Gō interactions described within each `[ moleculetype ]`
+directive. This issue is discussed more extensively in a recent paper by `Korshunova et al. <https://pubs.acs.org/doi/10.1021/acs.jctc.4c00677>`_.
+
+
+
+Visualising Gō networks
 ----------------------------
 
 If you want to look at your Gō network in VMD to confirm that it's been constructed in the