Brill Academic Publishers
P.O. Box 9000, 2300 PA Leiden,
The Netherlands

Lecture Series on Computerand Computational Sciences
Volume 4, 2005, pp. 663-666

Charge Response Kernel Theory based on Ab Initio andDensity Functional Calculations
Akihiro Morita1 and Tateki Ishida
Department of Computational Molecular Science
Institute for Molecular ScienceMyodaiji, Okazaki 444-8585, Japan
Received 9 July, 2005; accepted in revised form 14 August, 2005
Abstract: The charge response kernel (CRK) of a molecule, ∂Qa /∂Vb , is formulated andcalculated via ab initio molecular orbital and density functional theory, where Qa is thepartial charge of the site a and Vb is the electrostatic potential at the site b. The nonempirical definition of CRK provides a general and rigorous scheme of polarizable molecularmodeling.
Keywords: charge response kernel, polarizability, density functional theory

1Introduction

The electronic polarization is of vital importance in intermolecular interactions in condensed
phases, in particular, to correctly describe solvation energetics and structures. Therefore, most
molecular models in common use incorporate the electronic polarization effects in either implicit orexplicit way. The former, implicit way employs a nonpolarizable molecular model with augmented
polarization, such as the SPC [1] or TIP4P [2] water model, while the latter uses a polarizable
molecular model. Although the latter is computationally more demanding than the former foruse of molecular simulation, it is capable of describing both spatial and temporal fluctuation ofmolecular polarization, which often plays crucial roles in local solvation or interface structures[3, 4].
The charge response kernel (CRK) is considered as a general model to explicitly represent
the molecular polarization [5]. It is defined as (∂Qa /∂Vb )N on the basis of the interaction site
representation, where Qa is the partial charge at the site a and Vb is the electrostatic potential
at the site b. This property is readily incorporated into the molecular dynamics simulation basedon the interaction site molecular models, introducing fluctuating partial charges on the interaction
sites.
We propose the ab initio formulation of the CRK, due to the following advantages. First, theab initio theory obviates empirical parameterization to determine the CRK, and therefore it isgenerally applicable to any molecule by ab initio calculations. Second, the accuracy and reliabilityof the CRK are warranted by those of the underlying electronic structure theory at any level ofaccuracy. Third, these general features of the ab initio CRK provide an equal footing to treat theelectronic polarization of solute and solvent molecules, including novel solute species that are often
difficult to deal with by an empirical model. As one example, the CRK model in combination
1 Corresponding author. E-mail: amorita@ims.ac.jp

664

A. Morita, T. Ishidaof molecular dynamics simulation was successfully utilized to elucidate the anomalous diffusion ofsome aromatic transient radicals in solutions [6].
In the present paper, first we briefly summarize the theoretical framework of the CRK modelin the ab initio molecular orbital theory. Then we extend the formulation to the density functional
theory, and demonstrate the accuracy of the calculated results.

2Ab Initio Theory

In our theoretical treatment, the partial charge Qa is defined as a electrostatic potential charge,
which is optimized to reproduce the surrounding electric field generated by the molecule. The
partial charge Qa is given by an one-electron operator Q̂a as Qa = hΨ|Q̂a |Ψi, where Ψ denotes thewavefunction of the many-electron system. The detailed form of Q̂a is presented in our previousworks [5].
The molecular Hamiltonian in the condensed phase, Ĥ, is given by Ĥ = Ĥ0 + Ĥ 0 , where Ĥ0is the molecular Hamiltonian in the gas-phase and Ĥ 0 is the perturbation Hamiltonian Ĥ 0 by thesolvation environment. The latter, Ĥ 0 , is given in the interaction site representation asĤ 0 =sites
XQ̂a Va ,(1)awhere Va is the electrostatic potential at the site a generated by the ambient solvent.
The CRK Kab = (∂Qa /∂Vb ), by definition, is represented via the first derivative of the wavefunction Ψb with respect to the external electrostatic potential Vb at site b,Kab ≡∂2E
∂Qa== hΨb |Q̂a |Ψi + hΨ|Q̂a |Ψb i.
∂Vb∂Va ∂Vb(2)Note that the CRK is equivalent to the second derivative of the total energy E. The derivative ofthe wavefunction is generally formulated by the coupled-perturbed Hartree-Fock equation in theab initio theory [7]. We first illustrate this scheme in a simple case of the closed-shell Hartree-Focktheory.
aSuppose the molecular orbital and its derivative by Va are denoted by Cpi and Cpi. In thefollowings, suffixes p, q, ... refer to atomic orbital (AO), and i, j, ... to molecular orbital (MO).
aThe derivative coefficient Cpi
is expressed in the MO representation using the transformation
matrix U ,MOXaaCpj Uji,(3)Cpi=jand U is formulated by the following coupled-perturbed Hartree-Fock equation,(l − i )Ulia +occ XvirXj(M O)(M O)

aHlikj Ukj= −Qa,li ,(4)kwhere j and k indicate the occupied and virtual MO’s, respectively. i is the orbital energy of thei-th canonical MO, and H is given with the two-electron integrals (in the MO representation) as
Hlikj = 4(li|kj) − (lk|ij) − (lj|ik).

(5)

Eq. (4) is a set of linear equations for U , which derives the CRK as follows,
Kab =occ XvirXij(M O)

b4Qa,ij Uji.

(6)

665

Charge Response Kernel

3Density Functional Version

