User Tools

Site Tools


howto:cdft

Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
howto:cdft [2018/07/24 07:38]
nholmber Update to CP2K 6.1
howto:cdft [2018/11/02 14:11] (current)
nholmber [Structure of input file]
Line 21: Line 21:
 A more exhaustive list of potential applications has been presented in this [[doi>​10.1021/​cr200148b | review article]]. 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]].+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 these two papers: ​[[doi>​10.1021/​acs.jctc.6b01085 | paper 1]] and [[doi>​10.1063/​1.5038959 | paper 2]].
  
-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+In this tutorial, only the main theoretical aspects ​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} \begin{equation}
Line 35: Line 35:
 \end{equation} \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 ​+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. 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$   * charge density constraint ($\rho^\uparrow + \rho^\downarrow$):​ $w^\uparrow = w^\downarrow = w$
   * magnetization 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$   * 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+The Becke and Hirshfeld space partitioning schemes can be used as constraint weight functions in CP2K. The main differences between these two constraints will be explained in a subsequent section. Please note that Becke constraints have been tested much more extensively.  
 + 
 +When CDFT is used in a molecular dynamics or a geometry optimization simulation, additional force terms arising from the constraints are calculated ​(currently only implemented for Becke constraints)
  
 \begin{equation} \begin{equation}
Line 75: Line 77:
 ===== Using the CDFT module ===== ===== Using the CDFT module =====
  
-The input sections ​[[inp>​FORCE_EVAL/​DFT/​QS/​CDFT]] ​and [[inp>​FORCE_EVAL/​DFT/​QS/​BECKE_CONSTRAINT]] are used for setting ​up a CDFT simulation. A brief description of these sections ​will be given in the next two subsections. Subsequently,​ various aspects of running CDFT simulations will be explored through example calculations.+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 ​ ===== ==== Defining CDFT SCF parameters ​ =====
Line 81: Line 83:
 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. 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>​+<note important> ​Please note that the CDFT input structure was significantly ​altered in CP2K version ​7.0. As the main difference to older CP2K versions, the Becke constraint section has been merged into the CDFT section. The example calculations ​below include different input files for CP2K versions 5.16.1, and 7.0 (and later). The use of the latest available version is strongly ​recommended. </​note>​
  
 <code cp2k> <code cp2k>
 &QS &QS
   ...   ...
-  &​BECKE_CONSTRAINT 
-    ... 
-  &END BECKE_CONSTRAINT 
   ! CDFT loop settings   ! CDFT loop settings
 +  ! Please note that prior to CP2K version 7.0,
 +  ! Becke constraints were separate from the CDFT section
   &CDFT   &CDFT
     TYPE_OF_CONSTRAINT BECKE     TYPE_OF_CONSTRAINT BECKE
 +    ! Compute CDFT charges?
 +    ATOMIC_CHARGES ​ TRUE
 +    ! Constraint strength and target values
 +    ! Give one value per constraint
 +    STRENGTH ​       ${BECKE_STR}
 +    TARGET ​         ${BECKE_TARGET}
 +    ! 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
 +    ! CDFT convergence and optimizer settings
     &​OUTER_SCF ON     &​OUTER_SCF ON
-      TYPE BECKE_CONSTRAINT+      TYPE CDFT_CONSTRAINT
       EXTRAPOLATION_ORDER 2       EXTRAPOLATION_ORDER 2
       MAX_SCF 10       MAX_SCF 10
       ! Convergence threshold       ! Convergence threshold
       EPS_SCF 1.0E-3       EPS_SCF 1.0E-3
-      ! Optimizer selection: ​now Newton'​s method with backtracking line search+      ! Optimizer selection: 
 +      ! Now Newton'​s method with backtracking line search
       OPTIMIZER NEWTON_LS       OPTIMIZER NEWTON_LS
-      ! Optimizer step size+      ! Optimizer ​(initial) ​step size
       STEP_SIZE -1.0       STEP_SIZE -1.0
       ! Note that the section CDFT_OPT exists in CP2K version >= 6.1       ! Note that the section CDFT_OPT exists in CP2K version >= 6.1
