|
| 1 | +Additional Settings |
| 2 | +=================== |
| 3 | + |
| 4 | +Restraints |
| 5 | +########## |
| 6 | + |
| 7 | +.. danger:: |
| 8 | + This section only applies if you are running your simulations with **openMM**. |
| 9 | + Should you run your simulations using CHARMM, it will not apply the restraints |
| 10 | + and give **no warning about it**. |
| 11 | + |
| 12 | + |
| 13 | + |
| 14 | +|trafo| supports two types of restraints: automatic and manual restraints. |
| 15 | + |
| 16 | +Automatic Restraints |
| 17 | +********************* |
| 18 | + |
| 19 | +To activate automatic restraints for your ligand, add |
| 20 | + |
| 21 | +.. code-block:: yaml |
| 22 | + |
| 23 | + simulation: |
| 24 | + restraints: "auto" |
| 25 | + |
| 26 | +to your `config.yaml`. This wil restrain your ligand in its original position using a harmonic potential of the form :math:`E=0.5k \cdot (r-r_0)²` (where :math:`r`` is the current distance,:math:`r_0`` the initial distance, and :math:`k`` the force constant), applied within the centers of mass of your ligand and the surrounding protein structure, keeping these two vaguely fixed. |
| 27 | + |
| 28 | +You may specify further keywords: |
| 29 | + |
| 30 | +.. rubric:: Options for automatic restraints |
| 31 | + |
| 32 | +.. object:: restraints: auto |
| 33 | + |
| 34 | + |
| 35 | + .. option:: k=[int] |
| 36 | + |
| 37 | + *optional:* Defines the spring constant used for force calculations. Default is 3 |
| 38 | + |
| 39 | + .. option:: extremities=[int] |
| 40 | + |
| 41 | + *optional:* If used, |trafo| will not restraint the entire common core but rather look for [int] extremities. These are then restrained to the surrounding protein carbon-alphas. |
| 42 | + |
| 43 | + .. option:: shape=["harmonic","flatbottom"]` |
| 44 | + |
| 45 | + *optional:* Defines the shape of the energy potential. |
| 46 | + |
| 47 | + The potential energy term for the two forms is given as: |
| 48 | + |
| 49 | + + harmonic: :math:`E=0.5k \cdot (r-r_0)²` |
| 50 | + + flatbottom: :math:`E=k \cdot step(|r-r_0|),thresh) \cdot (r-r_0)²,` |
| 51 | + |
| 52 | + where :math:`r`` is the current distance, :math:`r_0`` the initial distance, and :math:`k`` the force constant. The step function returns 0 until the threshold :math:`thresh` (defined as :code:`wellsize` by the config.yaml) has been surpassed, after which it returns 1. As such, |
| 53 | + the flat-bottom potential is essentially a harmonic potential, where the energy within the \"well\" is zero. Note that the :code:`wellsize` acts like a *radius*, meaning it extends in all directions from its origin. |
| 54 | + |
| 55 | + .. figure:: assets/images/harm_vs_fb.png |
| 56 | + :alt: harmonic and flatbottom potentials |
| 57 | + |
| 58 | + Figure: The difference between a harmonic (*A*) and flat-bottom (*B*) potential. While both allow more movement closer to the restraint origin, a flat-bottom potential effectively allows defining |
| 59 | + a fixed sampling area that the ligand can move in without hindrance, but cannot leave (assuming a strong enough value for *k*). |
| 60 | + |
| 61 | + .. option:: scaling |
| 62 | + |
| 63 | + *optional:* If present, the k - Value of the force is linearly scaled in the first four intermediate states (``intst1: 0; intst2: 0.25; intst3: 0.5; intst4: 0.75; intst5: 1.0``) |
| 64 | + |
| 65 | + .. option:: wellsize [float] |
| 66 | + |
| 67 | + *optional*: Only takes effect in a flat-bottom potential. Defines the wellsize in nanometers. Default is 0.1 nm. |
| 68 | + |
| 69 | +A full entry in your config.yaml might thus look like this: |
| 70 | + |
| 71 | + |
| 72 | + |
| 73 | +.. code-block:: yaml |
| 74 | + |
| 75 | + restraints: "auto k=10 extremities=3 shape=harmonic scaling" |
| 76 | + |
| 77 | + |
| 78 | + |
| 79 | +.. caution:: Be somewhat sure of what your structure looks like, and do a sanity check on the generated restraints before production. As all restraints only act on the common core, setting an arbitrarily high number of extermities can lead to strange results |
| 80 | + |
| 81 | +It should be noted that this means that a small file called `restraints.yaml` is created in your `intst*` - folders. |
| 82 | +These have the following structure: |
| 83 | + |
| 84 | + |
| 85 | +.. code-block:: yaml |
| 86 | + |
| 87 | + system: |
| 88 | + structure: |
| 89 | + tlc: LIG # same as in the config.yaml, but only one structure (as only one relevant) |
| 90 | + |
| 91 | + simulation: |
| 92 | + restraints: "auto" # same as in config.yaml |
| 93 | + ccs: # this represents an array of your common core, upon which restraints can be applied |
| 94 | + - C1 |
| 95 | + - C2 |
| 96 | + - H2 |
| 97 | + intst: |
| 98 | + scaling:0.8 # for non-immediate switches, how far along the scaling is. Only relevant for harmonic potentials. |
| 99 | + |
| 100 | + |
| 101 | +It is not recommended to manually edit these files, as they are automatically created for each intermediate state. |
| 102 | + |
| 103 | +Manual Restraints |
| 104 | +***************** |
| 105 | + |
| 106 | +To activate manual restraints for your ligand, add |
| 107 | + |
| 108 | +.. code-block:: yaml |
| 109 | + |
| 110 | + simulation: |
| 111 | + restraints: "manual" |
| 112 | + |
| 113 | +to your config.yaml. Below, you may now specify an arbitrary number |
| 114 | +of restraints using the |
| 115 | +`MDAnalysis selection syntax <https://docs.mdanalysis.org/stable/documentation_pages/selections.html#simple-selections>`_ : |
| 116 | + |
| 117 | +.. code-block:: yaml |
| 118 | + |
| 119 | + simulation: |
| 120 | + restraints: "manual" |
| 121 | + manualrestraints: |
| 122 | + restraint1: |
| 123 | + shape: "harmonic" |
| 124 | + group1: "resname LIG and type C" |
| 125 | + group2: "protein and type CA" |
| 126 | + k: 30 |
| 127 | + r0: 2.41 |
| 128 | + |
| 129 | +You may define as many restraints as you like: |
| 130 | + |
| 131 | +Code example with multiple restraints: |
| 132 | + |
| 133 | +.. code-block:: yaml |
| 134 | + |
| 135 | + simulation: |
| 136 | + restraints: "manual" |
| 137 | + manualrestraints: |
| 138 | + restraint1: |
| 139 | + shape: "harmonic" |
| 140 | + group1: "resname LIG and type C" |
| 141 | + group2: "protein and type CA" |
| 142 | + restraint2: |
| 143 | + shape: "flatbottom" |
| 144 | + group1: "resname LIG and type C" |
| 145 | + group2: "protein and type CA" |
| 146 | + restraint3: |
| 147 | + shape: "harmonic" |
| 148 | + group1: "resname LIG and name C14" |
| 149 | + group2: "sphlayer 5 15 name C14 and protein and type CA" |
| 150 | + |
| 151 | +Note that the individual restraints all need to have distinct names (restraint1, restraint2 etc.). It is common that they are numbered, but not required - they simply need to adhere to the yaml syntax. |
| 152 | + |
| 153 | +.. rubric:: Options for manual restraints |
| 154 | + |
| 155 | + |
| 156 | +.. object:: restraints: manual |
| 157 | + |
| 158 | + .. object:: manual-restraint |
| 159 | + |
| 160 | + You may freely choose the name of the restraint here. It may be useful to provide a 'speaking' name, as this will allow identification later. |
| 161 | + |
| 162 | + .. option:: shape=["harmonic","flatbottom"]' |
| 163 | + |
| 164 | + Shape of the energy potential. Default is "harmonic". See automatic restraints for details. |
| 165 | + |
| 166 | + .. option:: group1,group2=[MDAnalysis selection string] |
| 167 | + |
| 168 | + Defines which Common Core atoms are members of group1 or group2. Please note that group1 **must** be the ligand, and group2 the protein. |
| 169 | + |
| 170 | + .. option:: k=[int] |
| 171 | + |
| 172 | + *(optional):* Defines the harmonic force constant. Default is 3. |
| 173 | + |
| 174 | + .. option:: wellsize=[float] |
| 175 | + |
| 176 | + *(optional):* Defines the wellsize for flat-bottom restraints (unit is nanometers). Defaults to 0.1 nm. |
| 177 | + |
| 178 | + |
| 179 | +As with automatic restraints, even manually specified restraints will never act on atoms not in the common core, as this would lead to nonsensical energy calculations. |
| 180 | + |
| 181 | + |
| 182 | +Hydrogen Mass Repartitioning (HMR) |
| 183 | +#################################### |
| 184 | + |
| 185 | +To further accelerate MD simulations it is possible to reweight the masses of the hydrogen atoms and thus |
| 186 | +reduce the vibrational frequency of the corresponding hydrogen-heavy atom bond. When using ``cons: HBonds`` |
| 187 | +one can safely increase the time step to 4 fs. |
| 188 | +To use HMR one can either check the box in the last step of the CHARMM-GUI solution builder while creating the |
| 189 | +system or one can use the `parmed <https://parmed.github.io/ParmEd/html/api/parmed/parmed.tools.actions.html#parmed.tools.actions.HMassRepartition>`_ |
| 190 | +action tool (also available in the ``fep`` conda environment) |
| 191 | + |
| 192 | +.. code-block:: python |
| 193 | + |
| 194 | + psf = pm.charmm.CharmmPsfFile("step3_input_orig.psf") |
| 195 | + pm.tools.actions.HMassRepartition(psf).execute() |
| 196 | + psf.save("step3_input.psf", overwrite = True) |
| 197 | + |
| 198 | +.. |trafo| replace:: :math:`\texttt{TRANSFORMATO}` |
| 199 | + |
0 commit comments