CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
howto:cdft

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
howto:cdft [2018/11/01 08:45] – [Using the CDFT module] nholmberhowto:cdft [2024/01/03 13:20] (current) – oschuett
Line 1: Line 1:
-====== Constrained DFT ====== +This page has been moved to: https://manual.cp2k.org/trunk/methods/dft/constrained.html
- +
-This tutorial demonstrates how to perform constrained DFT (CDFT) simulations with CP2K. No previous experience with CDFT simulations is required to complete this tutorial. However, a good understanding of running DFT simulations with CP2K/QS is recommended before proceeding. +
- +
-This tutorial is divided as follows. First, a brief overview of the underlying theory behind CDFT will be presented. Typical applications where CDFT simulations have been used will also be highlighted. The CDFT implementation in CP2K will then be described in detail in the next section with the aid of realistic example calculations. The last part of this tutorial covers how to calculate properties involving multiple CDFT states. +
- +
-<note important> +
- +
-CP2K version 5.1 or higher is needed to perform CDFT simulations. +
- +
-</note> +
- +
----- +
-===== CDFT in summary ===== +
- +
-CDFT is a tool for constructing charge and/or spin localized states. Such localized states are needed in a number of applications. These include for example the following +
-  * studying charge transfer phenomena and calculating electronic couplings (e.g. using the Marcus theory approach) +
-  * correcting spurious charge delocalization due to self-interaction error +
-  * parametrizing model Hamiltonians (e.g. the Heisenberg spin Hamiltonian) +
- +
-A more exhaustive list of potential applications has been presented in this [[doi>10.1021/cr200148b | review article]]. +
- +
-The charge and spin localized states are created by enforcing electron and spin density localization within atom centered regions of space. The relevant theory has been derived by Wu and Van Voorhis in a series of key papers: [[doi>10.1103/PhysRevA.72.024502 | paper 1]],  [[doi>10.1063/1.2360263 | paper 2]],  [[doi>10.1063/1.2360263 | paper 3]]. Further useful references can be found in the aforementioned review article. The CDFT implementation of CP2K has been throughly described in this [[doi>10.1021/acs.jctc.6b01085 | publication]]. +
- +
-In this tutorial, only the main results needed to understand what is happening during a CDFT simulation will be summarized. The charge/spin localized states can be generated by augmenting the Kohn-Sham energy functional, $E_\mathrm{KS}$, by additional constraint potentials +
- +
-\begin{equation} +
-E_\mathrm{CDFT}[\rho, \vec{\lambda}] = \max_\vec{\lambda} \min_\rho \left( E_\mathrm{KS}[\rho] + \sum_c \lambda_c \left[\sum_{i = \uparrow, \downarrow} \int w_c^i(\mathbf{r})\rho^i(\mathbf{r})d\mathbf{r} - N_c \right] \right) +
-\end{equation} +
- +
-where $\vec\lambda = [\lambda_1, \lambda_2, \cdots]^T$ are the constraint Lagrangian multipliers ("strength(s) of the constraint potential(s)"), $w^i(\mathbf{r})$ is an atom centered weight function, and $N_c$ is the target value of the constraint. Multiple constraints can be included in a CDFT simulation (the sum over $c$ above). The weight function is constructed as a normalized sum over select constraint atoms $\mathcal{C}$ +
- +
-\begin{equation} +
-w^i(\mathbf{r}) = \frac{\sum_{j \in \mathcal{C}}c_jP_j(\mathbf{r})}{\sum_{j \in \mathcal{N}}P_j(\mathbf{r})} +
-\end{equation} +
-  +
-where $c_j$ are atomic coefficients which determine how each atom is included in the constraint (more on this later), $P_j$ is the so-called cell function which determines the volume occupied by atom $j$ according to some population analysis method, and $\mathcal{N}$ is the set of all atoms in a system. Currently, only the Becke partitioning scheme is fully supported in CP2K, which will be elaborated in the following section. Different types of constraints can be constructed by modifying the weight function according to the following conventions  +
-  * charge density constraint ($\rho^\uparrow + \rho^\downarrow$): $w^\uparrow = w^\downarrow = w$ +
-  * magnetization density constraint ($\rho^\uparrow - \rho^\downarrow$): $w^\uparrow = -w^\downarrow = w$ +
-  * spin specific constraint ($\rho^{\uparrow/\downarrow}$): $w^{\uparrow/\downarrow} = w, w^{\downarrow/\uparrow} = 0$ +
- +
-When CDFT is used in a molecular dynamics or a geometry optimization simulation, additional force terms arising from the constraints are calculated +
- +
-\begin{equation} +
-\mathbf{F}_{c,i} = -\lambda_c \int \frac{\partial w(\mathbf{r})}{\partial \mathbf{R}_i}\rho(\mathbf{r})d\mathbf{r} +
-\end{equation}  +
- +
-The CDFT energy expression, $E_\mathrm{CDFT}$, is solved self-consistently using a two-tiered approach: one external optimization loop for the constraints, and an inner loop to converge the electronic structure. In practice, three SCF loops are needed to integrate CDFT with the OT method, which uses its own outer loop to reset the OT preconditioner. This process has been schematically illustrated in Figure 1.  +
- +
- +
-{{howto:cdftscf.png | Schematic of the CDFT SCF procedure}} +
-**Figure 1.** Schematic of the CDFT SCF procedure. The constraint Lagrangians $\vec\lambda$ are first optimized in the outer CDFT loop, their values are subsequently fixed, and the electron density corresponding to these fixed values is solved like in traditional CP2K DFT simulations. The control is then returned to the outer CDFT loop where convergence of the constraints is checked. This iteration process is repeated until convergence is achieved or until the number of maximum CDFT SCF steps is reached. The structure of the CDFT loop will be further described below.  +
- +
-By definition, all constraints are satisfied when +
- +
-\begin{equation} +
-\vec c(\vec\lambda) = \left[ \sum_{i = \uparrow, \downarrow} \int w_1^i(\mathbf{r})\rho^i(\mathbf{r})d\mathbf{r} - N_1, \cdots \right]^T = \vec 0 +
-\end{equation} +
- +
-The constraint Lagrangian multipliers $\vec\lambda$ can therefore be optimized by minimizing the constraint error expression $\max |\vec c(\vec\lambda)|$ until the largest element decreases below a threshold $\varepsilon$. Root-finding algorithms are used to optimize $\lambda$. For Newton and quasi-Newton class optimizers, a new guess for $\vec\lambda$ at step $n$ is generated according to the following iteration formula +
- +
-\begin{equation} +
-\vec\lambda_n = \vec\lambda_{n-1} - \alpha \mathbf{J}_n^{-1}\vec c(\vec\lambda_{n-1}) +
-\end{equation} +
- +
-where $\alpha \in (0, 1]$ is the step size and $\mathbf{J}^{-1}$ is the inverse Jacobian matrix. The step size $\alpha$ can be fixed or its value can be optimized with backtracking line search, where the value of $\lambda$ is successively reduced if it decreases the constraint error function. The Jacobian matrix is approximated with finite differences, e.g. using a first order forward difference stencil, by perturbing each element of $\vec\lambda$ slightly and re-evaluating the value of $\vec c$ self-consistently (an SCF energy optimization) +
- +
-\begin{equation} +
-\mathbf{J}_{ij} = \frac{\partial \vec c_i(\vec\lambda)}{\partial \lambda_j} \approx \frac{\vec c_i(\vec\lambda+\vec\delta_j)-\vec c_i(\vec\lambda)}{\left|\vec\delta_j\right|} +
-\end{equation}  +
- +
-where $\vec\delta_j$ is a small perturbation of the $j$th component of $\vec\lambda$.  +
- +
----- +
-===== Using the CDFT module ===== +
- +
-The input section [[inp>FORCE_EVAL/DFT/QS/CDFT]] is used to set up a CDFT simulation. A brief description of this input section will be given in the next two subsections. Subsequently, various aspects of running CDFT simulations will be explored through example calculations. +
- +
-==== Defining CDFT SCF parameters  ===== +
- +
-Settings for the CDFT SCF loop are controlled by the input section [[inp>FORCE_EVAL/DFT/QS/CDFT]]. An example of a typical CDFT input is given below. These parameter selections should be suitable for most systems. +
- +
-<note important>The CDFT input structure was altered slightly in CP2K version 6.1. The changes are indicated below. Input files for subsequent example calculations are provided for both versions 5.1 and 6.1, though using the latest available version is always recommended. </note> +
- +
-<code cp2k> +
-&QS +
-  ... +
-  &BECKE_CONSTRAINT +
-    ... +
-  &END BECKE_CONSTRAINT +
-  ! CDFT loop settings +
-  &CDFT +
-    TYPE_OF_CONSTRAINT BECKE +
-    &OUTER_SCF ON +
-      TYPE BECKE_CONSTRAINT +
-      EXTRAPOLATION_ORDER 2 +
-      MAX_SCF 10 +
-      ! Convergence threshold +
-      EPS_SCF 1.0E-3 +
-      ! Optimizer selection: now Newton's method with backtracking line search +
-      OPTIMIZER NEWTON_LS +
-      ! Optimizer step size +
-      STEP_SIZE -1.0 +
-      ! Note that the section CDFT_OPT exists in CP2K version >= 6.1 +
-      ! Remove section for CP2K version 5.1 (keywords are unchanged) +
-      &CDFT_OPT ON +
-        ! Line search settings +
-        MAX_LS 5 +
-        CONTINUE_LS +
-        FACTOR_LS 0.5 +
-        ! Finite difference settings for calculation of Jacobian matrix +
-        JACOBIAN_STEP 1.0E-2 +
-        JACOBIAN_FREQ 1 1 +
-        JACOBIAN_TYPE FD1 +
-        JACOBIAN_RESTART FALSE +
-      &END CDFT_OPT +
-    &END OUTER_SCF +
-  &END CDFT +
-&END QS +
-</code> +
- +
-The structure of this input section is quite straightforward. The keyword [[inp>FORCE_EVAL/DFT/QS/CDFT/OUTER_SCF#EPS_SCF|EPS_SCF]] defines the CDFT constraint convergence threshold $\varepsilon$ and [[inp>FORCE_EVAL/DFT/QS/CDFT/OUTER_SCF#OPTIMIZER|OPTIMIZER]] selects the CDFT optimizer. Using Newton or quasi-Newton optimizers (Broyden methods) is recommended for most applications. These optimizers accept additional control settings that define how the Jacobian matrix is calculated (keywords JACOBIAN_*) and how to optimize the step size $\alpha$ (keywords *_LS). These keywords are available in the [[inp>FORCE_EVAL/DFT/QS/CDFT/OUTER_SCF/CDFT_OPT|CDFT_OPT]] section. MD simulations with a single constraint might benefit from using the bisect optimizer, which avoids building the Jacobian matrix, in case a considerable amount of the total time per MD step is spent in building the Jacobian. Notice, however, that the frequency of Jacobian rebuilds [[inp>FORCE_EVAL/DFT/QS/CDFT/OUTER_SCF/CDFT_OPT#JACOBIAN_FREQ|JACOBIAN_FREQ]] can be controlled on a per MD step and per CDFT SCF step basis. The Broyden optimizers require less frequent rebuilds of the Jacobian matrix because the matrix is [[https://en.wikipedia.org/wiki/Broyden%27s_method|rank-one updated]] every iteration, although the stability of the method with respect to the rebuild frequency needs to be carefully studied. +
- +
-Above, for instance, the Jacobian is explicitly calculated every CDFT SCF iteration and MD step by perturbing each constraint Lagragian using a first order forward difference stencil with a step size of $10^{-2}$. The Newton step size is optimized with backtracking line search using the update formula $\alpha_n = 0.5*\alpha_{n-1}$ for a maximum of 5 steps as long as the CDFT constraint error decreases. +
- +
-==== Defining constraints ===== +
- +
-As explained in the previous section, the CDFT module in CP2K currently only supports using a Becke population based constraint, which is controlled through the section [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT|BECKE_CONSTRAINT]]. The Becke density partitioning method can be considered as a smoothed Voronoi scheme. In Voronoi partitioning, the volume occupied by each atom is the set of real space grid points $\mathbf{r}$ which are closer to that particular atom than to any other atom in the system. An example Voronoi diagram is given below in Figure 2. The line segments in this figure define real space points which are equidistant from two atoms, while vertices correspond to grid points which are equidistant from three or more atoms. The Becke cell function $P_i$ is overlayed on top of the Voronoi diagram and it decays smoothly from 1 to 0 across the Voronoi polyhedron boundary. Using a smooth density partitioning function improves numerical stability in simulations. +
- +
-{{ howto:cdft-becke-atomicsize.png?400}}{{ howto:cdft-becke.png?400 }}\\ +
-**Figure 2.** Comparison of the Voronoi (lines) and Becke partitioning (contours) schemes. At left, the Becke partitioning is performed without atomic size information. At right, the size of the red atom is 30 % larger than the black atoms, and the contours of the red atom extend farther than without atomic size adjustments.  +
- +
-The Voronoi and, by extension, the Becke partitioning methods treat each element equally. This leads to unphysical partial charges in most systems. For example, the Becke scheme predicts a positive charge on oxygen and a negative charge on hydrogen in water (see examples for input files). This problem can be remedied by accounting for atomic radii during the partitioning. This behavior is activated by the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#ADJUST_SIZE|ADJUST_SIZE]] and the atomic radii are defined with the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#ATOMIC_RADII|ATOMIC_RADII]]. The atomic radii should be set to values that reflect the system under simulation, e.g. using additive covalent radii for covalent molecules or Shannon's ionic radii for ionic compounds. An example on how atomic size adjustments affect the Becke cell functions has been visualized above in Figure 2 at right, where the size of the red atom is set to a value 30 % larger than the black atoms causing the red atom's contours to extend farther than without atomic size adjustments. +
- +
-The algorithmic implementation of the Becke density partitioning method has been detailed [[doi>10.1021/acs.jctc.6b01085|here]]. In brief, this involves iterating over each atom pair permutation $\{\mathbf{R}_i, \mathbf{R}_j\}, j\neq i$ at every real space grid point $\mathbf{r}$. This leads to a poor scaling with respect to the system size (cell size and planewave cutoff) and the number of atoms within the system, and is particularly troublesome for  solvated system simulations. The computational cost of the Becke method can be considerably decreased by noting that only the grid points within a cutoff distance $R_{cutoff}$ ([[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#CUTOFF_TYPE|CUTOFF_TYPE]]) of atoms involved in constraints actually need to be considered. The other grid points can be efficiently screened with constraint atom centered spherical Gaussian functions, activated by the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#CAVITY_CONFINE|CAVITY_CONFINE]] and controlled by other keywords of the form CAVITY_*. The exact details of this confinement scheme and why it can be used are explained in the implementation paper. +
- +
-The actual constraints are defined in the section [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT/ATOM_GROUP|ATOM_GROUP]]. Each repetition of this section defines a new constraint. The constraint atoms are selected with the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT/ATOM_GROUP#ATOMS|ATOMS]] and the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT/ATOM_GROUP#COEFF|COEFF]] determines how the atoms are summed up to form the constraint. Usually all coefficients are set to +1, but mixing +1 and -1 coefficients would define the constraint as the difference between two groups of atoms. The keywords [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#TARGET|TARGET]] and [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT#STRENGTH|STRENGTH]] define the constraint target values and the initial constraint strengths $\vec\lambda$, respectively. The constraint target value should be the desired number of valence electrons on the constraint atoms, suitably multiplied by atomic coefficients in case a relative constraint between two atom groups has been used. The constaint type is selected with the keyword [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT/ATOM_GROUP#CONSTRAINT_TYPE|CONSTRAINT_TYPE]]. It is also possible to use fragment based constraints [[inp>FORCE_EVAL/DFT/QS/BECKE_CONSTRAINT/ATOM_GROUP#FRAGMENT_CONSTRAINT|FRAGMENT_CONSTRAINT]], in which case the constraint target value is calculated from the superposition of isolated fragment densities according to the scheme in Figure 3.  +
- +
-{{ howto:cdft-fragment-constraint.png?350 }} +
-**Figure 3.** Using a fragment based Becke constraint. The system is first divided into two fragments with atomic positions fixed in the same configuration as in the full system. The electron and spin densities of the fragment systems are then saved to cube files and subsequently used as input files for the CDFT calculation, where the constraint target value is calculated from the superimposed fragment densities.  +
- +
-An example of a Becke constraint input section is given below. This choice of parameters should be reasonable for most systems, ignoring the atomic radii and constraint definitions which are system dependent. Decreasing the partitioning cutoff might be useful for solvated system MD simulations, but extensive testing is always necessary before starting production simulations. +
- +
-<code cp2k> +
-&QS +
-  ... +
-  &BECKE_CONSTRAINT +
-    ! Take atomic radii into account? +
-    ADJUST_SIZE     FALSE +
-    ATOMIC_RADII    0.63 0.32 +
-    ! Compute Becke charges? +
-    ATOMIC_CHARGES  TRUE +
-    ! Constraint strength and target values +
-    ! Give one value per constraint +
-    STRENGTH        ${BECKE_STR} +
-    TARGET          ${BECKE_TARGET} +
-    ! Cutoff scheme +
-    CUTOFF_TYPE     ELEMENT +
-    ELEMENT_CUTOFF  6.0 +
-    ! Perform Becke partitioning only within the space +
-    ! spanned by constraint atom centered spherical Gaussians +
-    ! (reduces cost for solvated systems) +
-    CAVITY_CONFINE  TRUE +
-    CAVITY_SHAPE    VDW +
-    EPS_CAVITY      1.0E-7 +
-    IN_MEMORY       TRUE +
-    SHOULD_SKIP     TRUE +
-    ! Constraint definitions, each repetition defines a new constraint +
-    &ATOM_GROUP +
-      ATOMS 1 +
-      COEFF 1 +
-      CONSTRAINT_TYPE CHARGE +
-    &END ATOM_GROUP +
-    ! No constraint applied but calculate charges +
-    &DUMMY_ATOMS +
-      ATOMS 2 +
-    &END DUMMY_ATOMS +
-    ! Print information about CDFT calculation +
-    &PROGRAM_RUN_INFO ON +
-      &EACH +
-        QS_SCF 1 +
-      &END EACH +
-      COMMON_ITERATION_LEVELS 2 +
-      ADD_LAST NUMERIC +
-      FILENAME ./${NAME} +
-    &END PROGRAM_RUN_INFO +
-  &END BECKE_CONSTRAINT +
-&END QS +
-</code> +
- +
-==== Selected examples  ===== +
-=== Zn dimer cation === +
- +
-In this example, we will perform two CDFT simulations for the Zn dimer cation $\mathrm{Zn}_2^+$. As the distance $R$ between the two atoms is increased, the excess charge in the system should localize onto one of the Zn atoms forming $\mathrm{Zn}^+ + \mathrm{Zn}$. With standard GGA and hybrid functionals, such as PBE, PBE0, BLYP or B3LYP, this is however not the case. The excess charge will instead be equally shared among the two atoms $\mathrm{Zn}^{0.5+} + \mathrm{Zn}^{0.5+}$ regardless of separation. We can force the charge to localize onto one of these atoms with CDFT. +
- +
-You can download the input files from {{:howto:cdft-tutorial-zn.zip|here}}. Unzip the folder and execute the file ''energy.bash'' (use flag ''-h'' for usage instructions) to run a standard DFT calculation, followed by two CDFT simulations where the first Zn atom is constrained to charges +1 and 0, respectively. The script will also run a mixed CDFT calculation, which will be analyzed in a [[cdft#exampleelectronic_coupling_of_zn_cation_dimer|subsequent section]].  Modify paths to basis sets and other CP2K data files in the file ''dft-common-params.inc'' in case they are in non-standard paths relative to the CP2K binary. The calculations will take a while to run. In the meanwhile, you can +
- +
-  * Check the calculated partial charges in the standard DFT output file +
-  * Study the script file to understand how it works +
-  * Study the generated output files from the CDFT simulations +
- +
-Assuming the CDFT simulations have finished, the following files were generated +
-  * *.out +
-    * This is the standard CP2K output file and now also contains all the output from the CDFT SCF iterations. Each iteration step starts a new DFT energy optimization using new values of the constraint Lagrangian multipliers $\vec{\lambda}$, as generated by the selected CDFT optimizer. +
-    * An example output from the end of one the CDFT simulations is provided below.  +
-    * Information about the CDFT SCF iteration process and the constraints and their convergence is printed alongside the usual CP2K SCF iteration information. +
-    * With CDFT optimizers that support backtracking line search, the density optimization process at each CDFT SCF step is restarted from the optimized constraint strength and density obtained previously during line search. Consequently, when the line search was successful (i.e. the energy optimization converged), the density optimization should converge terminate in exactly 1 step each CDFT SCF iteration (refer to the output files, as this should now be the case).    +
-  * *-LineSearch.out +
-    * The progress of the optimization of the Newton step size $\alpha$ using backtracking line search is reported in this file. Observe that the step size is halved on each iteration if the CDFT constraint error ("Deviation from target" in the output) decreases. +
-  * *.cdftLog +
-    * CDFT parameters (atomic radii, constraint definitions and cutoffs) and Becke partial charges are printed in these files. +
-  * *-JacobianInfo.out +
-    * The output from the calculation of the Jacobian matrix $\mathbf{J}$ is reported in this file. +
-  * *.inverseJacobian +
-    * This is a restart file for the inverse Jacobian matrix. +
- +
-Using the above list, study the generated output files to understand how the CDFT SCF loop is integrated with the standard CP2K DFT SCF process. +
- +
-<file> +
-  SCF WAVEFUNCTION OPTIMIZATION +
- +
-  ----------------------------------- OT --------------------------------------- +
-  Minimizer      : DIIS                : direct inversion +
-                                         in the iterative subspace +
-                                         using   7 DIIS vectors +
-                                         safer DIIS on +
-  Preconditioner : FULL_ALL            : diagonalization, state selective +
-  Precond_solver : DEFAULT +
-  stepsize       :    0.15000000                  energy_gap     :    0.08000000 +
-  ortho_irac     : CHOL                           irac_degree        :             4 +
-  max_irac       :            50                  eps_irac           :   0.10000E-09 +
-  eps_irac_switch:   0.10000E-01                  eps_irac_quick_exit:   0.10000E-04 +
-  on_the_fly_loc     : F +
-  ----------------------------------- OT --------------------------------------- +
- +
-  Step     Update method      Time    Convergence         Total energy    Change +
-  ------------------------------------------------------------------------------ +
-qs_ot_get_orbitals_ref  0: ||P-I||= 0.10493E-10, ortho_irac = POLY +
-qs_ot_ref_poly  1: quick exit! +
-qs_ot_get_orbitals_ref  0: ||P-I||= 0.11959E-12, ortho_irac = POLY +
-qs_ot_ref_poly  1: quick exit! +
-     1 OT DIIS     0.15E+00    2.5     0.00000022      -120.6126709217 -1.21E+02 +
- +
-  *** SCF run converged in     1 steps *** +
- +
- +
-  Electronic density on regular grids:        -22.9999999253        0.0000000747 +
-  Core density on regular grids:               24.0000000000       -0.0000000000 +
-  Total charge density on r-space grids:        1.0000000746 +
-  Total charge density g-space grids:           1.0000000746 +
- +
-  Overlap energy of the core charge distribution:               0.00000000000000 +
-  Self energy of the core charge distribution:               -159.30058829583706 +
-  Core Hamiltonian energy:                                     50.63132167014943 +
-  Hartree energy:                                               5.66516329750697 +
-  Exchange-correlation energy:                                -17.60803221234604 +
-  Dispersion energy:                                           -0.00058223840320 +
- +
-  Total energy:                                              -120.61267092172808 +
- +
-  outer SCF iter =    1 RMS gradient =   0.22E-06 energy =       -120.6126709217 +
-  outer SCF loop converged in   1 iterations or    1 steps +
- +
- +
-  CDFT SCF iter =     5 RMS gradient =   0.13E-03 energy =       -120.6126709217 +
-  CDFT SCF loop converged in   5 iterations or   37 steps +
- +
- +
-  --------------------- Becke constraint information --------------------- +
-  Atomic group                :                                          1 +
-  Type of constraint          :                  Charge density constraint +
-  Target value of constraint  :                            11.000000000000 +
-  Current value of constraint :                            11.000126158558 +
-  Deviation from target       :                                  1.262E-04 +
-  Strength of constraint      :                             0.371415167271 +
-  ------------------------------------------------------------------------ +
-</file> +
- +
-=== Charge transfer energy in water dimer === +
- +
-In this example, we will calculate the charge transfer energy, $-\Delta E_\mathrm{CT}$, of water dimer. This quantity is conveniently defined in CDFT as (see [[doi>10.1039/C6CP07475H|here]])  +
- +
-\begin{equation} +
--\Delta E_\mathrm{CT} = E_\mathrm{CDFT}-E_\mathrm{DFT} +
-\end{equation} +
- +
-where $E_\mathrm{DFT}$ is the DFT total energy of the system, and $E_\mathrm{CDFT}$ is the CDFT energy of the system when charge transfer between the two molecules is prevented. +
- +
-We will calculate the charge transfer energy with four different constraints: default Becke constraint, Becke constraint with atomic size adjustments using covalent radii from [[doi>10.1002/chem.200800987|this publication]], and a fragment based Becke constraint with and without the same atomic size adjustments. The input files can be downloaded from {{:howto:cdft-tutorial-water.zip|here}}. Execute the file ''energy.bash'' to run all the simulations: standard DFT simulations for the full and two fragment systems, and the aforementioned CDFT simulations with different constraints. Modify the include (''.inc'') file if necessary as you would have before. Study the input files while the calculations are running. +
- +
-After the calculations have finished, answer the following questions +
-  - Compare the partial charges of unconstrained PBE water as predicted by the Becke population analysis method with and without atomic size adjustments (Hint: the CDFT calculations were restarted from the DFT wavefunction)  +
-  - How much charge is transferred between the two water molecules according to the different constraint methods? Note that net charges for fragment based constraints are not reported with respect to the core charge, but they can be recovered by post-processing the reported absolute populations. Instead, each atomic charge is referenced to the number of electrons per atom in the system where the isolated densities are superimposed. +
-  - Compare the calculated charge transfer energies to the [[doi>10.1039/C6CP07475H|reference value 1.7 mHa]], calculated at the PBE0/def-QZVP/CDFT level of theory using a different code and constraint. Which is closest to the reference value? Why? (Hint: Look at previous question) +
- +
- +
----- +
-===== Using the mixed CDFT module ===== +
- +
-Additional properties can be calculated from the interactions between CDFT states. In CP2K, these types of simulations are called mixed CDFT simulations because the module leverages the [[inp>FORCE_EVAL/MIXED|MIXED FORCE_EVAL]] type to efficiently treat multiple CDFT states in parallel. Mixed CDFT calculations are useful in a number of applications including +
-  * Calculating charge transfer kinetics parameters +
-  * Performing configuration interaction calculations within the basis of CDFT states +
- +
-In this part of the tutorial, the theoretical basis for mixed CDFT will first be established. The quantities accessible through such simulations will also be introduced. The structure of a mixed CDFT input file will then be discussed. The tutorial is concluded with a walk through of an example calculation.  +
- +
-==== Theoretical basis ==== +
- +
-The theoretical concepts related to mixed CDFT calculations are best introduced through an example. Consider the following one electron transfer processs <chem>X^- + Y -> X + Y^-</chem>. Denote the initial and final states of this reaction as A and B, respectively. Now, according to the Marcus theory of electron transfer, the charge transfer rate of this reaction is given by the rate equation +
- +
-\begin{equation} +
-k_\mathrm{ab}=\frac{2\pi}{\hbar}\frac{\left<\left|\mathbf{H}_\mathrm{ab}\right|^2\right>_T}{\sqrt{4\pi k_bT\xi}}\exp\left(-\frac{(\xi+\Delta A)^2}{4\pi k_bT\xi} \right) +
-\end{equation} +
- +
-where $\Delta A$ is the reaction free energy, $\xi$ is the solvent reorganization energy, and $\left|\mathbf{H}_\mathrm{ab}\right|$ is the electronic coupling. The first two quantities can be obtained from free energy simulations as discussed e.g. in [[doi>10.1021/acs.jctc.6b01085|here]]. The electronic coupling is rigorously defined as the interaction energy between wavefunctions $\Psi$ representing the two reaction states +
- +
-\begin{equation} +
-\mathbf{H}_\mathrm{ab} = \left<\Psi_\mathrm{a}\left| \mathcal{H}\right|\Psi_\mathrm{b}\right> +
-\end{equation} +
- +
-where $\mathcal{H}$ is the many-electron Hamiltonian operator. The usefulness of the electronic coupling quantity is not limited to the Marcus rate equation, but it also a central quantity in other charge transfer theories as well as in CDFT based [[doi>10.1063/1.2800022|configuration interaction]]. +
- +
-The true, interacting many-electron wavefunctions or the Hamiltonian are not available in CDFT simulations. The electronic coupling is instead approximated using the CDFT surrogates +
- +
-\begin{equation} +
-\mathbf{H}_\mathrm{AB} \approx \left<\Phi_\mathrm{A}\left| \mathcal{H}_\mathrm{KS}\right|\Phi_\mathrm{B}\right> = E_\mathrm{B}S_\mathrm{AB}-\sum_c \lambda_c^\mathrm{B} \mathbf{W}_c^\mathrm{AB} \\ +
-\mathbf{H}_\mathrm{BA} \approx \left<\Phi_\mathrm{B}\left| \mathcal{H}_\mathrm{KS}\right|\Phi_\mathrm{A}\right> = E_\mathrm{A}S_\mathrm{BA}-\sum_c \lambda_c^\mathrm{A} \mathbf{W}_c^\mathrm{BA} +
-\end{equation} +
- +
-where $\Phi$ are the CDFT Kohn-Sham determinants, $S_\mathrm{AB}= \left<\Phi_\mathrm{A}|\Phi_\mathrm{B}\right>$, $E_\mathrm{I}$ is the CDFT energy of state $\mathrm{I}$, and $\mathbf{W}_c^\mathrm{AB}$ are the weight function matrices defined by +
- +
-\begin{equation} +
-\mathbf{W}_c^\mathrm{AB} = \left<\Phi_\mathrm{A}\left| w_c^\mathrm{B}(\mathbf{r}) \right|\Phi_\mathrm{B}\right>  +
-\end{equation} +
-  +
- +
-In the above expressions, capital subscripts have been used to emphasize the fact that the CDFT determinants are in general nonorthogonal. The electronic couplings and overlaps are collected into matrices $\mathbf{H}$ and $\mathbf{S}$, respectively. The off-diagonal elements of $\mathbf{H}$ are not symmetric. The matrix is converted to symmetric form by setting +
- +
-\begin{equation} +
-\mathbf{H'}_\mathrm{AB}=\frac{\mathbf{H}_\mathrm{AB}+\mathbf{H}_\mathrm{BA}}{2} = \frac{E_\mathrm{A}+E_\mathrm{B}}{2}S_\mathrm{AB}-\sum_c \left<\Phi_\mathrm{A}\left|\frac{\lambda_c^\mathrm{A}w_c^\mathrm{A}(\mathbf{r})+\lambda_c^\mathrm{B}w_c^\mathrm{B}(\mathbf{r})}{2} \right|\Phi_\mathrm{B}\right>   +
-\end{equation} +
- +
-and setting $\mathbf{H'}_\mathrm{BA}=\mathbf{H'}_\mathrm{AB}$.  +
- +
-The resulting matrix $\mathbf{H'}$ is then orthogonalized to yield the final electronic coupling. The following orthogonalization methods are available: +
-  * Rotate CDFT states to eigenstates of the weight matrix $\mathbf{W}$. This is the default behavior for systems with only one constraint that is identically defined across all CDFT states, not applicable otherwise. +
-  * Löwdin's symmetrical orthogonalization $\mathbf{H} = \mathbf{S}^{-1/2}\mathbf{H}'\mathbf{S}^{-1/2}$. This is the default behavior for systems with multiple constraints and can always be activated with the keyword [[inp>FORCE_EVAL/MIXED/MIXED_CDFT#LOWDIN|LOWDIN]]. +
-  * The so-called wavefunction overlap method where the ground state Kohn-Sham solution is represented as the linear combination of CDFT states, see keyword [[inp>FORCE_EVAL/MIXED/MIXED_CDFT#WFN_OVERLAP|WFN_OVERLAP]]. +
- +
-==== Structure of input file ==== +
- +
-Mixed CDFT calculations are activated through the [[inp>FORCE_EVAL/MIXED|MIXED FORCE_EVAL]] section by setting [[inp>FORCE_EVAL/MIXED#MIXING_TYPE|MIXING_TYPE]] to MIXED_CDFT and providing an appropriate [[inp>FORCE_EVAL/MIXED/MIXED_CDFT|MIXED_CDFT]] input section. The individual CDFT states involved in a mixed CDFT calculation should correspond to different localizations of charge and/or spin in the same system. The constraint definitions do not have to be identical (i.e. defined using the same sets of atoms) in all states as long as the number of constraints is the same.  +
- +
-The CDFT states are included as their own [[inp>FORCE_EVAL|FORCE_EVAL]] sections. It is highly recommended that the CDFT states are first converged in separate simulations, the [[inp>FORCE_EVAL|FORCE_EVAL]] sections from these simulations are then copy-pasted into the mixed CDFT input file, and the mixed CDFT method is used as post-processing analysis tool. The converged wavefunctions and constraint strengths $\vec\lambda$ should be supplied as restart quantities for the mixed CDFT calculation. Using templates and the CP2K ''@include'' and ''@set'' directives is strongly encouraged to keep the mixed CDFT input tidy, see the example input file below +
- +
-<file> +
-&FORCE_EVAL +
-  METHOD MIXED +
-  &MIXED +
-    MIXING_TYPE MIXED_CDFT +
-    NGROUPS  1 +
-    &MIXED_CDFT +
-      ! Calculate mixed CDFT properties every COUPLING step +
-      COUPLING     1 +
-      ! Settings determining how forces are mixed +
-      FORCE_STATES 1 2 +
-      LAMBDA       1.0 +
-      ! Orthogonalize CDFT states using Lowdin's method +
-      ! in addition to standard method +
-      LOWDIN      TRUE +
-      ! Configuration interaction? +
-      CI          FALSE +
-      ! Turn on printing +
-      &PRINT +
-        &PROGRAM_RUN_INFO ON +
-        &END +
-      &END PRINT +
-    &END MIXED_CDFT +
-  &END MIXED +
-  @include subsys.inc +
-&END FORCE_EVAL +
-# Zn+ Zn +
-&FORCE_EVAL +
-  @SET WFN_FILE       ${WFN_FILE_1} +
-  @SET RESTART        ${RESTART_1} +
-  @SET NAME           ${PROJECT_NAME}-state1 +
-  @SET BECKE_TARGET   ${BECKE_TARGET_1} +
-  @SET BECKE_STR      ${BECKE_STR_1} +
-  METHOD QS +
-  @include ${DFT_FILE} +
-&END FORCE_EVAL +
-# Zn Zn+ +
-&FORCE_EVAL +
-  @SET WFN_FILE       ${WFN_FILE_2} +
-  @SET RESTART        ${RESTART_2} +
-  @SET NAME           ${PROJECT_NAME}-state2 +
-  @SET BECKE_TARGET   ${BECKE_TARGET_2} +
-  @SET BECKE_STR      ${BECKE_STR_2} +
-  METHOD QS +
-  @include ${DFT_FILE} +
-&END FORCE_EVAL +
-</file> +
- +
-In the above example input file, a common file ''${DFT_FILE}}'' is used as a template for the [[inp>FORCE_EVAL/DFT|DFT]] subsection. The CDFT state specific constraint settings and the wavefunction filename are passed through variables. The electronic coupling is calculated using the default weight function matrix and Löwdin orthogonalization methods. If molecular dynamics were performed with the above input file, the forces would be mixed according to linear mixing scheme $F = \lambda F_1 + (1-\lambda) F_2$, where the states $F_i$ are selected with the keyword [[inp>FORCE_EVAL/MIXED/MIXED_CDFT#FORCE_STATES|FORCE_STATES]] and the mixing parameter $\lambda$ with [[inp>FORCE_EVAL/MIXED/MIXED_CDFT#LAMBDA|LAMBDA]]. No configuration interaction calculation is performed. The [[inp>FORCE_EVAL/MIXED/MIXED_CDFT|MIXED_CDFT]] section accepts some additional keywords, which have been described in the manual.  +
- +
-The keyword [[inp>FORCE_EVAL/MIXED#NGROUPS|NGROUPS]] is set to 1, which implies that the two CDFT states are treated sequentially utilizing the full set of $N$ MPI processes for the simulation. The CDFT weight function and its gradients are copied from state to state if the constraints definitions are identical in each CDFT state, because construction of these terms might be expensive in large systems. The keyword [[inp>FORCE_EVAL/MIXED#NGROUPS|NGROUPS]] could also be set to 2 in which case the CDFT weight function and gradients would first be built in parallel on $N$ MPI processes, and the two CDFT states would then calculated in parallel with $N/2$ MPI processes. This operating mode is limited to two CDFT states with one identically defined total charge density constraint. This operating mode is an advanced feature which might be useful for large scale MD simulations to save computational wallclock time at the expense of higher CPU core usage. The operating mode should only be used in conjuction with [[inp>FORCE_EVAL/MIXED/MIXED_CDFT#DLB|dynamic load balancing]] if possible, and will likely require tweaking the load balancing parameters. +
-  +
-==== Example: Electronic coupling of Zn cation dimer ====  +
- +
-In this example, we will calculate the electronic coupling for the reaction $\mathrm{Zn}^+ +\mathrm{Zn} \rightleftharpoons \mathrm{Zn}+ \mathrm{Zn}^+$. The initial and final states of this reaction were already converged with CDFT in a [[howto:cdft#zn_dimer_cation|previous section of this tutorial]]. That part of the tutorial must be completed before proceeding. The full input files, in particular the mixed CDFT input file ''energy_mixed.inp'', were also given in that section. +
- +
-The converged CDFT states are used as input for the mixed CDFT calculation. The calculation does not take long to run as a result. The mixed CDFT input file uses template files to keep the input tidy, as was discussed in the previous section. Find and study the corresponding section in the ''energy.bash'' script file to see how variables in the main mixed CDFT template ''energy_mixed.inp'' are initialized. +
- +
-A number of files are generated by the mixed CDFT calculation. The main output from the calculation can be found in the file ''Zn-5A-mixed-cdft.out''. The relevant part of the output is included below. The mixed CDFT analysis is printed after the header lines ''MIXED_CDFT|''. For each unique CDFT state permutation $\{i,j\}, i<j$, the constraint information is first summarized, the overlap and charge transfer energies are printed, and the calculated electronic coupling(s) (and possibly other quantities) are outputted. Here, both tested orthogonalization methods yield an electronic coupling of 5.67 mHa, in agreement with the [[doi>10.1063/1.4867077|5.49 mHa estimate ]] from the more expensive wavefunction based method CASSCF/MRCI+Q. +
- +
-<file> +
- MIXED_CDFT| Activating mixed CDFT calculation +
- MIXED_CDFT| Number of CDFT states:                                            2 +
- MIXED_CDFT| CDFT states calculation mode: serial +
- MIXED_CDFT| Becke constraint is built before the SCF procedure of the first +
-             CDFT state and subsequently copied to other states +
- MIXED_CDFT| Calculating electronic coupling between states:                   T +
- MIXED_CDFT| Calculating electronic coupling reliability metric:               F +
- MIXED_CDFT| Configuration interaction (CDFT-CI) was requested:                F +
- MIXED_CDFT| Block diagonalizing the mixed CDFT Hamiltonian:                   F +
- MIXED_CDFT| Dynamic load balancing enabled:                                   F +
- MIXED_CDFT| Matrix inversions calculated with LU decomposition. +
- +
-  ------------------------- CDFT coupling information -------------------------- +
-  Information at step (fs):                                                 0.00 +
- +
-  ############################################ +
-  ###### CDFT states I =  1 and J =   2 ###### +
-  ############################################ +
-  Atomic group:                                                                1 +
-  Strength of constraint I:                                       0.371415167271 +
-  Strength of constraint J:                                      -0.378315361740 +
-  Final value of constraint I:                                   11.000125488935 +
-  Final value of constraint J:                                   11.999828674539 +
- +
-  Overlap between states I and J:                                 0.030261294466 +
-  Charge transfer energy (J-I) (Hartree):                         0.000739539045 +
- +
-  Diabatic electronic coupling (rotation, mHartree):              5.674875246867 +
-  Diabatic electronic coupling (Lowdin, mHartree):                5.674714192287 +
-  ------------------------------------------------------------------------------ +
- NO FORCE_EVAL section calculated the dipole +
- +
- ENERGY| Total FORCE_EVAL ( MIXED ) energy (a.u.):          -120.612670921735003 +
-</file> +
- +
-Other files created during the execution are related to the individual CDFT states. The ''*-r-1.out, *-r-2.out, ...'' files are the main output files for the CDFT simulations of the studied states. Because preconverged solutions were employed, these CDFT simulations terminate immediately after the first SCF step, and any matrices that are subsequently needed in the mixed CDFT analysis are stored in memory. In the input file, the full project name and CDFT state ID number were stored in the variable ''${NAME}'' on a per state basis. This variable was used to prepend the name of any other output files (e.g. the cdftLog files) that are created during the CDFT simulation of the individual CDFT states. The output files from different CDFT states are therefore straightforward to distinguish. The content of the additional files was discussed in a [[howto:cdft#zn_dimer_cation|previous section]]. +
howto/cdft.1541061955.txt.gz · Last modified: 2020/08/21 10:15 (external edit)