Table of Contents

# Matt Watkins for the CP2K-UK summer school 2016

Multiscale modelling - mixed hamiltonians and external fields in CP2K

Matt Watkins

School of Mathematics and Physics, University of Lincoln, UK

https://www.cp2k.org

critcat_logo.jpg

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.

MgO cluster on bulk Ag?320

We could consider

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.

 CP2K_INPUT / FORCE_EVAL / DFT / EXTERNAL_POTENTIAL
(CP2K_INPUT / FORCE_EVAL / DFT / EFIELD (for real time propagation))
 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.

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

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.

 3D system 2D sandwich 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 …

$$ 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$

$$ 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

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.

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

CP2K_INPUT / FORCE_EVAL / QMMM / PERIODIC

section. The default is non periodic embedding.

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

CP2K_INPUT / FORCE_EVAL / QMMM / PERIODIC / MULTIPOLE 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

  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

$$ 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 exercise

Input files

Example setup for KCl that we used here.

We need to define the whole system as normal

  &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

We need a normal section for the QM part

  &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

The QMMM section is

 &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

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

&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

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 exercise

Example - subtractive QM/MM

We can implement very simple subractive QMMM using a mixed force_env that would look schematically like this

&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

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/