From 73c26ac6e248876de4796c6885bd653fb2c21e33 Mon Sep 17 00:00:00 2001 From: Dmitry Morozov Date: Tue, 5 Oct 2021 08:23:00 +0000 Subject: [PATCH] QMMM description and release notes update --- docs/reference-manual/references.rst | 23 +++ docs/reference-manual/special/qmmm.rst | 185 ++++++++++++++------- docs/release-notes/2022/major/features.rst | 12 ++ 3 files changed, 164 insertions(+), 56 deletions(-) diff --git a/docs/reference-manual/references.rst b/docs/reference-manual/references.rst index 0d2aa5f93a..2f32912f79 100644 --- a/docs/reference-manual/references.rst +++ b/docs/reference-manual/references.rst @@ -2636,6 +2636,29 @@ structures into cryoelectron microscopy maps using biased molecular dynamics sim :sup:`187` Lundborg M., Lidmar J. and Hess B., "The accelerated weight histogram method for alchemical free energy calculations", *J. Chem. Phys.*, **154**, 204103 (2021). +.. raw:: html + + + +
+ +.. _refcp2k2020: + +:sup:`188` Kühne T., Iannuzzi M., Del Ben M. and Hutter J. *et al.*, +"CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations", +*J. Chem. Phys.*, **152**, 194103 (2020). + +.. raw:: html + +
+ +
+ +.. _refLaino2005: + +:sup:`189` Laino T., Mohamed F., Laio A. and Parrinello M., +"An Efficient Real Space Multigrid QM/MM Electrostatic Coupling", *J. Chem. Theory Comput.*, **1**, 1176 (2005). + .. raw:: html
diff --git a/docs/reference-manual/special/qmmm.rst b/docs/reference-manual/special/qmmm.rst index a7f8a2a28f..a0bf71f208 100644 --- a/docs/reference-manual/special/qmmm.rst +++ b/docs/reference-manual/special/qmmm.rst @@ -1,7 +1,7 @@ .. _qmmm: -Mixed Quantum-Classical simulation techniques ---------------------------------------------- +Hybrid Quantum-Classical simulations (QM/MM) with CP2K interface +---------------------------------------------------------------- In a molecular mechanics (MM) force field, the influence of electrons is expressed by empirical parameters that are assigned on the basis of @@ -20,71 +20,144 @@ One approach to the simulation of chemical reactions in solution, or in enzymes, is to use a combination of quantum mechanics (QM) and molecular mechanics (MM). The reacting parts of the system are treated quantum mechanically, with the remainder being modeled using the force field. -The current version of |Gromacs| provides interfaces to several popular -Quantum Chemistry packages (MOPAC :ref:`150 `, -GAMESS-UK \ :ref:`151 `, Gaussian \ :ref:`152 ` and -CPMD \ :ref:`153 `). - -|Gromacs| interactions between the two subsystems are either handled as -described by Field et al. :ref:`154 ` or within -the ONIOM approach by Morokuma and coworkers \ :ref:`155 `, -:ref:`156 `. +The current version of |Gromacs| provides an interface to the popular +Quantum Chemistry package CP2K :ref:`188 `. Overview ^^^^^^^^ -Two approaches for describing the interactions between the QM and MM -subsystems are supported in this version: - -#. **Electronic Embedding** The electrostatic interactions between the - electrons of the QM region and the MM atoms and between the QM nuclei - and the MM atoms are included in the Hamiltonian for the QM - subsystem: +|Gromacs| interactions between the QM and the MM subsystems are handled using +the GEEP approach as described by Laino et al. :ref:`189 `. +This method of evaluating interactions between the QM and MM subsystems +is a variant of the "electrostatic embedding" scheme. The electrostatic +interactions between the electrons of the QM region and the MM atoms +and between the QM nuclei and the MM atoms are explicitly included into +the Hamiltonian for the QM subsystem: .. math:: H^{QM/MM} = H^{QM}_e-\sum_i^n\sum_J^M\frac{e^2Q_J}{4\pi\epsilon_0r_{iJ}}+\sum_A^N\sum_J^M\frac{e^2Z_AQ_J}{e\pi\epsilon_0R_{AJ}}, - where :math:`n` and :math:`N` are the number of electrons and nuclei - in the QM region, respectively, and :math:`M` is the number of - charged MM atoms. The first term on the right hand side is the - original electronic Hamiltonian of an isolated QM system. The first - of the double sums is the total electrostatic interaction between the - QM electrons and the MM atoms. The total electrostatic interaction of - the QM nuclei with the MM atoms is given by the second double sum. - Bonded interactions between QM and MM atoms are described at the MM - level by the appropriate force-field terms. Chemical bonds that - connect the two subsystems are capped by a hydrogen atom to complete - the valence of the QM region. The force on this atom, which is - present in the QM region only, is distributed over the two atoms of - the bond. The cap atom is usually referred to as a link atom. - -#. **ONIOM** In the ONIOM approach, the energy and gradients are first - evaluated for the isolated QM subsystem at the desired level of *ab - initio* theory. Subsequently, the energy and gradients of the total - system, including the QM region, are computed using the molecular - mechanics force field and added to the energy and gradients - calculated for the isolated QM subsystem. Finally, in order to - correct for counting the interactions inside the QM region twice, a - molecular mechanics calculation is performed on the isolated QM - subsystem and the energy and gradients are subtracted. This leads to - the following expression for the total QM/MM energy (and gradients - likewise): - - .. math:: - - E_{tot} = E_{I}^{QM} - +E_{I+II}^{MM}-E_{I}^{MM}, - - where the subscripts I and II refer to the QM and MM subsystems, - respectively. The superscripts indicate at what level of theory the - energies are computed. The ONIOM scheme has the advantage that it is - not restricted to a two-layer QM/MM description, but can easily - handle more than two layers, with each layer described at a different - level of theory. +where :math:`n` and :math:`N` are the number of electrons and nuclei +in the QM region, respectively, and :math:`M` is the number of +charged MM atoms. The first term on the right hand side is the +original electronic Hamiltonian of an isolated QM system. The first +of the double sums is the total electrostatic interaction between the +QM electrons and the MM atoms. The total electrostatic interaction of +the QM nuclei with the MM atoms is given by the second double sum. +An important advantage of using the CP2K/GEEP combination is that it allows +evaluation of forces for both QM-QM and QM-MM interactions, +in the case of systems with periodic boundary conditions (PBC). +To avoid double accounting for electrostatic interactions and LJ, +classical MM charge on the QM atoms are zeroed out as well as LJ +interactions between QM-QM atoms are excluded. It should be noted that +LJ interactions between QM-MM atoms are kept and still calculated by |Gromacs|. +Bonded interactions between QM and MM atoms are described at the MM +level by the appropriate force-field terms. All bonds, +consisting of 2 QM atoms, angles and settles containing 2 or 3 QM atoms, +dihedrals containing 3 or 4 QM atoms are excluded from the forcefield +evaluation. Broken chemical bonds between QM and MM subsystems needs to be capped +in the QM calculation. This is done within CP2K by adding a hydrogen atom to +complete the valence of the QM region. The force on this atom, which is present +in the QM region only, is distributed over the two atoms of the bond. +The cap atom is usually referred to as a link atom. Within the interface +all described topology modifications are performed automatically during :ref:`gmx grompp` pre-processing. + +Software prerequisites +^^^^^^^^^^^^^^^^^^^^^^ + +CP2K version 8.1 (or later) should be linked into |Gromacs| as libcp2k. +For a specific installation instructions please follow the :ref:`installing with CP2K` guide. + +Limitations in simulation techniques +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The QM/MM interface limits simulations in two ways. +First, no topology modifications are possible during the simulations in the QM region. +Second, interface completely ignores "B" state parameters in the topology, making +double topology setups impossible, e.g. free-energy perturbation simulations (:ref:`dgimplement`). + +In addition it should be noted that the contribution of forces from QM/MM to the system +virial are not accounted for. The size of the effect on the pressure-coupling algorithm +grows with the total summed force due to QM-MM interactions and might produce artifacts +in simulations with the NPT ensemble. Usage ^^^^^ -QMMM is currently not supported in GROMACS. +QM/MM simulations with CP2K interface are controlled by setting :ref:`mdp` file options and, +in some cases, providing an additional input file for :ref:`gmx grompp` with the ``-qmi`` +command-line option. All options that are related to QM/MM simulations with CP2K +are prefixed with ``qmmm-cp2k``. + +Setting :mdp-value:`qmmm-cp2k-active=true` will trigger a QM/MM simulation using the whole +system as QM part and default parameters for all other options. + +Choosing atoms for QM calculation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The QM part of your system is chosen with a name that corresponds to an atom group +in the index file of |Gromacs| to the :mdp:`qmmm-cp2k-qmgroup` option in :ref:`mdp` file. +The typical QM part should consist of atoms that are interesting from the chemical point of view, +i.e. part of the system where reaction happens. To make computation of the +QM part feasible, it should be small and +as compact as possible in a space. DFT simulations often scale as 3rd order of +the number of atoms in the QM part. This means increasing number of atoms in the QM part +by a factor of 2 will slow down the simulation by a factor of 8. + +In addition user should provide total charge of your QM subsystem with +:mdp:`qmmm-cp2k-qmcharge` option and spin-state (multiplicity) with :mdp:`qmmm-cp2k-qmmultiplicity` +option. + +Supported QM methods +^^^^^^^^^^^^^^^^^^^^ + +The QM method is chosen with :mdp:`qmmm-cp2k-qmmethod` in the :ref:`mdp` file. +Currently the following QM methods are supported: + +#. :mdp-value:`qmmm-cp2k-qmmethod=PBE` - DFT using PBE functional and DZVP-MOLOPT basis set. +#. :mdp-value:`qmmm-cp2k-qmmethod=BLYP` - DFT using BLYP functional and DZVP-MOLOPT basis set. + +That list will be updated with a new methods once they are tested and included into the +interface. + +Providing your own CP2K input file +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In addition it is possible to use custom external CP2K input file with +:mdp-value:`qmmm-cp2k-qmmethod=INPUT` and providing file with +:ref:`gmx grompp` with ``-qmi`` option. The external file will be incorporated into the +:ref:`tpr` file of the simulation and are subject to the following restrictions: + +#. ``RUN_TYPE`` option in the CP2K input should be equal to ``ENERGY_FORCE``. +#. ``CHARGE`` option should be present. +#. ``MULTIPILICTY`` option should be present. +#. ``COORD_FILE_NAME`` option should be present pointing towards :ref:`pdb` file. +#. Both ``CHARGE_EXTENDED TRUE`` and ``COORD_FILE_FORMAT PDB`` options should be present. +#. Incremental includes (``@INCLUDE`` directive) are not allowed in the CP2K input file . + +Changing names of CP2K files +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +During :ref:`gmx mdrun` simulation additional files will be produced with ``.inp``, ``.out`` and +``.pdb``. They contain CP2K input, CP2K output and :ref:`pdb` file with point charges of MM atoms +in the extended beta field. By default all CP2K related files names will be deduced from :ref:`tpr` +simulation file name by adding ``_cp2k`` suffix. In order to change it manually +:mdp:`qmmm-cp2k-qmfilenames` option should be used. + +Output +^^^^^^ + +The energy output file will contain an additional "Quantum En." term. +This is the energy that is added to the system from the QM/MM interactions. +In addition, a file containing CP2K output will appear in the simulation directory +with the ``.out`` extension. + +Future developments +^^^^^^^^^^^^^^^^^^^ + +support of additional DFT methods will be added in the future, as well as semi-empirical and +DFTB description of the QM subsystem will be allowed. Support of the multiple +time-stepping approach to speed-up simulation will be added. Excited state simulations +will be implemented with TD-DFT description of the wavefunction. diff --git a/docs/release-notes/2022/major/features.rst b/docs/release-notes/2022/major/features.rst index 8c982042a2..acf78a2804 100644 --- a/docs/release-notes/2022/major/features.rst +++ b/docs/release-notes/2022/major/features.rst @@ -7,6 +7,18 @@ New and improved features Also, please use the syntax :issue:`number` to reference issues on GitLab, without the a space between the colon and number! +Hybrid Quantum-Classical simulations (QM/MM) with CP2K interface +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +Simulations of chemical reactions pathways can provide an atomistic insight into many +biological and chemical processes. To perform such kind of modelling in complex systems, +that includes solvent and/or proteins Multi-scale Quantum Mechanics / Molecular Mechanics +(QM/MM) approaches are often used. Here we introduce a whole new interface to perform QM/MM +simulations in fully periodic systems using MDModule that couples |Gromacs| with CP2K +quantum chemistry package. This enables hybrid simulations of systems in systems +where chemical reactions occurs. The interface supports most of the simulations techniques +available in |Gromacs| including energy minimization, classical MD and enhanced sampling methods +such as umbrella sampling and accelerated weight histogram method. Transformation pull coordinate for mathematical transformations of pull coordinates """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" -- 2.22.0