Theor Chem Acc (2016) 135:6
1 3
naphthalene example (for comparison, the calculation of
the nuclear Hessian or force constant matrix requires only
54 CPHF equations). While the effort is still manageable in
this case, the method of Ref. [ 29 ] has a steep formal scaling of O( N
2 ) solutions of the CPHF/CPKS equations, each
of which scales formally as O( N
4 ). This leads to an overall
formal scaling of O( N
6 ) and becomes impractical for large
molecules and basis sets. This method does not lend itself
to truncation because it is formulated in delocalized canonical MO basis.
The large computational effort associated with the
solution of the CPHF/CPKS equations can be avoided by
neglecting the change of the effective one-electron (Fock)
operator caused by the orbital perturbations. This is known
in the chemistry literature as the uncoupled Hartree–Fock
(UCHF) approximation, and in physics as the independent
particle approximation. With this approximation, the zerofrequency propagator is given, assuming real orbitals, as
Unlike the true propagator, the UCHF approximation is
given by a simple closed formula and requires only minimum computational effort to evaluate “on the fl y” if the
orbitals are available. The uncoupled Hartree–Fock/Kohn–
Sham approximation has almost completely vanished from
the chemistry literature about 40 years ago when modern
derivative techniques became available because of the poor
results it produced for second-order properties. Some systematic expositions of analytical derivative methods still
use it as a starting point, but it is in our opinion pedagogically inappropriate, as it requires considerable effort to
recover the coupled-perturbed Hartree–Fock results which
can be derived in a simpler way. UCHF/UCKS is still used
in some approximate theories, but we suspect that its only
merit is easy computability. According to Geerlings et al.
[ 29 ], the polarizabilities derived from the uncoupled density response function correlate well with accurate results
but can be off by up to a factor of 2, and thus they are only
qualitatively useful. Our results in Table 1 confi rm this.
2 Generalized polarizabilities for ultrafast
quantum/molecular mechanics simulations
Recently we developed a quantum/molecular mechanics
(QM/MM) method for the accurate statistical simulation of
molecules in polar environments such as in aqueous solutions [ 30 , 31 ], which speeds up the calculations by 4 or
more orders of magnitude or more, even for relatively small
solutes and at the inexpensive DFT level. The generalized
polarizabilities in our protocol can be used for the effi cient
(2)
χ
0
(r, r
) =
δρ(r)
δu(r )
= 4
i,a
φ i (r)φ i (r )φ a (r)φ a (r )
ε i − ε a
.
evaluation of the density–density response function. In this
section, we give a brief description of our ultrafast QM/
MM method to clarify the connection between the two very
different methods. However, only the generalized polarizabilities of the quantum system (and not the models for the
solvent or the statistical sampling) are relevant for the calculation of the DRF.
In our simulations, a small molecule or region, which we
will call the “solute,” is described by a quantum mechanical
method, and its much larger environment (the “solvent”) is
modeled by simpler and much faster molecular mechanics. The most important and computationally most expensive interaction between the two subsystems is long-range
electrostatics. Van der Waals interaction has a much shorter
range. In general, the environment is very fl exible and must
be statistically sampled for thermodynamic properties like
free energy. For reliable statistical predictions, particularly of entropic properties, 1–50 million solvent confi gurations have to be sampled [ 32 ]. This is straightforward if
the polarization of the quantum subsystem by the solvent is
neglected, as it reduces to the evaluation of the electrostatic
interaction between the solute and the solvent (usually represented by point charges), and does not require quantum
mechanical calculations for each solvent confi guration.
However, the polarizability of the solute is important, particularly in polar solvents. Statistical (for instance, Monte
Carlo) simulation becomes extremely time-consuming if
the polarization energy of the quantum system in the fi eld
of the solvent molecules is calculated by quantum mechanics, even if a single calculation takes only a fraction of a
minute.
We eliminate the need for separate QM calculations for
each solvent confi guration by approximating the electric
potential of the solvent within the solute molecule as linear
combination of predefi ned basis functions for the potential :
The potential within the volume of the solute molecule
is a smooth function and can be expanded in an appropriate set of expansion functions g k ( r) . This is a generalization
of the usual multipole expansion. The latter uses Cartesian
monomials 1, x , y , z , x
2 , xy , y
2 , xz , yz , z
2 , x
3 , …, or the corresponding solid spherical harmonics to expand the potential. However, for reasons explained in [ 17 ], the origin-centered multipole expansion is unsuitable for most systems.
We experimented with several expansion sets. One important requirement is that the expansion functions should not
diverge to infi nity like the solid harmonics do. We fi nally
settled on a sine function expansion of the potential. It is
conveniently defi ned in an outer box that is larger than the
extent of the molecular electronic density, to avoid problems with the periodic nature of the sine expansion. Other
(3)
U(r) ≈
k
c k g k (r) k = 1, . . . M
11
Reprinted from the journal
1 3
naphthalene example (for comparison, the calculation of
the nuclear Hessian or force constant matrix requires only
54 CPHF equations). While the effort is still manageable in
this case, the method of Ref. [ 29 ] has a steep formal scaling of O( N
2 ) solutions of the CPHF/CPKS equations, each
of which scales formally as O( N
4 ). This leads to an overall
formal scaling of O( N
6 ) and becomes impractical for large
molecules and basis sets. This method does not lend itself
to truncation because it is formulated in delocalized canonical MO basis.
The large computational effort associated with the
solution of the CPHF/CPKS equations can be avoided by
neglecting the change of the effective one-electron (Fock)
operator caused by the orbital perturbations. This is known
in the chemistry literature as the uncoupled Hartree–Fock
(UCHF) approximation, and in physics as the independent
particle approximation. With this approximation, the zerofrequency propagator is given, assuming real orbitals, as
Unlike the true propagator, the UCHF approximation is
given by a simple closed formula and requires only minimum computational effort to evaluate “on the fl y” if the
orbitals are available. The uncoupled Hartree–Fock/Kohn–
Sham approximation has almost completely vanished from
the chemistry literature about 40 years ago when modern
derivative techniques became available because of the poor
results it produced for second-order properties. Some systematic expositions of analytical derivative methods still
use it as a starting point, but it is in our opinion pedagogically inappropriate, as it requires considerable effort to
recover the coupled-perturbed Hartree–Fock results which
can be derived in a simpler way. UCHF/UCKS is still used
in some approximate theories, but we suspect that its only
merit is easy computability. According to Geerlings et al.
[ 29 ], the polarizabilities derived from the uncoupled density response function correlate well with accurate results
but can be off by up to a factor of 2, and thus they are only
qualitatively useful. Our results in Table 1 confi rm this.
2 Generalized polarizabilities for ultrafast
quantum/molecular mechanics simulations
Recently we developed a quantum/molecular mechanics
(QM/MM) method for the accurate statistical simulation of
molecules in polar environments such as in aqueous solutions [ 30 , 31 ], which speeds up the calculations by 4 or
more orders of magnitude or more, even for relatively small
solutes and at the inexpensive DFT level. The generalized
polarizabilities in our protocol can be used for the effi cient
(2)
χ
0
(r, r
) =
δρ(r)
δu(r )
= 4
i,a
φ i (r)φ i (r )φ a (r)φ a (r )
ε i − ε a
.
evaluation of the density–density response function. In this
section, we give a brief description of our ultrafast QM/
MM method to clarify the connection between the two very
different methods. However, only the generalized polarizabilities of the quantum system (and not the models for the
solvent or the statistical sampling) are relevant for the calculation of the DRF.
In our simulations, a small molecule or region, which we
will call the “solute,” is described by a quantum mechanical
method, and its much larger environment (the “solvent”) is
modeled by simpler and much faster molecular mechanics. The most important and computationally most expensive interaction between the two subsystems is long-range
electrostatics. Van der Waals interaction has a much shorter
range. In general, the environment is very fl exible and must
be statistically sampled for thermodynamic properties like
free energy. For reliable statistical predictions, particularly of entropic properties, 1–50 million solvent confi gurations have to be sampled [ 32 ]. This is straightforward if
the polarization of the quantum subsystem by the solvent is
neglected, as it reduces to the evaluation of the electrostatic
interaction between the solute and the solvent (usually represented by point charges), and does not require quantum
mechanical calculations for each solvent confi guration.
However, the polarizability of the solute is important, particularly in polar solvents. Statistical (for instance, Monte
Carlo) simulation becomes extremely time-consuming if
the polarization energy of the quantum system in the fi eld
of the solvent molecules is calculated by quantum mechanics, even if a single calculation takes only a fraction of a
minute.
We eliminate the need for separate QM calculations for
each solvent confi guration by approximating the electric
potential of the solvent within the solute molecule as linear
combination of predefi ned basis functions for the potential :
The potential within the volume of the solute molecule
is a smooth function and can be expanded in an appropriate set of expansion functions g k ( r) . This is a generalization
of the usual multipole expansion. The latter uses Cartesian
monomials 1, x , y , z , x
2 , xy , y
2 , xz , yz , z
2 , x
3 , …, or the corresponding solid spherical harmonics to expand the potential. However, for reasons explained in [ 17 ], the origin-centered multipole expansion is unsuitable for most systems.
We experimented with several expansion sets. One important requirement is that the expansion functions should not
diverge to infi nity like the solid harmonics do. We fi nally
settled on a sine function expansion of the potential. It is
conveniently defi ned in an outer box that is larger than the
extent of the molecular electronic density, to avoid problems with the periodic nature of the sine expansion. Other
(3)
U(r) ≈
k
c k g k (r) k = 1, . . . M
11
Reprinted from the journal