Line 109: Line 128:
         CONTINUE_LS         CONTINUE_LS
         FACTOR_LS 0.5         FACTOR_LS 0.5
-        ! Finite difference settings for calculation of Jacobian matrix+        ! Finite difference settings for Jacobian matrix
         JACOBIAN_STEP 1.0E-2         JACOBIAN_STEP 1.0E-2
         JACOBIAN_FREQ 1 1         JACOBIAN_FREQ 1 1
Line 115: Line 134:
         JACOBIAN_RESTART FALSE         JACOBIAN_RESTART FALSE
       &END CDFT_OPT       &END CDFT_OPT
-    &​END ​OUTER_SCF+    &END 
 +    ! Settigs specific to Becke constraints 
 +    &​BECKE_CONSTRAINT 
 +      ... 
 +    &END BECKE_CONSTRAINT 
 +    ! 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 CDFT   &END CDFT
 &END QS &END QS
 </​code>​ </​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.+The structure of this input section is quite straightforward ​and consists of three parts: 
 +  - Constraint definitions (type, which atoms to include, constraint target, etc) 
 +  - CDFT SCF loop settings (solver, convergence criterion, etc) 
 +  - Constraint weight function specific settings (Becke/​Hirshfeld subsections)  
 + 
 +In the above example, a Becke constraint is selected using the keyword [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT#​TYPE_OF_CONSTRAINT|TYPE_OF_CONSTRAINT]]. The actual constraints are defined using the section [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​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/​CDFT/​ATOM_GROUP#​ATOMS|ATOMS]] and the keyword [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​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/​CDFT#​TARGET|TARGET]] and [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT#​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/​CDFT/​ATOM_GROUP#​CONSTRAINT_TYPE|CONSTRAINT_TYPE]]. It is also possible to use fragment based constraints [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​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 2.  
 + 
 +{{ howto:​cdft-fragment-constraint.png?​350 }} 
 +**Figure 2.** Using a fragment based CDFT 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.  
 + 
 +The OUTER_SCF section within the CDFT section defines settings for the CDFT SCF loop. 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. 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 =====+==== Available ​constraints =====
  
-As explained in the previous section, the CDFT module in CP2K currently ​only supports using 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 atomsThe 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.+The CDFT module in CP2K currently supports using Becke or Hirshfeld (no forces) ​based constraints. The main aspects of these weight functions and their use as CDFT constraints will be explained in this section. Weight function specific settings are defined in the sections ​[[inp>​FORCE_EVAL/​DFT/​QS/CDFT/​BECKE_CONSTRAINT|BECKE_CONSTRAINT]] [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​HIRSHFELD_CONSTRAINT|HIRSHFELD_CONSTRAINT]]. If you want to visualize ​the weight functions ​to e.g. see the effects ​of using different parameters, you can use the section [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​PROGRAM_RUN_INFO/​WEIGHT_FUNCTION|WEIGHT_FUNCTION]] ​to print the CDFT weight ​function ​to a cube file
  
-{{ 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.+== Becke constraints ==
  
-The algorithmic implementation of the Becke density partitioning method ​has been detailed [[doi>​10.1021/​acs.jctc.6b01085|here]]. In briefthis 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 Becke density partitioning method ​can be considered as a smoothed Voronoi scheme. In Voronoi partitioningthe 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 3. 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.
  
-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 constraintUsually all coefficients are set to +1but mixing +1 and -1 coefficients would define ​the constraint as the difference between two groups of atomsThe 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-becke-atomicsize.png?​400}}{{ howto:​cdft-becke.png?400 }}\\ 
 +**Figure 3.** Comparison ​of the Voronoi (lines) ​and Becke partitioning (contours) schemesAt left, the Becke partitioning is performed without atomic size informationAt right, the size of the red atom is 30 % larger than the black atomsand the contours ​of the red atom extend farther than without atomic size adjustments
  
