Differences

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

--- events:2016_summer_school:qmmm [2016/08/25 09:47] – 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 19: / Line 19: @@
 It is a very common problem that we need to describe an imhomogeneous system, and we may be more interested in some parts than others. However, the rest of the system is still important.
-{{ exercises:2016_summer_school:MgOAg.png |MgO cluster on bulk Ag}}
+{{ exercises:2016_summer_school:MgOAg.png |MgO cluster on bulk Ag?320}}
 We could consider
@@ 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}}
+{{exercises:2016_summer_school:rhodopsin.png?320}}
 typically the active site was surrounded by a finite number of classical point charges
-{{ exercises:2016_summer_school:materials/qmmmcartoon.png}}
+{{ exercises:2016_summer_school:qmmmcartoon.png?320}}
 and the surface terms (boundary of the MM) was just hydrogen terminated, or extra point charges were added to make electrostatics well behaved, or a continuum field model was added.
@@ 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/