Reimplement constant acceleration groups
[alexxy/gromacs.git] / docs / user-guide / mdp-options.rst
index 36c2ee104ea6a717433862e7c0377d82d45e3f94..a1b96fa84ac135e94694637e732a0e86f0c1580d 100644 (file)
@@ -231,6 +231,43 @@ Run control
          same simulation. This option is generally useful to set only
          when coping with a crashed simulation where files were lost.
 
+.. mdp:: mts
+
+   .. mdp-value:: no
+
+      Evaluate all forces at every integration step.
+
+   .. mdp-value:: yes
+
+      Use a multiple timing-stepping integrator to evaluate some forces, as specified
+      by :mdp:`mts-level2-forces` every :mdp:`mts-level2-factor` integration
+      steps. All other forces are evaluated at every step. MTS is currently
+      only supported with :mdp-value:`integrator=md`.
+
+.. mdp:: mts-levels
+
+        (2)
+       The number of levels for the multiple time-stepping scheme.
+       Currently only 2 is supported.
+
+.. mdp:: mts-level2-forces
+
+   (longrange-nonbonded)
+   A list of one or more force groups that will be evaluated only every
+   :mdp:`mts-level2-factor` steps. Supported entries are:
+   ``longrange-nonbonded``, ``nonbonded``, ``pair``, ``dihedral``, ``angle``,
+   ``pull`` and ``awh``. With ``pair`` the listed pair forces (such as 1-4)
+   are selected. With ``dihedral`` all dihedrals are selected, including cmap.
+   All other forces, including all restraints, are evaluated and
+   integrated every step. When PME or Ewald is used for electrostatics
+   and/or LJ interactions, ``longrange-nonbonded`` can not be omitted here.
+
+.. mdp:: mts-level2-factor
+
+      (2) [steps]
+      Interval for computing the forces in level 2 of the multiple time-stepping
+      scheme
+
 .. mdp:: comm-mode
 
    .. mdp-value:: Linear
@@ -359,24 +396,30 @@ Output control
    (0) [steps]
    number of steps that elapse between writing coordinates to the output
    trajectory file (:ref:`trr`), the last coordinates are always written
+   unless 0, which means coordinates are not written into the trajectory
+   file.
 
 .. mdp:: nstvout
 
    (0) [steps]
    number of steps that elapse between writing velocities to the output
    trajectory file (:ref:`trr`), the last velocities are always written
+   unless 0, which means velocities are not written into the trajectory
+   file.
 
 .. mdp:: nstfout
 
    (0) [steps]
    number of steps that elapse between writing forces to the output
-   trajectory file (:ref:`trr`), the last forces are always written.
+   trajectory file (:ref:`trr`), the last forces are always written,
+   unless 0, which means forces are not written into the trajectory
+   file.
 
 .. mdp:: nstlog
 
    (1000) [steps]
    number of steps that elapse between writing energies to the log
-   file, the last energies are always written
+   file, the last energies are always written.
 
 .. mdp:: nstcalcenergy
 
@@ -401,7 +444,8 @@ Output control
 
    (0) [steps]
    number of steps that elapse between writing position coordinates
-   using lossy compression (:ref:`xtc` file)
+   using lossy compression (:ref:`xtc` file), 0 for not writing
+   compressed coordinates output.
 
 .. mdp:: compressed-x-precision
 
@@ -448,7 +492,9 @@ Neighbor searching
       a minimum value and :ref:`gmx mdrun` might increase it, unless
       it is set to 1. With parallel simulations and/or non-bonded
       force calculation on the GPU, a value of 20 or 40 often gives
-      the best performance.
+      the best performance. With energy minimization this parameter
+      is not used as the pair list is updated when at least one atom
+      has moved by more than half the pair list buffer size.
 
    .. mdp-value:: 0
 
@@ -803,10 +849,10 @@ Tables
 
    Currently unsupported.
    When user tables are used for electrostatics and/or VdW, here one