-{{ howto:​cdft-fragment-constraint.png?350 }} +The Voronoi and, by extension, the Becke partitioning methods treat each element equallyThis leads to unphysical partial charges in most systemsFor 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/​CDFT/​BECKE_CONSTRAINT#​ADJUST_SIZE|ADJUST_SIZE]] and the atomic radii are defined ​with the keyword [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​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 3 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. 
-**Figure 3.** Using a fragment based Becke constraintThe 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 calculationwhere the constraint target value is calculated from the superimposed fragment densities+ 
 +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/​CDFT/​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 functionsactivated by the keyword [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​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.
  
 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. 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> <code cp2k>
-&QS+&CDFT
   ...   ...
   &​BECKE_CONSTRAINT   &​BECKE_CONSTRAINT
Line 149: Line 191:
     ADJUST_SIZE ​    FALSE     ADJUST_SIZE ​    FALSE
     ATOMIC_RADII ​   0.63 0.32     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 scheme
     CUTOFF_TYPE ​    ​ELEMENT     CUTOFF_TYPE ​    ​ELEMENT
Line 166: Line 202:
     IN_MEMORY ​      TRUE     IN_MEMORY ​      TRUE
     SHOULD_SKIP ​    TRUE     SHOULD_SKIP ​    TRUE
-    ! Constraint definitions,​ each repetition defines a new constraint +&END CDFT
-    &​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>​ </​code>​
 +
 +== Hirshfeld constraints ==
 +
 +Hirshfeld constraints are cheaper to construct than Becke constraints in large systems because Hirshfeld constraints are essentially just weighted sums of spherical Gaussian functions. The keywords [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​HIRSHFELD_CONSTRAINT#​SHAPE_FUNCTION|SHAPE_FUNCTION]] and [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​HIRSHFELD_CONSTRAINT#​GAUSSIAN_SHAPE|GAUSSIAN_SHAPE]] define which type of Hirshfeld constraint to apply to the system. The shape function keyword accepts two values: Gaussian or Density. ​
 +  * The first choice implies that the CDFT weight function for each atom is a single Gaussian function whose radius is controlled by the GAUSSIAN_SHAPE keyword. By default, tabulated covalent radii are used as the radii of the Gaussian, but it is also possible to select van der Waals radii or to define custom radii. ​
 +  * The latter choice implies that the atomic weight function are constructed from isolated atomic densities which are expanded in terms of multiple spherical Gaussians. This choice is affected by the selected basis set and pseudopotential but not by the Gaussian shape keyword.  ​
  
 ==== Selected examples ​ ===== ==== Selected examples ​ =====
Line 293: Line 315:
   - 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.   - 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)   - 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)
 +
 +
 +=== Zn dimer cation with Hirshfeld constraints ===
 +
 +<note important>​This simulation requires CP2K version 7.0 or later.</​note>​
 +
 +This tutorial is exactly the same as the Zn dimer example above but using Hirshfeld partitioning based constraints instead of Becke constraints. You can find the input files {{:​howto:​cdft-tutorial-hirshfeld.zip|here}}. ​
 +
 +It might be instructive to visualize how the Becke and Hirshfeld weight function schemes differ, in particular, how the methods assign a volume to each atom in the system. You can activate the section [[inp>​FORCE_EVAL/​DFT/​QS/​CDFT/​PROGRAM_RUN_INFO/​WEIGHT_FUNCTION|WEIGHT_FUNCTION]] to output the weight function as a cube file which you can visualize with e.g. VMD. Feel free to modify the water tutorial above to look at the differences between Becke and Hirshfeld constraints in a system with different chemical elements.
  
  
