VSCRF: Vertical Excitation Self-Consistent Reaction Field¶

Introduction

by Donald Bashford, St. Jude Children’s Research Hospital Memphis, October 26, 2011.

VSCRF is a method of accounting for the effect of the solvent environment (and possibly a macromolecular environment) on the energy of absorption/emission transitions between electronic states. It is built on the Delta SCF method of calculating electronic transition energies and a continuum dielectric model of the environment in a way that is consistent with the Franck–Condon principle. The original development and application of the procedure can be found in the papers of Liu et al. [1] and Han et al. [2].

Suppose one is calculating the energy for a transition between an initial state i with total density (including nuclei) \(\rho\)_i and a final state f with total density \(\rho\)_f . State i has existed long enough that the dielectric environment surrounding it has had time to fully respond, including changes in the distributions of nuclear coordinates and solvent-molecule dipole moments. The density \(\rho\)_i has therefore given rise to a reaction field \(\phi\)^(r) _i,eq and is self-consistent with it according to the ideas described in the SCRF part of the documentation, therefore the energy in solvent of the initial system is

where E₀ is the vacuum density functional. During the excitation process, the change in density \(\Delta\) \(\rho\)_fi = \(\rho\)_f - \(\rho\)_i must work against the non-optical component of the pre-existing equilibrium reaction field as well as the optical reaction field arising in response to \(\Delta\) \(\rho\) itself. The full theoretical development is given in the paper of Liu et al. [1] and the result is that the excitation energy is:

where \(\Delta\) \(\phi\)^(r) _op is the reaction field that arises from the optical dielectric constant which is related to the index of refraction: \(\epsilon\)_op = n² \(\approx\) 2 for organic liquids.

In the following description of the VSCRF algorithm, we adopt the language of an optical absorption calculation. The initial state is therefore referred to as the ground state, and the final state as the excited state. For emission calculations, these roles would be reversed. If there is a macromolecular environment surrounding the protein that modeled by macroscopic dielectric theory as described in the SCRF documentation in this manual, this is referred to as the ‘protein’. With these definitions in mind, the VSCRF algorithm is:

• An SCRF calculation for the ground state, either with or without a surrounding protein, is a prerequisite for a VSCRF calculation.
• The electronic structure of the excited state (with occupations initially selected by the user) is calculated in vacuum.
• The ground-state reaction field, the protein field (if any) and electron density are read in from the results of the previous ground-state calculation. It is essential that the excited-state and ground state calculations use exactly the same geometry, ADF grids, basis functions, XC functionals, MEAD grids and so on.
• The excited-state electron density calculated in step 1 is subtracted from the stored ground-state density to obtain \(\Delta\)\(\rho\).
• Nucleus-centered partial charges are derived from \(\delta\)\(\rho\) using a least-squares potential-fitting SVD algorithm [3] which is somewhat similar to the CHELPG algorithm [4].
• Set the exterior dielectric environment to the optical dielectric constant of the solvent, and the optical dielectric constant of the protein. These are equal to the square of the index of refraction of these materials, and are typically \(\approx\) 2 (the default value).
• Using the nucleus-centered partial charges for \(\Delta\)\(\rho\), solve the Poisson equation using the finite-difference method in the above dielectric environment to obtain \(\Delta\)\(\phi\)_sol .
• Set the exterior dielectric constant to 1 everywhere and solve the Poisson equation by the finite-difference method to obtain \(\Delta\)\(\phi\)_vac
• The optical reaction field to the excitation, \(\Delta\)\(\phi\)^(r) _op is the difference between \(\Delta\)\(\phi\)_sol and \(\Delta\)\(\phi\)_vac . In this subtraction, artifacts related to the finite difference method cancel out.
• Recalculate the electronic structure in the presence of the optical reaction field, as well as the previously stored ground-state reaction field and the protein field (if any). This is done by placing these potentials on the ADF numerical integration grid when calculating Fock matrix elements.
• Check changes in energy for convergence, and if not converged, return to step 6.

It should be noted that some additional calculations may be needed to obtain a proper excited-state energy. For example, if the first excited singlet state is prepared by promoting one of the \(\alpha\) electrons from the HOMO to the LUMO, the resulting (single-determinant) wave function is not a pure singlet state but actually a mixture of singlet and triplet states [5]. Calling this mixed state S’₁ , the resolution of the problem is to also calculate the energy of a corresponding triplet state, T₁ (prepared, for example, by adding one \(\alpha\) electron to the LUMO and removing one \(\beta\) electron from the HOMO) and using the formula [5]:

See Liu et al. [1] or Han et al. [2]. for specific examples of this in the ADF/VSCRF context.

VSCRF Input

The VSCRF input is contained in a VSCRF input block as shown below, optional keywords being surrounded by curly brackets.

