This is an old revision of the document!
Table of Contents
# Matt Watkins for the CP2K-UK summer school 2016
Multiscale modelling - mixed hamiltonians and external fields in CP2K
Flexible modelling
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.
We could consider
- QM/MM a Quantum Mechanical system (Quantum, probably DFT) is coupled to a Molecular Mechanics (classical force field)
- QM/QM a Quantum Mechanical system embedded into a larger Quantum mechanical system described using a cheaper method (smaller basis set, non-hybrid functional, orbital free DFT …)
- External fields, perhaps the influence of the outside world can just be approximated as an external electric field?
Also, we might need to collect a lot of data when dealing with amorphous or disordered systems. Within CP2K there is quite reasonable task farming capability allowing complex jobs to be automated internally. Practically this can mean only a single job being submitted into a queue, rather than thousands…
External fields
If we have our QM system experiencing an applied field we can use several methods in CP2K - this can be confusing.
External fields can be applied in two general ways, reminding us of the GPW method.
- the field / potential can be applied on a grid
CP2K_INPUT / FORCE_EVAL / DFT / EXTERNAL_POTENTIAL (CP2K_INPUT / FORCE_EVAL / DFT / EFIELD (for real time propagation))
- the field can be applied analytically to the GTO
CP2K_INPUT / FORCE_EVAL / DFT / PERIODIC_EFIELD CP2K_INPUT / FORCE_EVAL / DFT / EFIELD
External fields on a grid
This works exactly in the same way that nuclear charge or electrostatics in general is dealt with in GPW.
&EXTERNAL_POTENTIAL
FUNCTION (A/B)*Z
VALUES [eV] 0.2 [angstrom] 1.0
PARAMETERS A B
&END EXTERNAL_POTENTIAL
by default it places an analytical function onto the grid.
Here a constant potential gradient 0.2 eV / Angstrom (electric field) in the $z$ direction.
It is alternatively possible to read the potential on a grid (the cube file read must be a valid cube file with points the fit the finest grid).
Note that this is not periodic! Can only be used with molecules and slabs.
Periodic electric field
The periodic_efield (and efield) keywords apply electric fields that are calculated analytically in the Gaussian basis.
$$ \big{<} \phi_a \big{|} \hat{O} \big{|} \phi_b \big{>} $$
Periodic electric field uses the Berry phase formalism of the Modern Theory of Polarizablility and can be used for periodic systems.
QM/MM
Well known method is QMMM, one main strand arose from the bio community - aiming to accurately model active sites in proteins
typically the active site was surrounded by a finite number of classical point charges
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.
Periodic QMMM
An attractive feature of CP2K's QMMM implementation is that it can be fully periodic, or anything from a cluster to a 3D 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.
[http://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
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
The short range part is put onto grids in much the same manner as in the GPW method.
Periodic embedding
this is the confusing bit to work with. Must be activated with the
CP2K_INPUT / FORCE_EVAL / QMMM / PERIODIC
section. The default is non periodic embedding.
