first draft of martinize2 command tutorials

doc/source/tutorials/basic_usage.rst

@@ -0,0 +1,89 @@
+===========
+Basic usage
+===========
+
+Overly basic usage
+==================
+
+Without any other additions, ``martinize2`` can take your protein, and make a ready coarse
+grained model with some martini
+forcefield:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb``
+
+This command will (try) and convert your protein from the atomistic input to one
+in the Martini3 force field.
+
+Other force fields are available to convert your protein to. To view them you
+can use ``martinize2 -list-ff``.
+
+Minimal physical reality usage
+==============================
+
+Our knowledge of proteins and Martini tells us that we need to add some more
+information to the topology to account for secondary structure.
+
+Let martinize2 deal with secondary structure for you
+----------------------------------------------------
+
+Martinize2 can deal with secondary structure intelligently using dssp in one of two ways:
+
+1) Use a dssp executable installed on your machine
+2) Use the dssp implementation in `mdtraj <https://mdtraj.org/1.9.4/api/generated/mdtraj.compute_dssp.html>`_
+
+To explain these more:
+
+1) If you have dssp installed locally, note that martinize2 is only validated for particular versions of dssp.
+   If a non-validated version is used, a warning will be raised and nothing is written.
+   If you know what you're doing and are confident with what's been produced, you can override the warning
+   with the ``-maxwarn`` flag. Otherwise, dssp can be used using the ``-dssp`` flag in martinize.
+
+   Where you have a local installation, replace ``/path/to/dssp`` in the following command with the
+   location of your installation:
+
+   ``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp /path/to/dssp``
+
+2) If you have `mdtraj <https://mdtraj.org/1.9.4/api/generated/mdtraj.compute_dssp.html>`_ installed in
+   your python environment, the ``-dssp`` flag can be left blank, and martinize2 will attempt to use
+   dssp from there:
+
+   ``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp``
+
+
+User knows best
+---------------
+
+If you already know the secondary structure of your protein and don't want to worry about
+dssp assigning it correctly, the ``-ss`` flag can be used instead.
+
+The ``-ss`` flag must be one of either:
+
+1)   The same length as the number of residues in your protein, with a dssp code for each residue
+2)   A single letter (eg. '``H``'), which will apply the same secondary structure parameters to your entire protein.
+
+For example:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ss HHHHHHHHHH``
+
+will read the ``HHHHHHHHHH`` dssp formatted string, indicating that there are exactly 10 residues in the
+input pdb which should all be treated as part of a helix. If the string provided does not contain the same
+number of residues as the input file, an error will be raised.
+
+Alternatively in this case, as we know everything should be a helix, the same result can be achieved through
+using a single letter as described above:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ss H``
+
+If instead we needed to assert that only the first five residues are in a helix, and the final five are coiled,
+we must we the full length string again:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ss HHHHHCCCCC``
+
+
+
+
+
+
+
+
+

doc/source/tutorials/elastic_networks.rst