VSCRF
  MEADGRID string integer real
  RADIUS string real
  INITIAL_STATE string
  {CYCLES integer}
  {TOLERANCE real}
  {ATOM_MAXR real}
  {CHGFIT_CONSTRAIN string}
  {DELATOM N Iat1 Iat2 ... IatN}
  {GRID_SPACING real}
  {SVD_CONDITION real}
  {SV_DELETE integer}
  {EPS_OPT_SOL real}
  {SOLRAD real}
  {PROTEIN string}
  {EPS_OPT_PROT real}
END

The VSCRF block contains three mandatory keys: MEADGRID, RADIUS and INITIAL_STATE. All other keys are optional. It is highly recommended that the parameters specifying the MEAD grid and the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations. These procedures involve some numerical error, and keeping parameters the same will promote cancellation of these errors as energy differences are taken as the final results.

MEADGRID string integer real

Specifies the centering style, dimension and spacing for the MEAD grid. Recognized centering styles are “ON_ORIGIN” and “ON_GEOM_CENT”. The grid dimension specifies the number of points on one edge of a cubic grid. The grid spacing is given in Angstroms. The edge length of the grid is the product of the dimension minus 1 and the spacing. Multiple records may be used to specify sequentially finer grid levels, but finer grids must fit within the coarsest grid. It is highly recommended that the parameters specifying the MEAD grid be the same for both the SCRF initial state and the VSCRF final state calculations.

RADIUS string real

Specifies the radius in Angstroms for an atom type. Used in fitting the ADF electronic structure to partial atomic charges and for defining the boundary between regions of low and high dielectric in MEAD. The atom types should be the same as those used in the ATOMS input block. There must be one RADIUS record for each atom type in the ATOMS input block.

INITIAL_STATE string

Specifies the filename of a binary KF file containing the ground state data for VSCRF. The KF file is created by a preliminary SCRF calculation using the SAVESTATE option. It is highly recommended that the parameters specifying the MEAD grid and the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations.

CYCLES integer

Specifies the maximum number of cycles of VSCRF to perform. Whether or not the VSCRF run has converged, it will terminate when the number of cycles exceeds the value specified by CYCLES.

TOLERANCE real

Specifies convergence criterion in kcal/mol for VSCRF. For each cycle of VSCRF the sum of ADF energy, Potential Term, and Response Term is calculated. If the difference in subsequent sums is less than the TOLERANCE value, VSCRF is considered to have converged. Defaults to 0.01.

ATOM_MAXR real

Specifies the outer atomic radius in Angstroms for the system. For each atom, grid points that lie between the atomic radius specified by the RADIUS keyword and the outer atomic radius specified here are included in charge fitting. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations. Defaults to 5.0.

CHGFIT_CONSTRAIN string

Specifies the type of constraints to be used in charge fitting. Recognized constraints are “NONE”, “CHARGE” or “DIPOLE”. NONE specifies no constraints will be applied, CHARGE specifies that only the molecular charge will be constrained and DIPOLE specifies that both the molecular charge and dipole will be constrained. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations. Default is DIPOLE.

DELATOM N Iat1 Iat2 ... IatN

Specifies which atoms should be excluded from the charge fitting procedure. N is the number of atoms to be excluded, followed by a space separated list of atoms numbers Iat1, .., IatN, that should be excluded from the charge fitting procedure. The input order in the ATOMS input block is used to identify the excluded atoms. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations.

GRID_SPACING real

Specifies the grid spacing in Angstroms for the charge fitting grid. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations. The default is 0.2.

SVD_CONDITION real

Specifies a condition number threshold for inclusion of singular values (SV) in singular value decomposition (SVD) during charge fitting. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations. The default is 0.000001.

SV_DELETE integer

Instead of using a condition number threshold for deciding which SV to include in charge fitting, SV_DELETE may be used to specify how many SV should be excluded. The smallest SV are excluded first. The default is to use a condition number threshold. If both SV_DELETE and SVD_CONDITION are specified, the SV_DELETE value will take precedence. It is highly recommended that the parameters specifying the details of the charge-fitting SVD procedure be the same for both the SCRF initial state and the VSCRF final state calculations.

EPS_OPT_SOL real

Specifies the dielectric constant of the solvent at optical frequencies for MEAD. This value is equal to the square of the index of refraction of the solvent. Defaults to the value for water: 2.0.

SOLRAD real

Specifies the radius in Angstroms of a solvent-sized probe that rolls along the surface of the molecular system to define the dielectric boundary. Defaults to a water-sized probe size of 1.4.

PROTEIN string

Use of this keyword turns on the VSCRF solinprot option by specifying the prefix of a pqr file containing the protein definition for MEAD solinprot. This option can only be used if the ground state SCRF calculation also used a protein. The filename suffix must be pqr. Pqr format contains one line per atom and begins with the ATOM keyword followed by 10 fields separated by white space and in the order: atom number, atom name, residue name, residue number, x, y and z coordinates, partial atomic charge and atomic radius. The atoms of the quantum mechanical system should NOT be included in the pqr file.

EPS_OPT_PROT real

Specifies the dielectric constant of the protein at optical frequencies for MEAD. This value is equal to the square of the index of refraction of the protein. Defaults to the value for organic liquids: 2.0.

References