-   can give pairs of energy groups for which seperate user tables
+   can give pairs of energy groups for which separate user tables
    should be used. The two energy groups will be appended to the table
    file name, in order of their definition in :mdp:`energygrps`,
-   seperated by underscores. For example, if ``energygrps = Na Cl
+   separated by underscores. For example, if ``energygrps = Na Cl
    Sol`` and ``energygrp-table = Na Na Na Cl``, :ref:`gmx mdrun` will
    read ``table_Na_Na.xvg`` and ``table_Na_Cl.xvg`` in addition to the
    normal ``table.xvg`` which will be used for all other energy group
@@ -892,9 +938,9 @@ Ewald
    .. mdp-value:: 3dc
 
       The reciprocal sum is still performed in 3D, but a force and
-      potential correction applied in the `z` dimension to produce a
+      potential correction applied in the ``z`` dimension to produce a
       pseudo-2D summation. If your system has a slab geometry in the
-      `x-y` plane you can try to increase the `z`-dimension of the box
+      ``x-y`` plane you can try to increase the ``z``-dimension of the box
       (a box height of 3 times the slab height is usually ok) and use
       this option.
 
@@ -969,8 +1015,10 @@ Temperature coupling
 
    (-1)
    The frequency for coupling the temperature. The default value of -1
-   sets :mdp:`nsttcouple` equal to :mdp:`nstlist`, unless
-   :mdp:`nstlist` <=0, then a value of 10 is used. For velocity
+   sets :mdp:`nsttcouple` equal to 10, or fewer steps if required
+   for accurate integration. Note that the default value is not 1
+   because additional computation and communication is required for
+   obtaining the kinetic energy. For velocity
    Verlet integrators :mdp:`nsttcouple` is set to 1.
 
 .. mdp:: nh-chain-length
@@ -1027,6 +1075,13 @@ Pressure coupling
       ensemble, but it is the most efficient way to scale a box at the
       beginning of a run.
 
+   .. mdp-value:: C-rescale
+
+      Exponential relaxation pressure coupling with time constant
+      :mdp:`tau-p`, including a stochastic term to enforce correct
+      volume fluctuations.  The box is scaled every :mdp:`nstpcouple`
+      steps. It can be used for both equilibration and production.
+
    .. mdp-value:: Parrinello-Rahman
 
       Extended-ensemble pressure coupling where the box vectors are
@@ -1088,14 +1143,14 @@ Pressure coupling
    .. mdp-value:: surface-tension
 
       Surface tension coupling for surfaces parallel to the
-      xy-plane. Uses normal pressure coupling for the `z`-direction,
-      while the surface tension is coupled to the `x/y` dimensions of
+      xy-plane. Uses normal pressure coupling for the ``z``-direction,
+      while the surface tension is coupled to the ``x/y`` dimensions of
       the box. The first :mdp:`ref-p` value is the reference surface
       tension times the number of surfaces ``bar nm``, the second
-      value is the reference `z`-pressure ``bar``. The two
+      value is the reference ``z``-pressure ``bar``. The two
       :mdp:`compressibility` values are the compressibility in the
-      `x/y` and `z` direction respectively. The value for the
-      `z`-compressibility should be reasonably accurate since it
+      ``x/y`` and ``z`` direction respectively. The value for the
+      ``z``-compressibility should be reasonably accurate since it
       influences the convergence of the surface-tension, it can also
       be set to zero to have a box with constant height.
 
@@ -1103,8 +1158,10 @@ Pressure coupling
 
    (-1)
    The frequency for coupling the pressure. The default value of -1
-   sets :mdp:`nstpcouple` equal to :mdp:`nstlist`, unless
-   :mdp:`nstlist` <=0, then a value of 10 is used. For velocity
+   sets :mdp:`nstpcouple` equal to 10, or fewer steps if required
+   for accurate integration. Note that the default value is not 1
+   because additional computation and communication is required for
+   obtaining the virial. For velocity
    Verlet integrators :mdp:`nstpcouple` is set to 1.
 
 .. mdp:: tau-p
