Introduction

by Donald Bashford, St. Jude Children's Research Hospital Memphis, July 8, 2009.

SCRF (Self-Consistent Reaction Field) is a method of accounting for the effect of a polarizable solvent on the quantum system. The solvent is modeled as a dielectric continuum with a dielectric constant, EPSSOL, that fills the space outside the quantum system. The boundary between the interior (where the dielectric constant is unity) and the higher-dielectric exterior is the molecular surface, as defined by Connolly [251]. The charges of the quantum system cause polarization of this continuum, giving rise to a reaction field which then acts back on the quantum system potentially altering its charge distribution. The algorithm calculates the reaction field through solutions to the Poisson or Poisson-Boltzmann equation, and iteratively obtains self-consistency between the reaction field and charge distribution of the quantum system. These ideas have their roots in Onsager's consideration of a dipole and a dipole polarizability inside a sphere [252], and have also been developed in the Polarizable Continuum Model of Tomasi and co-workers [253, 254]. The first use of these ideas in DFT calculations using ADF was by Chen et al. [255] and Mouesca et al. [256]. At around the same time the Tomasi group also used the PCM model with DFT [257]. In the past, SCRF calculations with ADF were done in the research groups of Noodleman, Bashford and Case using custom modifications of ADF. Now the method is available in a standardized form.

The specific algorithm is:

1. Solve for the electronic structure in vacuum by the usual QM methods.
2. Derive a nucleus-centered partial charges from the electronic structure using the CHELPG algorithm [258].
3. Using the CHELPG-calculated charges: (a) Solve the Poisson equation (or, if an ionic strength is specified, the Poisson-Boltzmann equation) for φ_sol, the electrostatic potential in the presence of the above-described solvent dielectric environment. (b) Solve the Poisson equation for φ_vac, the electrostatic potential in a uniform vacuum environment. (c) The reaction field potential φ_rf is the difference between these two potentials: φ_rf = φ_sol - φ_vac.
4. Recalculate the electronic structure in the presence of the reaction field. This is done by adding φ_rf to the total potential evaluated on the numerical integration grid when calculating Fock matrix elements.
5. Check changes in energy for convergence. If not converged, return to step 2.

The algorithm also includes a correction for the small difference between the fit and true electron densities.

Several user-settable options can affect the SCRF procedure.

The choice of atomic radii and the probe radius determines the location of the dielectric boundary. This molecular surface comprises spherical patches of the contact surface generated by a solvent-sized probe rolling over the atomic radii, the toroidal surfaces swept out as the probe roles in grooves between pairs of atoms, and the inverse spherical patches generated when the probe simultaneously touches three or more atoms [251, 259]. In volumetric terms, the molecular surface divides the space of all points that are accessible to any part of a probe sphere that cannot penetrate into any of the atomic spheres, from the space that is not accessible [260]. The solvent accessible volume is assigned the solvent dielectric constant, while the inaccessible volume is assigned a dielectric constant of unity (the vacuum dielectric constant). Smaller radii move the dielectric surface closer to the atomic nuclei which typically leads to stronger calculated solvent effects.

The CHELPG routine for calculating atomic partial charges chooses charges that best reproduce the potential outside the molecule that is generated by the nuclei and the electron density. It sets up a grid of potential-sampling points in a region outside the molecule, calculates the potential on this grid due to the electron density and nuclei, and then finds a set of nucleus-centered charges that provides the best fit, in a least-squares sense, to the potential on the sampling points. The charge optimization is done using a singular value decomposition (SVD) method described by Mouesca et al. [256]. These calculations can be affected by user options concerning constraints on total charge and dipole, charge-fitting grid spacing and SVs to be deleted.

The solution of the Poisson or Poisson-Boltzmann equation utilizes libraries from Donald Bashford's MEAD programming suite. These use a finite-difference method that involves setting up cubic lattices around the molecule. Finer grids can be nested inside coarser ones to help manage trade-offs between accuracy and computational cost. The finest grid should cover the entire quantum system (that is, regions of significant electron density), and for good accuracy of φ_rf should be no coarser than about 0.15 Å. The outermost grid should extend 10 to 15 Å into the space beyond the model so that boundary conditions are accurate. A reasonable scheme for grid selection based on atomic coordinates is implemented as the default, but the MEAD grids are also user-adjustable.

The next section describes each input option in detail.