Differences

This shows you the differences between two versions of the page.

--- events:2016_summer_school:qmmm [2016/08/25 09:48] – mwatkins
+++ events:2016_summer_school:qmmm [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 9: / Line 9: @@
 School of Mathematics and Physics, University of Lincoln, UK
-http://www.cp2k.org
+https://www.cp2k.org
 {{https://www.cp2k.org/lib/tpl/cp2kwiki/images/logo.png }}
@@ Line 74: / Line 74: @@
 $$
-Periodic electric field uses the Berry phase formalism of the [[http://journals.aps.org/prb/abstract/10.1103/PhysRevB.47.1651|Modern Theory of Polarizablility]] and can be used for periodic systems.
+Periodic electric field uses the Berry phase formalism of the [[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.47.1651|Modern Theory of Polarizablility]] and can be used for periodic systems.
-{{http://pubs.acs.org/appl/literatum/publisher/achs/journals/content/jpclcd/2016/jpclcd.2016.7.issue-14/acs.jpclett.6b01127/20160725/images/medium/jz-2016-01127a_0005.gif}}
+{{https://pubs.acs.org/appl/literatum/publisher/achs/journals/content/jpclcd/2016/jpclcd.2016.7.issue-14/acs.jpclett.6b01127/20160725/images/medium/jz-2016-01127a_0005.gif}}
-[[http://pubs.acs.org/doi/abs/10.1021/acs.jpclett.6b01127|Computing the Kirkwood g-Factor by Combining Constant Maxwell Electric Field and Electric Displacement Simulations: Application to the Dielectric Constant of Liquid Water, Chao Zhang, Jürg Hutter, and Michiel Sprik, J. Phys. Chem. Lett., 2016, 7 (14), pp 2696–2701]]
+[[https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.6b01127|Computing the Kirkwood g-Factor by Combining Constant Maxwell Electric Field and Electric Displacement Simulations: Application to the Dielectric Constant of Liquid Water, Chao Zhang, Jürg Hutter, and Michiel Sprik, J. Phys. Chem. Lett., 2016, 7 (14), pp 2696–2701]]
 ==== QM/MM ====
@@ Line 84: / Line 84: @@
 Well known method is QMMM, one main strand arose from the bio community - aiming to accurately model active sites in proteins
-{{ exercises:2016_summer_school:rhodopsin.png?320}}
+{{exercises:2016_summer_school:rhodopsin.png?320}}
 typically the active site was surrounded by a finite number of classical point charges
@@ Line 96: / Line 96: @@
 An attractive feature of CP2K's QMMM implementation is that it can be fully periodic, or anything from a cluster to a 3D system.
-{{ exercises:2016_summer_school:qmmm_island.png?320}}
+{{exercises:2016_summer_school:qmmm_island.png?320 | 3D system }}
+{{ exercises:2016_summer_school:qmmm_sandwich.png?320|2D sandwich system}}
-D system
+==== QMMM Hamiltonian ====
+Generally CP2K works with an additive QMMM Hamiltonian:
+$$
+E_{tot}(\mathbf{R}_\alpha , \mathbf{R}_a) = E_{QM}(\mathbf{R}_\alpha) + E_{MM}( \mathbf{R}_a) + E_{QMMM}(\mathbf{R}_\alpha , \mathbf{R}_a)
+$$
+Total energy is just the QM energy + the MM energy + the interaction between them.
+Where the system is partitioned into QM atoms, at positions $(\mathbf{R}_\alpha)$ and MM atoms at position $(\mathbf{R}_a)$.
+It is also possible to use 'subtractive schemes' (ONIOMM in Gaussian code for instance):
+$$
+E_{tot}(\mathbf{R}_\alpha , \mathbf{R}_a) = E_{QM}(\mathbf{R}_\alpha) - E_{MM}( \mathbf{R}_\alpha) + E_{MM}(\mathbf{R}_\alpha , \mathbf{R}_a)
+$$
+or for a QM in QM embedding:
+$$
+E_{tot}(\mathbf{R}_\alpha , \mathbf{R}_a) = E_{QM^1}(\mathbf{R}_\alpha) - E_{QM^2}( \mathbf{R}_\alpha) + E_{QM^2}(\mathbf{R}_\alpha , \mathbf{R}_a)
+$$
+=== Additive QMMM in CP2K ===
+$$
+E_{tot}(\mathbf{R}_\alpha , \mathbf{R}_a) = E_{QM}(\mathbf{R}_\alpha) + E_{MM}( \mathbf{R}_a) + E_{QMMM}(\mathbf{R}_\alpha , \mathbf{R}_a)
+$$
+where the coupling term is mainly electrostatic
+$$
+E_{QMMM}(\mathbf{R}_\alpha , \mathbf{R}_a) = \sum_{a \in MM} q_a \int_r \frac{n_{tot} (\mathbf{r})}{\mid \mathbf{r} - \mathbf{R}_a \mid} \text{d}\mathbf{r}
+$$
+where $n_{tot}$ is the total electronic and nuclear charge density of the QM system and $q_a$ is the charge of the MM atom at location $\mathbf{R}_a$
+== Gaussian Expansion of the Electrostatic Potential (GEEP) ==
+As always in CP2K, we try and use Gaussians ...
+  * The point charge MM atoms can be replaced with Gaussian charge distributions
+$$
+n(|\mathbf{r}-\mathbf{R}_a|) = \left( \frac{1}{\sqrt \pi r_{c,a}}\right) exp \left( \frac{|\mathbf{r}-\mathbf{R}_a|^2}{r_{c,a}^2}\right) \Rightarrow
+v_a(\mathbf{r},\mathbf{R}_a) = \frac{erf(\frac{|\mathbf{r}-\mathbf{R}_a|}{r_{c,a}})}{|\mathbf{r}-\mathbf{R}_a|}
+$$
+where the error function is $erf(x) = \frac{2}{\sqrt \pi} \int_0^x e^{-t^2}\text{d}t$
+  * expand the error function as a linear combination of Gaussians with different exponents
+$$
+v_a(\mathbf{r},\mathbf{R}_a) = \frac{erf(\frac{|\mathbf{r}-\mathbf{R}_a|}{r_{c,a}})}{|\mathbf{r}-\mathbf{R}_a|} =
+\sum_{N_g} A_g exp \big(\frac{|\mathbf{r}-\mathbf{R}_a|^2}{r_{c,a}^2} \big) + R_{low} (|\mathbf{r}-\mathbf{R}_a|)
+$$
+the final term $R_{low} (|\mathbf{r}-\mathbf{R}_a|)$ is the residual part of the function not represented by the Gaussians, and should be rather smooth.
+The number of terms in the sum $N_g$ is set by the input variable ''USE_GEEP_LIB''
+{{https://pubs.acs.org/appl/literatum/publisher/achs/journals/content/jctcce/2005/jctcce.2005.1.issue-6/ct050123f/production/images/medium/ct050123ff00001.gif}}
+[[https://pubs.acs.org/doi/full/10.1021/ct050123f|An Efficient Real Space Multigrid QM/MM Electrostatic Coupling, Teodoro Laino, Fawzi Mohamed , Alessandro Laio , and Michele Parrinello, J. Chem. Theory Comput., 2005, 1 (6), pp 1176–1184]]
+== Short range electrostatic coupling - collocating the potential ==
+<code>
+  METHOD QMMMM
+  @include QS.inc
+  @include MM.inc
+  &QMMM
+    #this defines the QS cell in the QMMM calc
+    &CELL
+      ABC 12.6 15.0 12.6
+      PERIODIC XZ
+    &END CELL
+    ECOUPL GAUSS # use GEEP method
+    NOCOMPATIBILITY
+    USE_GEEP_LIB 6  # use GEEP method
+</code>
+{{https://pubs.acs.org/appl/literatum/publisher/achs/journals/content/jctcce/2005/jctcce.2005.1.issue-6/ct050123f/production/images/medium/ct050123ff00002.gif}}
+The short range part is put onto grids in much the same manner as in the GPW method.
+[[https://pubs.acs.org/doi/full/10.1021/ct050123f|An Efficient Real Space Multigrid QM/MM Electrostatic Coupling, Teodoro Laino, Fawzi Mohamed , Alessandro Laio , and Michele Parrinello, J. Chem. Theory Comput., 2005, 1 (6), pp 1176–1184]]
+== Periodic embedding ==
+this is the confusing bit to work with. Must be activated with the
+[[https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/QMMM/PERIODIC.html|CP2K_INPUT / FORCE_EVAL / QMMM / PERIODIC]]
+section. The default is non periodic embedding.
+{{https://pubs.acs.org/appl/literatum/publisher/achs/journals/content/jctcce/2006/jctcce.2006.2.issue-5/ct6001169/production/images/medium/ct6001169f00001.gif}}
+[[https://pubs.acs.org/doi/abs/10.1021/ct6001169|An Efficient Linear-Scaling Electrostatic Coupling for Treating Periodic Boundary Conditions in QM/MM Simulations, Teodoro Laino, Fawzi Mohamed, A. Laio, M. Parrinello, JCTC, 2, 1370 (2006)]]
+== Long range coupling ==
+has two components
+  * QM/QM interactions (probably small and maybe not critical)
+[[https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/QMMM/PERIODIC.html | CP2K_INPUT / FORCE_EVAL / QMMM / PERIODIC / MULTIPOLE]]
+this is on be default if the periodic keyword is activated
+  * Residual potential $R_{low}$ is long ranged and can be periodically summed using Ewald techniques. This is on be default if the periodic keyword is activated.
+== Coulomb coupling ==
+Alternatively for Semi Empirical Hamiltonians or DFTB it is possible to use "coulomb" embedding.
+Here the field from the classical ions acts on the Gaussian basis functions, much like the efield talked about earlier
+<code>
+  METHOD QMMMM
+  @include QS.inc
+  @include MM.inc
+  &QMMM
+    #this defines the QS cell in the QMMM calc
+    &CELL
+      ABC 12.6 15.0 12.6
+      PERIODIC XZ
+    &END CELL
+    ECOUPL COULOMB # use classical point charge method
+</code>
+$$
+V_{ab,QMMM}^{coulomb} = -\sum_{I \in MM atoms} \big{<} \phi_a \big{|} \frac{Z_I}{\mid \mathbf{R_I}-\mathbf{r} \mid} \big{|} \phi_b \big{>}
+$$
+see this [[https://www.cp2k.org/exercises:2015_cecam_tutorial:urea|exercise]]
+===== Input files =====
+Example setup for KCl that we used [[https://onlinelibrary.wiley.com/doi/10.1002/jcc.23904/full| here]].
+{{https://onlinelibrary.wiley.com/store/10.1002/jcc.23954/asset/image_n/jcc23954-toc-0001.png?v=1&s=6f79087c34654bcc15eda422e9c03888b4ee9550}}
+We need to define the whole system as normal
+<code>
+  &SUBSYS
+    #this defines the cell of the whole system
+    #must be orthorhombic, I think
+    &CELL
+      ABC 12.6 100.0 12.6
+    &END CELL
+    &TOPOLOGY
+      COORD_FILE_NAME kcl.xyz
+      COORD_FILE_FORMAT XYZ
+      &GENERATE
+         &ISOLATED_ATOMS
+         #ignores bonds dihedrals etc in classical part
+            LIST 1..48
+         &END
+      &END
+    &END
+    &KIND K
+      ELEMENT K
+      BASIS_SET DZVP-MOLOPT-SR-GTH
+      POTENTIAL GTH-PBE-q9
+    &END KIND
+    &KIND Cl
+      BASIS_SET DZVP-MOLOPT-GTH
+      POTENTIAL GTH-PBE-q7
+    &END
+  &END SUBSYS
+</code>
+We need a normal section for the QM part
+<code>
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_MOLOPT
+    POTENTIAL_FILE_NAME GTH_POTENTIALS
+    &MGRID
+      COMMENSURATE # this keyword is required for QMMM with GEEP
+      CUTOFF 150
+    &END MGRID
+    &QS
+      EPS_DEFAULT 1.0E-12
+    &END QS
+    &SCF
+      EPS_SCF 1.0E-06
+      MAX_SCF 26
+      SCF_GUESS RESTART
+      &OT
+        MINIMIZER CG
+        PRECONDITIONER FULL_SINGLE_INVERSE
+      &END OT
+      &OUTER_SCF
+        EPS_SCF 1.0E-05
+      &END OUTER_SCF
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE
+      &END XC_FUNCTIONAL
+    &END XC
+    &PRINT
+       &MO_CUBES
+           NLUMO 10
+           WRITE_CUBE F
+       &END MO_CUBES
+    &END PRINT
+  &END DFT
+<\code>
+A MM section
+<code>
+ &MM
+    &FORCEFIELD
+      &CHARGE
+         ATOM K
+         CHARGE 1.0
+      &END CHARGE
+      &CHARGE
+         ATOM Cl
+         CHARGE -1.0
+      &END CHARGE
+      &NONBONDED
+        &WILLIAMS
+          atoms K   Cl
+          A [eV] 4117.9
+          B [angstrom^-1] 3.2808
+          C [eV*angstrom^6] 0.0
+          RCUT [angstrom] 3.0
+        &END WILLIAMS
+        &WILLIAMS
+          atoms Cl  Cl
+          A [eV] 1227.2
+          B [angstrom^-1] 3.1114
+          C [eV*angstrom^6] 124.0
+          RCUT [angstrom] 3.0
+        &END WILLIAMS
+        &WILLIAMS
+          atoms K   K
+          A [eV] 3796.9
+          B [angstrom^-1] 3.84172
+          C [eV*angstrom^6] 124.0
+          RCUT [angstrom] 3.0
+        &END WILLIAMS
+      &END NONBONDED
+    &END FORCEFIELD
+    &POISSON
+      &EWALD
+        EWALD_TYPE spme
+        ALPHA .44
+        GMAX  40
+      &END EWALD
+    &END POISSON
+  &END MM
+</code>
+The QMMM section is
+<code>
+ &QMMM
+    #this defines the QS cell in the QMMM calc
+    &CELL
+      ABC 12.6 15.0 12.6
+      PERIODIC XZ
+    &END CELL
+    ECOUPL GAUSS # use GEEP method
+    NOCOMPATIBILITY
+    USE_GEEP_LIB 6  # use GEEP method
+    &PERIODIC # apply periodic potential
+      #in this case QM box = MM box in XZ so turn
+      #off coupling/recoupling of the QM multipole
+      &MULTIPOLE OFF
+      &END
+    &END PERIODIC
+    #these are just the ionic radii of K Cl
+    #but should be treated as parameters in general
+    #fit to some physical property
+    &MM_KIND K
+      RADIUS 1.52
+    &END MM_KIND
+    &MM_KIND Cl
+      RADIUS 1.67
+    &END MM_KIND
+    #define the model
+    &QM_KIND K
+      MM_INDEX 25..32 41..48
+    &END QM_KIND
+    &QM_KIND Cl
+      MM_INDEX 17..24 33..40
+    &END QM_KIND
+  &END QMMM
+</code>
+Note the CELL in the QMMM section is not the same size as in the main `&SUBSYS` section. We only need a cell large enough to contain the electron density of the QM region.
+==== Multiple force environments ====
+it is possible to create rather interesting effects by combining results from several calculations in some way:
+For instance there is an example in `$CP2K/cp2k/tests/QS/regtest-meta/H2O-IP-meta.inp` that performs metadynamics using the ionisation energy of a molecule as a collective variable.
+A mixed calculation in CP2K will have multiple `FORCE_EVAL` sections
+<code>
+&MULTIPLE_FORCE_EVALS
+  FORCE_EVAL_ORDER 2 3
+&END
+&FORCE_EVAL
+ METHOD MIXED
+ &MIXED
+   MIXING_TYPE GENMIX
+   &GENERIC
+     MIXING_FUNCTION X+Y
+     VARIABLES X Y
+   &END GENERIC
+ &END
+&END FORCE_EVAL
+&FORCE_EVAL
+  METHOD FIST
+&END FORCE_EVAL
+&FORCE_EVAL
+  METHOD QS
+&END FORCE_EVAL
+</code>
+The default is to have a mapping 1-1 between atom index (i.e. all force_eval share the same geometrical structure).
+This can be changed by providing a mapping between atoms in the different force_evals.
+See this [[https://www.cp2k.org/exercises:2015_cecam_tutorial:neb|exercise]]
+=== Example - subtractive QM/MM ===
+We can implement very simple subractive QMMM using a mixed force_env that would look schematically like this
+<code>
+&MULTIPLE_FORCE_EVALS
+  FORCE_EVAL_ORDER 2 3
+&END
+&FORCE_EVAL
+ METHOD MIXED
+ &MIXED
+   MIXING_TYPE GENMIX
+   &GENERIC
+     MIXING_FUNCTION X+Y-Z
+     VARIABLES X Y Z
+   &END GENERIC
+ &END
+&END FORCE_EVAL
+&FORCE_EVAL
+  METHOD FIST
+&END FORCE_EVAL
+&FORCE_EVAL
+  METHOD QS
+&END FORCE_EVAL
+&FORCE_EVAL
+  METHOD FIST
+&END FORCE_EVAL
+</code>
+==== Task farming ====
+A final note is that CP2K has quite reasonable task farming capability
+There are some examples in the test directories $CP2K/cp2k/tests/