@@ -1462,6 +1519,7 @@ Walls
 COM pulling
 ^^^^^^^^^^^
 
+Sets whether pulling on collective variables is active.
 Note that where pulling coordinates are applicable, there can be more
 than one (set with :mdp:`pull-ncoords`) and multiple related :ref:`mdp`
 variables will exist accordingly. Documentation references to things
@@ -1631,7 +1689,8 @@ pull-coord2-vec, pull-coord2-k, and so on.
       Center of mass pulling using a constraint between the reference
       group and one or more groups. The setup is identical to the
       option umbrella, except for the fact that a rigid constraint is
-      applied instead of a harmonic potential.
+      applied instead of a harmonic potential. Note that this type is
+      not supported in combination with multiple time stepping.
 
    .. mdp-value:: constant-force
 
@@ -1734,6 +1793,29 @@ pull-coord2-vec, pull-coord2-k, and so on.
       then defined as the angle between two planes: the plane spanned by the
       the two first vectors and the plane spanned the two last vectors.
 
+   .. mdp-value:: transformation
+
+      Transforms other pull coordinates using a mathematical expression defined by :mdp:`pull-coord1-expression`.
+      Pull coordinates of lower indices can be used as variables to this pull coordinate.
+      Thus, pull transformation coordinates should have a higher pull coordinate index
+      than all pull coordinates they transform.
+
+.. mdp:: pull-coord1-expression
+
+   Mathematical expression to transform pull coordinates of lower indices to a new one.
+   The pull coordinates are referred to as variables in the equation so that
+   pull-coord1's value becomes 'x1', pull-coord2 value becomes 'x2' etc.
+   The mathematical expression are evaluated using muParser.
+   Only relevant if :mdp:`pull-coord1-geometry` is set to :mdp-value:`transformation`.
+
+.. mdp:: pull-coord1-dx
+
+   (1e-9)
+   Size of finite difference to use in numerical derivation of the pull coordinate
+   with respect to other pull coordinates.
+   The current implementation uses a simple first order finite difference method to perform derivation so that
+   f'(x) = (f(x+dx)-f(x))/dx
+   Only relevant if :mdp:`pull-coord1-geometry` is set to :mdp-value:`transformation`.
 
 .. mdp:: pull-coord1-groups
 
@@ -1837,7 +1919,8 @@ AWH adaptive biasing
       The applied biasing potential is the convolution of the bias function and a
       set of harmonic umbrella potentials (see :mdp-value:`awh-potential=umbrella` below). This results
       in a smooth potential function and force. The resolution of the potential is set
-      by the force constant of each umbrella, see :mdp:`awh1-dim1-force-constant`.
+      by the force constant of each umbrella, see :mdp:`awh1-dim1-force-constant`. This option is not
+      compatible with using the free energy lambda state as an AWH reaction coordinate.
 
    .. mdp-value:: umbrella
 
@@ -1845,8 +1928,9 @@ AWH adaptive biasing
       using Monte-Carlo sampling.  The force constant is set with
       :mdp:`awh1-dim1-force-constant`. The umbrella location
       is sampled using Monte-Carlo every :mdp:`awh-nstsample` steps.
-      There are no advantages to using an umbrella.
-      This option is mainly for comparison and testing purposes.
+      This is option is required when using the free energy lambda state as an AWH reaction coordinate.
+      Apart from that, this option is mainly for comparison
+      and testing purposes as there are no advantages to using an umbrella.
 
 .. mdp:: awh-share-multisim
 
@@ -1999,7 +2083,7 @@ AWH adaptive biasing
       The file name can be changed with the ``-awh`` option.
       The first :mdp:`awh1-ndim` columns of
       each input file should contain the coordinate values, such that each row defines a point in
