Theor Chem Acc (2016) 135:6
1 3
choices such as distributed (atom-centered) Gaussiandamped spherical harmonics would also work. The latter
provide an alternative defi nition of distributed polarizabilities [ 3 ]. A plane wave expansion of the electric potential
has been used successfully for atoms a long time ago by
Koide [ 33 ].
In the ultrafast QM/MM method, we precalculate generalized polarizabilities , i.e., second derivatives with respect
to electrostatic perturbations representing the potential.
For single-determinant wavefunctions, this requires the
solution of the coupled-perturbed Hartree–Fock (or Kohn–
Sham) (CPHF/CPKS) equations for each perturbation g k ,
i.e., for each function used to expand the potential. The
solution is simplifi ed by the fact that the perturbation is a
one-electron operator. The CPHF procedure [ 27 , 34 ] is routinely included in most quantum chemistry programs; for
a review, see [ 35 ]. As the potential is irregular within the
molecule, a fairly large number of expansion (in our case,
sine) functions (several hundred to a few thousand for a
drug-sized molecule) are required to represent the electrostatic potential accurately. The result of the second derivative calculation in the fi eld-free ( λ = 0) case is a matrix of
generalized polarizabilities:
This is analogous to the usual dipole polarizability tensor which has only 3 × 3 components, corresponding to the
perturbations g ( r) = x , y , z. The second-order part of the
energy in the electric fi eld is given as
We calculate the electric potential within the molecular box on a grid, and these values are fi tted by a linear
combination of the potential basis functions, Eq. ( 3 ), by
a weighted least-squares procedure. This gives the coeffi -
cients of the expansion functions, c k in Eq. ( 3 ), as
Here, c is the (column) vector of the expansion
coeffi cients, A is a rectangular matrix, A μk = g k ( r μ ),
μ = 1,… N , i.e., the value of the expansion (in our case,
sine) function g k at grid point r μ . To ensure a unique
solution, the number of grid points N should exceed the
number of expansion functions M . W is a weight matrix,
see [ 31 ] for a discussion about appropriate weights. For
the calculation of DRFs, no weights were used. The vector u contains the values of the electrostatic potential at
the grid points within the molecule, u μ = U (r μ ), and L is
the N × M matrix ( A
T WA )
−1 A
T W . Calculating the coeffi cients c and the second-order energy, Eq. ( 5 ) involves
(4)
α kl = −
∂ 2 E
∂c k ∂c l
=0
(5)
E
(2)
= −
1
2
kl
c k c l α kl = −
1
2
c
T
αc
(6)
c = [(A
T WA)
−1 A
T W]u = Lu
only linear algebra manipulations and is very fast on modern CPUs.
3 Calculation of the density response function
from generalized polarizabilities
We will discuss in this section the determination of the
density response function from generalized polarizabilities. This has lower formal scaling, O( N
5 ), than the method
of Yang et al. [ 28 ], and its spatial resolution can be easily
adjusted. One possibility is taking the 3-D Fourier (actually,
sine) transform of the generalized polarizability expressed
by sine waves. We have implemented an alternative, more
general method. Substituting Eq. ( 6 ) into Eq. ( 5 ) gives
where L is defi ned in Eq. ( 6 ) and in the text. Recall that u
contains the values of the electrostatic potential at the grid
points and α is the matrix of generalized polarizabilities
(second derivatives with respect to modulated perturbations in the potential). Equation ( 7 ) gives the second-order
(polarization) energy as a discretized sum over the values
of the potential on the grid points. By considering this discrete sum as a quadrature formula for the two-dimensional
integral
where u ( r ) is the (continuous) electric potential, we can
identify the density response function at the grid points as
where the Δ v ’s are the weights, i.e., grid cell volumes corresponding to the points μ and ν . We use a regular grid
(with grid points where the electron density is practically
zero removed), with all cells of the same volume for our
sine function expansion, and thus the density response
function at the pair of points ( r μ , r ν ) is simply (Δ v )
−2 times
the μν matrix element in Eq. ( 9 ).
The computational cost of the transformation in Eq. ( 9 )
is negligible compared to the cost of calculating the generalized polarizabilities, α . Compared to the method used by
Geerlings et al. [ 29 ], the advantage of our method is that its
resolution, i.e., the size of the sine basis, can be defi ned by
the user and is not given by the atomic orbital basis set. The
main features of the DRF do not require high resolution.
The number of the expansion (sine) basis is proportional
to the molecular volume and at a constant resolution scales
linearly with molecular size. This quantity determines the
main computational cost, the calculation of the fi rst-order
perturbed wavefunctions by the CPHF or CPKS procedure
(7)
E
(2)
= −
1
2
u
T
(L
T
αL)u
(8)
E
(2)
= −
1
2
dr
3
dr
3 u(r)χ(r, r
)u(r
)
(9)
χ(r μ , r ν ) = −[L
T
αL
T
] μν ((v μ )
−1
((v ν )
−1
12
Reprinted from the journal
1 3
choices such as distributed (atom-centered) Gaussiandamped spherical harmonics would also work. The latter
provide an alternative defi nition of distributed polarizabilities [ 3 ]. A plane wave expansion of the electric potential
has been used successfully for atoms a long time ago by
Koide [ 33 ].
In the ultrafast QM/MM method, we precalculate generalized polarizabilities , i.e., second derivatives with respect
to electrostatic perturbations representing the potential.
For single-determinant wavefunctions, this requires the
solution of the coupled-perturbed Hartree–Fock (or Kohn–
Sham) (CPHF/CPKS) equations for each perturbation g k ,
i.e., for each function used to expand the potential. The
solution is simplifi ed by the fact that the perturbation is a
one-electron operator. The CPHF procedure [ 27 , 34 ] is routinely included in most quantum chemistry programs; for
a review, see [ 35 ]. As the potential is irregular within the
molecule, a fairly large number of expansion (in our case,
sine) functions (several hundred to a few thousand for a
drug-sized molecule) are required to represent the electrostatic potential accurately. The result of the second derivative calculation in the fi eld-free ( λ = 0) case is a matrix of
generalized polarizabilities:
This is analogous to the usual dipole polarizability tensor which has only 3 × 3 components, corresponding to the
perturbations g ( r) = x , y , z. The second-order part of the
energy in the electric fi eld is given as
We calculate the electric potential within the molecular box on a grid, and these values are fi tted by a linear
combination of the potential basis functions, Eq. ( 3 ), by
a weighted least-squares procedure. This gives the coeffi -
cients of the expansion functions, c k in Eq. ( 3 ), as
Here, c is the (column) vector of the expansion
coeffi cients, A is a rectangular matrix, A μk = g k ( r μ ),
μ = 1,… N , i.e., the value of the expansion (in our case,
sine) function g k at grid point r μ . To ensure a unique
solution, the number of grid points N should exceed the
number of expansion functions M . W is a weight matrix,
see [ 31 ] for a discussion about appropriate weights. For
the calculation of DRFs, no weights were used. The vector u contains the values of the electrostatic potential at
the grid points within the molecule, u μ = U (r μ ), and L is
the N × M matrix ( A
T WA )
−1 A
T W . Calculating the coeffi cients c and the second-order energy, Eq. ( 5 ) involves
(4)
α kl = −
∂ 2 E
∂c k ∂c l
=0
(5)
E
(2)
= −
1
2
kl
c k c l α kl = −
1
2
c
T
αc
(6)
c = [(A
T WA)
−1 A
T W]u = Lu
only linear algebra manipulations and is very fast on modern CPUs.
3 Calculation of the density response function
from generalized polarizabilities
We will discuss in this section the determination of the
density response function from generalized polarizabilities. This has lower formal scaling, O( N
5 ), than the method
of Yang et al. [ 28 ], and its spatial resolution can be easily
adjusted. One possibility is taking the 3-D Fourier (actually,
sine) transform of the generalized polarizability expressed
by sine waves. We have implemented an alternative, more
general method. Substituting Eq. ( 6 ) into Eq. ( 5 ) gives
where L is defi ned in Eq. ( 6 ) and in the text. Recall that u
contains the values of the electrostatic potential at the grid
points and α is the matrix of generalized polarizabilities
(second derivatives with respect to modulated perturbations in the potential). Equation ( 7 ) gives the second-order
(polarization) energy as a discretized sum over the values
of the potential on the grid points. By considering this discrete sum as a quadrature formula for the two-dimensional
integral
where u ( r ) is the (continuous) electric potential, we can
identify the density response function at the grid points as
where the Δ v ’s are the weights, i.e., grid cell volumes corresponding to the points μ and ν . We use a regular grid
(with grid points where the electron density is practically
zero removed), with all cells of the same volume for our
sine function expansion, and thus the density response
function at the pair of points ( r μ , r ν ) is simply (Δ v )
−2 times
the μν matrix element in Eq. ( 9 ).
The computational cost of the transformation in Eq. ( 9 )
is negligible compared to the cost of calculating the generalized polarizabilities, α . Compared to the method used by
Geerlings et al. [ 29 ], the advantage of our method is that its
resolution, i.e., the size of the sine basis, can be defi ned by
the user and is not given by the atomic orbital basis set. The
main features of the DRF do not require high resolution.
The number of the expansion (sine) basis is proportional
to the molecular volume and at a constant resolution scales
linearly with molecular size. This quantity determines the
main computational cost, the calculation of the fi rst-order
perturbed wavefunctions by the CPHF or CPKS procedure
(7)
E
(2)
= −
1
2
u
T
(L
T
αL)u
(8)
E
(2)
= −
1
2
dr
3
dr
3 u(r)χ(r, r
)u(r
)
(9)
χ(r μ , r ν ) = −[L
T
αL
T
] μν ((v μ )
−1
((v ν )
−1
12
Reprinted from the journal