Applications of Constrained Density Functional Theory (CDFT) by cp2k

We here explain how to apply constrained density functional theory (CDFT) to charge-transfer reactions using cp2k program package. We here focus on the proton transfer (PT) and electron/hole (polaron) formation. The tutorial on the cp2k website (https://www.cp2k.org/howto:cdft) offers the basic theory and information.

In protic solvents the transferring H atom moves between water molecules and defining one H atom is impossible. Thus, an OH distance is a poor reaction coordinate. On the other hand, in aprotic solvents, an OH distance should be a good reaction coordinate because the transferring proton stays between the reactant and product. Proton transfer in HF acid in water[2,3] is an exceptional example in that the product state is the contact-ion pair and the HF distance is a good reaction coordinate.

A good candidate which has been employed for a long time is the vertical energy gap coordinate between the reactant and product states. This has been first employed for electron transfer, and later extended to proton transfer.[4-6] To this end, CDFT implemented in cp2k provides a straightforward computation of energetics and dynamics. We here explain how to perform such calculation efficiently.

Gas-Phase Benchmark Calculation

Before moving on to full-scale computations in condensed phases, it is advisable to check the computational details using a small cluster because you have to choose the target value in CDFT. In the case of charge/magnetization constraint, the target is constrained number of valence electrons/spin density on a specified atom. We here employ the charge constraint in the following.

The following is the xyz coordinate file of a silicic-water cluster.

Starting from such a structure, we optimize the geometry for each mixing parameter LAMBDA to obtain the potential energy surface.

In the case of aqueous silicic acid, the target values of 5.9 (reactant) and 6.6 (product) by the following input

efficiently sample the proton transfer as shown in Figure 2 of ref [1].

Liquid-Phase Calculation

After determining the computational setup in the gas phase, we can move on to the liquid-phase calculation. The input file for MD simulation can be made by changing the "runtype" from "geo_opt" to "md" and adding the computational conditions of MD simulation. CDFT MD simulation produces the trajectory of diabatic energies in {Project name}-mix-1.ener. Putting the results into the free-energy perturbation formula such as in ref [1], the free-energy profile can be computed.

Charge Localization for Hole/Electron Polaron Formation

In photocatalysis, an irradiated light excites electrons in the valence band to the conduction band, leaving holes in the valence band. The CDFT method can also localize the photoexcited holes and/or electrons (polaron). The example is the hole formation of anatase TiO2 in ref [7].

The CDFT method provides a workaround for charge delocalization due to the effects of self-interaction errors in exchange-correlation functionals with the generalized gradient approximation (GGA). Thus, computationally-demanding hybrid functionals may not be necessary and the GGA functionals which are computationally less demanding will suffice.

Another advantage is straightforwardness. The CDFT method can specify atoms on which the charges and/or magnetizations are constrained. Therefore, any specification of arbitrary geometric coordinates such as bond lengths[8] and atomic replacement (https://wiki.wpi.edu/deskinsgroup/Localizing_Charge_-_VASP) are unnecessary.

In practice, the computational procedure for geometry optimization is as follows:

1. DFT calculation for the positively (+1)/negatively (-1) charged system.
   This step relaxes the surface structure and the positive charge is delocalized over the entire slab.

2. CDFT calculation for charge localization.
   CDFT optimization restarting from the structure and wavefunction optimized in the first step.

3. DFT calculation without constraint
   Confirm the charge is localized even without constraints.

Here we demonstrate the trapping of a photogenerated hole on hydroxyl groups in the hydroxylated anatase TiO2(112) surface which was found to have high photocatalytic activity.[7] Panel (a) of the upper figure shows the slab model of hydroxylated anatase TiO2(112) surface. The following is the xyz file of hydroxylated anatase TiO2(112). One dissociated water molecule was added to the surface of anatase TiO2(112).

The following is the excerpted CDFT section of the cp2k input file.

In the CDFT section, the keyword TYPE_OF_CONSTRAINT specifies the type of constraint to be used. In the example above, the Becke constraint is selected. The keyword TARGET defines the constraint target value. The constraint target value is the number of valence electrons on the constraint atom. For example, the number of valence electrons on an oxygen ion is 6, but in the above example, it is set to 5.4 so that the O atom is positively charged and a hole is localized. The keyword ATOMS in the ATOM_GROUP section specifies a constraint atom. In the example above, the 31th atom (the oxygen atom of the hydroxyl group, the line with a * in the xyz file) is the constraint atom. If you have atoms that you want to calculate their charges without applying any constraints, specify them in the keyword ATOMS in the DUMMY_ATOMS section. The keyword CONSTRAINT_TYPE specifies the type of constraint to be applied. In the example above, the total charge density is constrained (CHARGE).

Spin Density and Density of States

If you want to output the spin density of each atom, add the following sections to the input file. In the following example, Hirshfeld and Löwdin spin population analysis are specified.

If you want to visualize the spin density, you can generate a cube file by adding the following section to the input file.

By reading the cube file with a visualization program such as VESTA, the spin density can be visualized. Panel (b) shows a visualization of the spin density (hole) of OH groups on the anatase TiO2 surface.

If you want to calculate the projected density of states (PDOS), you can generate a pdos file by adding the following section to the input file.

The pdos file is output for each atom and has a name like md-k1-1.pdos. The pdos file outputs the results for each atomic orbital. For the convoluted DOS, it is better to use the Python script for convolution provided on the website. [https://wiki.wpi.edu/deskinsgroup/Density_of_States_-_CP2K] Download the two files pdos.py and get-smearing-pdos.py and, on the folder with the .py files, execute the Python script as follows

Panel (c) shows the PDOS of holes trapped at the surface OH of the hydroxylated (112) surface. The defect levels are generated in the band gap and a hole is trapped within the deep levels.

References

[1] "Constrained Density Functional Theory Molecular Dynamics Simulation of Deprotonation in Aqueous Silicic Acid", Tatsuya Joutsuka and Koji Ando, J. Phys. Chem. B, 124, 38, 8323–8330 (2020).

[2] “Hydration Structure in Dilute Hydrofluoric Acid”, Tatsuya Joutsuka and Koji Ando, J. Phys. Chem. A, 115 (5), 671–677 (2011).

[3] “Dynamics of Proton Transfer and Vibrational Relaxation in Dilute Hydrofluoric Acid”, Tatsuya Joutsuka and Koji Ando, J. Phys. Chem. A, 115 (5), 678–684 (2011).

[4] "Molecular Mechanism of HCl Acid Ionization in Water: ab initio Potential Energy Surfaces and Monte Carlo Simulations", Koji Ando and James T. Hynes, J. Phys. Chem. B, 101, 49, 10464–10478 (1997).

[5] “Molecular Mechanism of HF Acid Ionization in Water: An Electronic Structure-Monte Carlo Study”, Koji Ando and James T. Hynes, J. Phys. Chem. A, 103 (49), 10398–10408 (1999).

[6] "Using the Constrained DFT Approach in Generating Diabatic Surfaces and Off Diagonal Empirical Valence Bond Terms for Modeling Reactions in Condensed Phases", Gongyi Hong, Edina Rosta, and Arieh Warshel, J. Phys. Chem. B, 110 (39), 19570–19574, (2006).

[7] "Facet Dependence of Photocatalytic Activity in Anatase TiO2: Combined Experimental and DFT Study", Tatsuya Joutsuka, Hiroto Yoshinari, and Satoshi Yamauchi, Bull. Chem. Soc. Jpn., in press.

[8] "Does Polaronic Self-Trapping Occur at Anatase TiO2 Surfaces?", John J. Carey and Keith P. McKenna. J. Phys. Chem. C, 122 (48), 27540–27553, (2018).