-      coordinate space. Column :mdp:`awh1-ndim` + 1 should contain the PMF value for each point.
+      coordinate space. Column :mdp:`awh1-ndim` + 1 should contain the PMF value (in kT) for each point.
       The target distribution column can either follow the PMF (column  :mdp:`awh1-ndim` + 2) or
       be in the same column as written by :ref:`gmx awh`.
 
@@ -2011,15 +2095,15 @@ AWH adaptive biasing
 
    .. mdp-value:: positive
 
-      Share the bias and PMF estimates within and/or between simulations.
-      Within a simulation, the bias will be shared between biases that have the
-      same :mdp:`awh1-share-group` index (note that the current code does not support this).
-      With :mdp-value:`awh-share-multisim=yes` and
-      :ref:`gmx mdrun` option ``-multidir`` the bias will also be shared across simulations.
+      Share the bias and PMF estimates between simulations. This currently
+      only works between biases with the same index. Note that currently
+      sharing within a single simulation is not supported.
+      The bias will be shared across simulations that specify the same
+      value for :mdp:`awh1-share-group`. To enable this, use
+      :mdp-value:`awh-share-multisim=yes` and the :ref:`gmx mdrun` option
+      ``-multidir``.
       Sharing may increase convergence initially, although the starting configurations
       can be critical, especially when sharing between many biases.
-      Currently, positive group values should start at 1 and increase
-      by 1 for each subsequent bias that is shared.
 
 .. mdp:: awh1-ndim
 
@@ -2033,8 +2117,16 @@ AWH adaptive biasing
 
    .. mdp-value:: pull
 
-      The module providing the reaction coordinate for this dimension.
-      Currently AWH can only act on pull coordinates.
+      The pull module is providing the reaction coordinate for this dimension.
+      With multiple time-stepping, AWH and pull should be in the same MTS level.
+
+   .. mdp-value:: fep-lambda
+
+      The free energy lambda state is the reaction coordinate for this dimension.
+      The lambda states to use are specified by :mdp:`fep-lambdas`, :mdp:`vdw-lambdas`,
+      :mdp:`coul-lambdas` etc. This is not compatible with delta-lambda. It also requires
+      calc-lambda-neighbors to be -1. With multiple time-stepping, AWH should
+      be in the slow level. This option requires :mdp-value:`awh-potential=umbrella`.
 
 .. mdp:: awh1-dim1-coord-index
 
@@ -2066,7 +2158,7 @@ AWH adaptive biasing
 
 .. mdp:: awh1-dim1-diffusion
 
-   (10\ :sup:`-5`) [nm\ :sup:`2`/ps] or [rad\ :sup:`2`/ps]
+   (10\ :sup:`-5`) [nm\ :sup:`2`/ps], [rad\ :sup:`2`/ps] or [ps\ :sup:`-1`]
    Estimated diffusion constant for this coordinate dimension determining the initial
    biasing rate. This needs only be a rough estimate and should not critically
    affect the results unless it is set to something very low, leading to slow convergence,
@@ -2247,7 +2339,7 @@ NMR refinement
 
    (1000) [kJ mol\ :sup:`-1` nm\ :sup:`-2`]
    force constant for distance restraints, which is multiplied by a
-   (possibly) different factor for each restraint given in the `fac`
+   (possibly) different factor for each restraint given in the ``fac``
    column of the interaction in the topology file.
 
 .. mdp:: disre-tau
@@ -2439,20 +2531,31 @@ Free energy calculations
    written out. For normal BAR such as with :ref:`gmx bar`, a value of
    1 is sufficient, while for MBAR -1 should be used.
 
+.. mdp:: sc-function
+
+   (beutler)
+
+   .. mdp-value:: beutler
+
+   Beutler *et al.* soft-core function
+
+   .. mdp-value:: gapsys
+
+   Gapsys *et al.* soft-core function
+
 .. mdp:: sc-alpha
 
    (0)