The density function theory (DFT) is in widespread use in the last decade in the field of quantum
chemistry [8]. Its amazing success hinges on the exchange-correlation functionals developed todate, such as B3LYP [9, 10], which can take account of the electron correlation effects fairly
accurately at a modest computational cost. Since the Kohn-Sham formalism is quite analogous to
the Hartree-Fock formalism, the extension of the CRK model to the density functional theory is
straightforward, as shown in the followings.
The Kohn-Sham matrix is represented asXKSxcFpq= hpq +(pq|rs)Drs + Vpq[ρ],(7)r,sxc
where Vpq[ρ] = δE xc [ρ]/δDpq is the exchange-correlation potential, which is the derivative of theexchange-correlation energy functional E xc [ρ] with respect to the density matrix Dpq .

The formulation of the Kohn-Sham matrix, (7), coincides with that of P
the Fock matrix byxc
replacing the exchange-correlation potential Vpq
with the exchange term, − r,s (pr|qs)Drs . Accordingly, eq. (4) is modified in the DFT case by replacing H with the Kohn-Sham equivalentas,xc
Hpqrs → 4(pq|rs) + δVpq/δDrs .
(8)
The DFT version of eq. (4) is often called coupled-perturbed Kohn-Sham equation. Details for thesecond term are given in the general analytical second derivative of the density functional theory[11].
We calculated the CRK via the DFT and compared them to those by the Hartree-Fock case.
Table 3 illustrates the results of NH3 . The program code for the CRK calculation was implemented
into the GAMESS-UK package [12]. To evaluate the reliability of the calculated CRK, the electronic
polarizability derived from the CRK is compared to the experimental value. The polarizability α ij
is derived from the CRK as
αij = −sites
Xa,b

∂Qa
ri (a) rj (b),∂Vbi, j = x, y, z,(9)where ri (a) is the i-th Cartesian coordinate of the site a. The calculated CRK derives an isotropic
polarizability of ᾱ = 12.13 a.u. (HF) and 13.59 a.u. (B3LYP), while the experimental value is 15.25
a.u. While both results are fairly consistent to the experiment, the density functional calculationshows some improvement in evaluating the electronic polarization. Computational cost of the CRK
is almost equivalent to that of the analytical calculation of the conventional dipole polarizability.
Table 1: Charge Response Kernel of NH3 molecule, calculated by (a) HF and (b) B3LYP. The
basis set is d-aug-cc-pVTZ for N and aug-cc-pVTZ for H, with s and p diffusion functions only,(12s7p2d1f/6s3d1d) / [6s5p2d1f/4s3p1d]. Unit: a.u.

1(N)
2(H)
3(H)
4(H)

1-23.98
7.99
7.99
7.99

(a) HF
2-4.38-1.81-1.81

3-4.38-1.81

4-4.38

1(N)
2(H)
3(H)
4(H)

(b) B3LYP
123-28.39
9.46 -5.03
9.46 -2.22 -5.03
9.46 -2.22 -2.22

4

-5.03

666

A. Morita, T. Ishida

4SummaryIn this work, the Charge Response Kernel (∂Qa /∂Vb ) formulated via the ab initio molecular orbital
theory by us is extended to the density functional theory. These definitions of CRK rigorously
based on the quantum chemical theories can provide a general scheme to construct polarizable
molecular models. We plan to adopt this CRK model to molecular simulation of the interfacial
nonlinear spectroscopy. Such work is now in progress in our group.

Acknowledgment
This work was supported by the Grant-in-Aid (no. 15550012) and National Research Grid Initiative(NAREGI) of the Ministry of Education, Japan.

References[1] H. J. C. Berendsen, J. P. M. Postma, M. F. van Gunsteren, and J. Hermans. IntermolecularForces, p. 331. Reidel, Dordrecht, 1981.
[2] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. J. Chem.
Phys., 79, 926 (1983).
[3] J. Caldwell, L. X. Dang, and P. A. Kollman. J. Am. Chem. Soc., 112, 9144 (1990).
[4] P. Jungwirth and D. J. Tobias. J. Phys. Chem. B, 106, 6361 (2002).
[5] A. Morita and S. Kato. J. Am. Chem. Soc., 119, 4021 (1997).
[6] A. Morita and S. Kato. J. Chem. Phys., 108, 6809 (1998).
[7] P. Pulay. Adv. Chem. Phys., 69, 241 (1987).
[8] R. G. Parr and W. Yang. Density-Functional Theory of Atoms and Molecules. Oxford Univ.,
New York, 1989.
[9] A. D. Becke. J. Chem. Phys., 98, 5468 (1993).
[10] C. Lee, W. Yang, and R. G. Parr. Phys. Rev. B, 37, 785 (1988).
[11] N. C. Handy, D. J. Tozer, G. J. Laming, C. W. Murray, and R. D. Amos. Israel J. Chem.,
33, 331 (1993).
[12] GAMESS-UK is a package of ab initio programs written by M.F. Guest, J.H. van Lenthe, J.
Kendrick, K. Schoffel, and P. Sherwood, with contributions from R.D. Amos, R.J. Buenker,H.J.J. van Dam, M. Dupuis, N.C. Handy, I.H. Hillier, P.J. Knowles, V. Bonacic-Koutecky, W.
von Niessen, R.J. Harrison, A.P. Rendell, V.R. Saunders, A.J. Stone, D.J. Tozer, and A.H. deVries. The package is derived from the original GAMESS code due to M. Dupuis, D. Spangler
and J. Wendoloski, NRCC Software Catalog, Vol. 1, Program No. QG01 (GAMESS), 1980.