Line 403: Line 434:
 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. ​ 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.+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 or a larger value if using more than 2 CDFT states. In this case, each CDFT state is solved ​in parallel using $N/​N_\mathrm{groups}$ processors. This will likely reduce the wall clock time of your simulation at the expense of more computing resources. However, you should note that the weight function and its gradients are computed separately for each state instead of copied from state to state (if possible). This can be costly for large solvated systems. 
 + 
 +A special run type is available for ''​NGROUPS 2''​ if the keyword [[inp>​FORCE_EVAL/​MIXED/​MIXED_CDFT#​PARALLEL_BUILD|PARALLEL_BUILD]] is activated. In this casethe CDFT weight function and gradients ​are first built in parallel on $N$ MPI processes, ​which are subsequently copied onto the two MPI processor groups of size $N/2$ which solve the CDFT states in parallel. 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 ====  ==== Example: Electronic coupling of Zn cation dimer ==== 
Line 450: Line 485:
  
 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]]. ​ 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]]. ​
 +
 +
 +==== Example: Configuration interaction calculations with CDFT (CDFT-CI): The case of $\mathrm{H}_2^+$ ==== 
 +
 +DFT exchange-correlation functionals suffer from varying degrees of self-interaction error (SIE). A pathological example of a system where SIE leads to unphysical results with DFT is the simple dissociation reaction
 +
 +$\mathrm{H}_2^+ \rightarrow \mathrm{H}^+ + \mathrm{H}$.
 +
 +Even though this system contains only 1 electron, the dissociation profile obtained with PBE notably deviates from the exact Hartree-Fock profile as shown in Figure 4 below.
 +
 +{{ howto:​cdftci-h2-dissociation.png?​450 }}
 +** Figure 4. ** Illustration of DFT self-interaction error for the reaction $\mathrm{H}_2^+ \rightarrow \mathrm{H}^+ + \mathrm{H}$. PBE notably deviates from the exact Hartree-Fock dissociation profile. The correct profile can be recovered with constrained DFT configuration interaction (CDFT-CI) using fragment constraint states $|\mathrm{H^+H}>​$ and $|\mathrm{HH^+}>​$ as the basis.
 +
 +We can use CDFT states as the basis of a configuration interaction (CI) simulation to correct for SIE in this system. As the figure above shows, CDFT-CI using the PBE functional is able to reproduce the exact dissociation profile. You can read up on the theory behind CDFT-CI simulations from the references given at the start of this tutorial. Very briefly, CDFT-CI simulations involve representing the system'​s wavefunction as a linear combination of multiple CDFT states where the charge/spin density is constrained differently in different states. The CI expansion coefficients and energies are then obtained by solving a generalized eigenvalue equation where the effective Hamiltonian matrix describes how the CDFT states interact with each other. ​
 +
 +In this tutorial, you will reproduce the DFT and CDFT results from the figure above. You can find the input files {{:​howto:​cdft-tutorial-h2.zip|here}}. The reference data used to plot Figure 4 are also included in the zip-folder. Please note that the reference results were obtained with a larger basis set and planewave cutoff as well as tighter convergence criteria than the settings you will be using in this tutorial.
 +
 +  - Start by examining the simulation script ''​energy.bash''​. This tutorial involves a rather large number of simulations so running them will take a while. You can use the flag ''​-x''​ to separately run the different types of simulations (DFT, CDFT, CDFT-CI) needed in this tutorial.
 +  - While the CDFT simulations are running, look at the results from the DFT simulations with PBE. Can you figure out the reason why PBE predicts an unphysical dissociation profile? (Hint. Compute the partial charges).
 +  - Inspect the output files produced by the CDFT-CI simulations once they are done. Find the output from the CDFT-CI module in the main output files. Look at the CI expansion coefficients in terms of the CDFT states. How would you characterize the CI wavefunction?​ Is the result sensible? What about the atomic partial charges in different CDFT states as a function of distance? ​
 +  - Plot the CDFT-CI and DFT dissociation profiles with your favorite plotting tool. Use the Hartree-Fock data from the provided data file as a reference. You can produce a similar data file from your simulations by calling ''​energy.bash''​ with the flag ''​-x results''​. ​
howto/cdft.1532410714.txt.gz · Last modified: 2018/07/24 07:38 by nholmber