-   the soft-core alpha parameter, a value of 0 results in linear
-   interpolation of the LJ and Coulomb interactions
+   for `sc-function=beutler` the soft-core alpha parameter,
+   a value of 0 results in linear interpolation of the
+   LJ and Coulomb interactions.
+   Used only with `sc-function=beutler`
 
 .. mdp:: sc-r-power
 
    (6)
    power 6 for the radial term in the soft-core equation.
-
-   (48)
-   (deprecated) power 48 for the radial term in the soft-core equation. 
-   Note that sc-alpha should generally be much lower (between 0.001 and 0.003).
+   Used only with `sc-function=beutler`
 
 .. mdp:: sc-coul
 
@@ -2465,18 +2568,47 @@ Free energy calculations
    states are used, not with :mdp:`couple-lambda0` /
    :mdp:`couple-lambda1`, and you can still turn off soft-core
    interactions by setting :mdp:`sc-alpha` to 0.
+   Used only with `sc-function=beutler`
 
 .. mdp:: sc-power
 
    (0)
    the power for lambda in the soft-core function, only the values 1
-   and 2 are supported
+   and 2 are supported. Used only with `sc-function=beutler`
 
 .. mdp:: sc-sigma
 
    (0.3) [nm]
-   the soft-core sigma for particles which have a C6 or C12 parameter
-   equal to zero or a sigma smaller than :mdp:`sc-sigma`
+   for `sc-function=beutler` the soft-core sigma for particles
+   which have a C6 or C12 parameter equal to zero or a sigma smaller
+   than :mdp:`sc-sigma`.
+   Used only with `sc-function=beutler`
+
+.. mdp:: sc-gapsys-scale-linpoint-lj
+
+   (0.85)
+   for `sc-function=gapsys` it is the unitless alphaLJ parameter.
+   It controls the softness of the van der Waals interactions
+   by scaling the point for linearizing the vdw force.
+   Setting it to 0 will result in the standard hard-core
+   van der Waals interactions.
+   Used only with `sc-function=gapsys`
+
+.. mdp:: sc-gapsys-scale-linpoint-q
+
+   (0.3) [nm/e^2]
+   For `sc-function=gapsys` the alphaQ parameter
+   with the unit of [nm/e^2] and default value of 0.3. It controls
+   the softness of the Coulombic interactions. Setting it to 0 will
+   result in the standard hard-core Coulombic interactions.
+   Used only with `sc-function=gapsys`
+
+.. mdp:: sc-gapsys-sigma-lj
+
+   (0.3) [nm]
+   for `sc-function=gapsys` the soft-core sigma for particles
+   which have a C6 or C12 parameter equal to zero.
+   Used only with `sc-function=gapsys`
 
 .. mdp:: couple-moltype
 
@@ -2895,7 +3027,10 @@ Non-equilibrium MD
 
    groups for constant acceleration (*e.g.* ``Protein Sol``) all atoms
    in groups Protein and Sol will experience constant acceleration as
-   specified in the :mdp:`accelerate` line
+   specified in the :mdp:`accelerate` line. Note that the kinetic energy
+   of the center of mass of accelarated groups contributes to the kinetic
+   energy and temperature of the system. If this is not desired, make
+   each accelerate group also a separate temperature coupling group.
 
 .. mdp:: accelerate
 
@@ -2918,8 +3053,8 @@ Non-equilibrium MD
 .. mdp:: freezedim
 
    dimensions for which groups in :mdp:`freezegrps` should be frozen,