@@ -0,0 +1,129 @@
+================
+Elastic Networks
+================
+
+The first method applied to maintain the secondary and tertiary structure
+of Martini proteins was `Elastic Networks <https://doi.org/10.1021/ct9002114>`_.
+
+Elastic networks are formed by finding contacts between protein backbone
+beads within a particular cutoff distance, and applying harmonic bonds between them,
+restraining them at the initial distance throughout your simulation.
+
+Elastic networks are applied in Martinize2 as per the section of the help::
+
+
+  -elastic              Write elastic bonds (default: False)
+  -ef RB_FORCE_CONSTANT
+                        Elastic bond force constant Fc in kJ/mol/nm^2 (default: 500)
+  -el RB_LOWER_BOUND    Elastic bond lower cutoff: F = Fc if rij < lo (default: 0)
+  -eu RB_UPPER_BOUND    Elastic bond upper cutoff: F = 0 if rij > up (default: 0.9)
+  -ermd RES_MIN_DIST    The minimum separation between two residues to have an RB the default value is set by the force-field. (default: None)
+  -ea RB_DECAY_FACTOR   Elastic bond decay factor a (default: 0)
+  -ep RB_DECAY_POWER    Elastic bond decay power p (default: 1)
+  -em RB_MINIMUM_FORCE  Remove elastic bonds with force constant lower than this (default: 0)
+  -eb RB_SELECTION      Comma separated list of bead names for elastic bonds (default: None)
+  -eunit RB_UNIT        Establish what is the structural unit for the elastic network. Bonds are only created within a unit. Options are molecule, chain,
+                        all, or aspecified region defined by resids, with followingformat: <start_resid_1>:<end_resid_1>, <start_resid_2>:<end_resid_2>...
+                        (default: molecule)
+
+Let's have a look at how to combine these in more detail.
+
+
+Basic usage
+-----------
+Without any further consideration, an elastic network can be added to your martinize2 command easily:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -elastic``
+
+which will apply the default harmonic bond constant (500 kJ/mol/nm^2) between non-successive BB beads
+which are < 0.9 nm apart in space.
+
+NOTE! For proteins in martini 3, the default constant should be 700 kJ/mol/nm^2. Changing the default
+elastic constant per force field will be fixed in future versions of martinize2.
+
+
+Customising cutoffs
+-------------------
+
+If your system requires more of a custom cutoff to better reproduce the dynamics of your protein,
+the region for the force to be applied in can be customised using the ``-el`` and ``-eu`` flags:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -elastic -el 0.1 -eu 0.5``
+
+In this example, the elastic network will only be applied between backbone beads which are between 0.1 and 0.5 nm
+apart.
+
+Using decays
+------------
+
+The strength of the elastic bond can be tuned with distance using an exponential decay function,
+which uses the ``-ea`` and ``-ep`` flags as input parameters:
+
+
+$$decay = e^{(- f * ((x - l) ^ p)}$$
+
+where:
+
+- ``l`` = lower bound  (``-el``)
+- ``f`` = decay factor (``-ea``)
+- ``p`` = decay power  (``-ep``)
+
+Combining parameters
+--------------------
+
+
+.. image:: elastic_examples.png
+
+Here we see several examples of how the strength of elastic network harmonic bonds can be tuned.
+
+The first example uses default parameters (albeit a shorter cutoff), such that a constant force is
+applied between all backbone beads within the cutoff.
+
+The second and third examples use a slightly longer cutoff, and apply a gentle decay function
+to the strength of the network. In the first case, the decay is applied naively, and as such its
+strength decays from 0 distance. In the second case, combining the decay with a lower cutoff means that
+for backbone beads that are close the elastic network strength is constand, but is lower between pairs slightly
+further away.
+
+The fourth example shows a similar function to the second example, but with a longer cutoff and a stronger decay.
+(note the form of the exponential decay above)
+
+The fifth example adds an additional parameter ``-em`` into the function. As described in the help, if forces are
+calculated to be lower than this force, they are removed and set to zero. Note how the input values are almost identical
+to the fourth example, which would otherwise get cutoff at 0.9 nm. Because the decay function reduces the force below
+the minimum before the cutoff, it overrides it and the force is zeroed before the upper cutoff anyway.
+
+
+Defining structural units
+-------------------------
+
+By default, martinize2 will look at the structure given in the input file, and construct a distance-based elastic
+network, filtered by each molecule. This behaviour is controlled by the `-eunit` flag. If you have multiple molecules
+within your input file and would like the way the elastic network is written to be changed, this can be achieved
+through different specifications as described in the help above.
+
+For example:
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -elastic -eunit chain``
+
+will limit the unit to individual chains in the system. *i.e.* chain A of your protein will *not* have any elastic
+bonds with chain B, and so on.
+
+Conversely,
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -elastic -eunit all``
+
+will write elastic bonds between every molecule in your system in the positions that have been found.
+
+Finally:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -elastic -eunit 1:100 150:200``
+
+Will write elastic networks internally between residues 1 to 100, and residues 150 to 200, but *not* between either of
+these domains, nor between either of these domains and residues 101 to 149.
+
+
+Visualising elastic networks
+----------------------------
+
+If you want to look at your elastic network in VMD to confirm that it's been constructed in the
+way that you're expecting, the `MartiniGlass <https://github.com/Martini-Force-Field-Initiative/MartiniGlass>`_
+package can help write visualisable topologies to view.

doc/source/tutorials/go_models.rst

@@ -0,0 +1,133 @@
+=========
+Go models
+=========
+
+GoMartini is a popular method for retaining the secondary and tertiary structures of proteins originating from the lab
+of `Adolfo Poma <https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00986>`_. In contrast to an elastic network, the Go
+model enforces interactions between specific pairs of beads within a protein based on residue overlap and restricted
+chemical structural unit criteria.
+
+The latest version of Martinize2 implements the newest version of the
+`Go model <https://www.biorxiv.org/content/10.1101/2024.04.15.589479v1>`_. In this version of the Go model, interactions
+are mediated through the addition of extra virtual sites on top of backbone beads in the protein. Interactions are in
+the form of Lennard-Jones interactions, which are written as an extra file to be included in the protein topology.
+
+The Go 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)
+   -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)
+
+To add a Go 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/>`_.
+
+The contact map is then used in your martinize2 command:
+
+``martinize2 -f protein.pdb -o topol.top -x cg_protein.pdb -ff martini3001 -dssp -go contact_map.out``
+
+
+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.
+ 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
+    molecular graph. These options are flexible through the ``-go-low`` and ``-go-up`` flags.
+ 4) If the contact map finds any atoms within contact range defined, but are *also* within 3 residues of each other,
+    then the contacts are removed. This is defined through the ``-go-res-dist`` flag.
+
+As a result, along with the standard output of martinize2 (*i.e.* itp files for your molecules, a generic .top file,
+a coarse grained structure file), you will get two extra files: ``go_atomtypes.itp`` and ``go_nbparams.itp``. The atomtypes
+file defines the new virtual sites as atoms for your system, and the nbparams file defines specific non-bonded
+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
+ ...
+
+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
+ ...
+
+To activate your Go model for use in Gromacs, the `martini_v3.0.0.itp` master itp needs the additional files included.
+The additional atomtypes defined in the ``go_atomtypes.itp`` file should be included at the end of the `[ atomtypes ]`
+directive as::
+
+
+ [ atomtypes ]
+ ...
+ TX1er 36.0 0.000 A 0.0 0.0
+ W  72.0 0.000 A 0.0 0.0
+ SW 54.0 0.000 A 0.0 0.0
+ TW 36.0 0.000 A 0.0 0.0
+ U  24.0 0.000 A 0.0 0.0
+
+ #ifdef GO_VIRT
+ #include "go_atomtypes.itp"
+ #endif
+
+ [ nonbond_params ]
+ P6    P6  1 4.700000e-01    4.990000e+00
+ P6    P5  1 4.700000e-01    4.730000e+00
+ P6    P4  1 4.700000e-01    4.480000e+00
+ ...
+
+Similarly, the nonbonded parameters should be included at the end of the `[ nonbond_params ]`
+directive::
+
+ ...
+ TX2er  SQ1n  1 3.660000e-01    3.528000e+00
+ TX2er  TQ1n  1 3.520000e-01    5.158000e+00
+ TX1er   Q1n  1 3.950000e-01    1.981000e+00
+ TX1er  SQ1n  1 3.780000e-01    3.098000e+00
+ TX1er  TQ1n  1 3.660000e-01    4.422000e+00
+
+ #ifdef GO_VIRT
+ #include "go_nbparams.itp"
+ #endif
+
+Then in the .top file for your system, simply include `#define GO_VIRT` along with the other files
+to be included to active the Go network in your model.
+
+As a shortcut for writing the include statements above, you can simply include these files in your master
+``martini_v3.0.0.itp`` file with the following commands::
+
+ sed -i "s/\[ nonbond_params \]/\#ifdef GO_VIRT\n\#include \"go_atomtypes.itp\"\n\#endif\n\n\[ nonbond_params \]/" martini_v3.0.0.itp
+
+ echo -e "\n#ifdef GO_VIRT \n#include \"go_nbparams.itp\"\n#endif" >> martini_v3.0.0.itp
+
+The Go model should then be usable in your simulations following the `general protein tutorial <https://pubs.acs.org/doi/10.1021/acs.jctc.4c00677>`_.
+But careful! While the Go model specifies nonbonded interactions, the interactions are only defined
+internally for each molecule. This means that if you have multiple copies of your Go model protein
+in the system, the Go bonds are still only specified internally for each copy of the molecule,
+not truly as intermolecular forces in the system as a whole. For more detail on this phenomenon,
+see the paper by `Korshunova et al. <https://pubs.acs.org/doi/10.1021/acs.jctc.4c00677>`_.
+
+
+Visualising Go networks
+----------------------------
+
+If you want to look at your Go network in VMD to confirm that it's been constructed in the
+way that you're expecting, the `MartiniGlass <https://github.com/Martini-Force-Field-Initiative/MartiniGlass>`_
+package can help write visualisable topologies to view.

doc/source/tutorials/index.rst

@@ -1,11 +1,19 @@
 Tutorials
 =========
-You can find some examples on how to use martinize 2 in the martinize-examples
+This page contains several examples of how martinize2 can be used to convert
+proteins from atomistic to a martini3 representation.
+
+You can find further examples on how to use martinize2 in the martinize-examples
 repository: https://github.com/marrink-lab/martinize-examples
 
 .. Organize the tutorials in directories, so it's easier to keep their files
     together
 
 .. toctree::
+    basic_usage
+    elastic_networks
+    go_models
+    water_biasing
+    mutations_and_modifications
     6_adding_residues_links/index
     7_adding_modifications/index