-   specify `Y` or `N` for X, Y and Z and for each group (*e.g.* ``Y Y
-   N N N N`` means that particles in the first group can move only in
+   specify ``Y`` or ``N`` for X, Y and Z and for each group (*e.g.*
+   ``Y Y N N N N`` means that particles in the first group can move only in
    Z direction. The particles in the second group can move in any
    direction).
 
@@ -2971,8 +3106,8 @@ Electric fields
    with units (respectively) V nm\ :sup:`-1`, ps\ :sup:`-1`, ps, ps.
 
    In the special case that ``sigma = 0``, the exponential term is omitted
-   and only the cosine term is used. If also ``omega = 0`` a static
-   electric field is applied.
+   and only the cosine term is used. In this case, ``t0`` must be set to 0.
+   If also ``omega = 0`` a static electric field is applied.
 
    Read more at :ref:`electric fields` and in ref. \ :ref:`146 <refCaleman2008a>`.
 
@@ -2980,101 +3115,15 @@ Electric fields
 Mixed quantum/classical molecular dynamics
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-.. MDP:: QMMM
-
-   .. mdp-value:: no
-
-      No QM/MM.
-
-   .. mdp-value:: yes
-
-      Do a QM/MM simulation. Several groups can be described at
-      different QM levels separately. These are specified in the
-      :mdp:`QMMM-grps` field separated by spaces. The level of *ab
-      initio* theory at which the groups are described is specified by
-      :mdp:`QMmethod` and :mdp:`QMbasis` Fields. Describing the
-      groups at different levels of theory is only possible with the
-      ONIOM QM/MM scheme, specified by :mdp:`QMMMscheme`.
-
 .. mdp:: QMMM-grps
 
-   groups to be descibed at the QM level (works also in case of MiMiC QM/MM)
-
-.. mdp:: QMMMscheme
-
-   .. mdp-value:: normal
-
-      normal QM/MM. There can only be one :mdp:`QMMM-grps` that is
-      modelled at the :mdp:`QMmethod` and :mdp:`QMbasis` level of
-      *ab initio* theory. The rest of the system is described at the
-      MM level. The QM and MM subsystems interact as follows: MM point
-      charges are included in the QM one-electron hamiltonian and all
-      Lennard-Jones interactions are described at the MM level.
+   groups to be descibed at the QM level for MiMiC QM/MM
 
-   .. mdp-value:: ONIOM
-
-      The interaction between the subsystem is described using the
-      ONIOM method by Morokuma and co-workers. There can be more than
-      one :mdp:`QMMM-grps` each modeled at a different level of QM
-      theory (:mdp:`QMmethod` and :mdp:`QMbasis`).
-
-.. mdp:: QMmethod
-
-   (RHF)
-   Method used to compute the energy and gradients on the QM
-   atoms. Available methods are AM1, PM3, RHF, UHF, DFT, B3LYP, MP2,
-   CASSCF, and MMVB. For CASSCF, the number of electrons and orbitals
-   included in the active space is specified by :mdp:`CASelectrons`
-   and :mdp:`CASorbitals`.
-
-.. mdp:: QMbasis
-
-   (STO-3G)
-   Basis set used to expand the electronic wavefuntion. Only Gaussian
-   basis sets are currently available, *i.e.* ``STO-3G, 3-21G, 3-21G*,
-   3-21+G*, 6-21G, 6-31G, 6-31G*, 6-31+G*,`` and ``6-311G``.
-
-.. mdp:: QMcharge
-
-   (0) [integer]
-   The total charge in `e` of the :mdp:`QMMM-grps`. In case there are
-   more than one :mdp:`QMMM-grps`, the total charge of each ONIOM
-   layer needs to be specified separately.
-
-.. mdp:: QMmult
-
-   (1) [integer]
-   The multiplicity of the :mdp:`QMMM-grps`. In case there are more
-   than one :mdp:`QMMM-grps`, the multiplicity of each ONIOM layer
-   needs to be specified separately.
-
-.. mdp:: CASorbitals
-
-   (0) [integer]
-   The number of orbitals to be included in the active space when
-   doing a CASSCF computation.
-
-.. mdp:: CASelectrons
-
-   (0) [integer]
-   The number of electrons to be included in the active space when
-   doing a CASSCF computation.
-
-.. MDP:: SH
+.. MDP:: QMMM
 
    .. mdp-value:: no
 
-      No surface hopping. The system is always in the electronic
-      ground-state.
-
-   .. mdp-value:: yes
-
-      Do a QM/MM MD simulation on the excited state-potential energy
-      surface and enforce a *diabatic* hop to the ground-state when
-      the system hits the conical intersection hyperline in the course
-      the simulation. This option only works in combination with the
-      CASSCF method.
-
+      QM/MM is no longer supported via these .mdp options. For MiMic, use no here.
 
 Computational Electrophysiology
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -3234,6 +3283,11 @@ electron-microscopy experiments. (See the `reference manual`_ for details)
       between reference density and simulated density as similarity measure.
       Negative density values are ignored.
 
+   .. mdp-value:: cross-correlation
+
+      Uses the Pearson correlation coefficient between reference density and
+      simulated density as similarity measure.
+
 .. mdp:: density-guided-simulation-atom-spreading-weight
 
    (unity) Determines the multiplication factor for the Gaussian kernel when
@@ -3300,6 +3354,71 @@ electron-microscopy experiments. (See the `reference manual`_ for details)
    (4) [ps] Couple force constant to increase in similarity with reference density
    with this time constant. Larger times result in looser coupling.
 
+.. mdp:: density-guided-simulation-shift-vector
+
+   (0,0,0) [nm] Add this vector to all atoms in the 
+   density-guided-simulation-group before calculating forces and energies for
+   density-guided-simulations. Affects only the density-guided-simulation forces
+   and energies. Corresponds to a shift of the input density in the opposite
+   direction by (-1) * density-guided-simulation-shift-vector.
+
+.. mdp:: density-guided-simulation-transformation-matrix
+
+   (1,0,0,0,1,0,0,0,1) Multiply all atoms with this matrix in the 
+   density-guided-simulation-group before calculating forces and energies for
+   density-guided-simulations. Affects only the density-guided-simulation forces
+   and energies. Corresponds to a transformation of the input density by the
+   inverse of this matrix. The matrix is given in row-major order.
+   This option allows, e.g., rotation of the density-guided atom group around the
+   z-axis by :math:`\theta` degress by using following input:
+   :math:`(\cos \theta , -\sin \theta , 0 , \sin \theta , \cos \theta , 0 , 0 , 0 , 1)` .
+
+QM/MM simulations with CP2K Interface 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+These options enable and control the calculation and application of additional
+QM/MM forces that are computed by the CP2K package if it is linked into |Gromacs|.
+For further details about QM/MM interface implementation follow :ref:`qmmm`. 
+
+.. mdp:: qmmm-cp2k-active
+
+   (false) Activate QM/MM simulations. Requires CP2K to be linked with |Gromacs|
+
+.. mdp:: qmmm-cp2k-qmgroup
+
+   (System) Index group with atoms that are treated with QM.
+
+.. mdp:: qmmm-cp2k-qmmethod
+
+   (PBE) Method used to describe the QM part of the system.
+
+   .. mdp-value:: PBE
+
+      DFT using PBE functional and DZVP-MOLOPT basis set.
+
+   .. mdp-value:: BLYP
+
+      DFT using BLYP functional and DZVP-MOLOPT basis set.
+
+   .. mdp-value:: INPUT
+
+      Provide an external input file for CP2K when running :ref:`gmx grompp` with the ``-qmi`` command-line option.
+      External input files are subject to the limitations that are described in :ref:`qmmm`.
+
+.. mdp:: qmmm-cp2k-qmcharge
+
+   (0) Total charge of the QM part.
+
+.. mdp:: qmmm-cp2k-qmmultiplicity
+
+   (1) Multiplicity or spin-state of QM part. Default value 1 means singlet state.
+
+.. mdp:: qmmm-cp2k-qmfilenames
+
+   () Names of the CP2K files that will be generated during the simulation. 
+   When using the default, empty, value the name of the simulation input file will be used 
+   with an additional ``_cp2k`` suffix.
+
 User defined thingies
 ^^^^^^^^^^^^^^^^^^^^^