In Q-Chem, the unrestricted and restricted GVB methods are implemented with a
resolution of the identity (RI) algorithm that makes them computationally very
J. Chem. Phys.
(2002), 117, pp. 9190. , 1052 J. Chem. Theory Comput.
(2006), 2, pp. 300. They can be applied to systems with more than 100 active electrons, and both energies and analytical gradients are available. These methods are requested via the standard CORRELATION keyword. If AUX_BASIS is not specified, the calculation uses four-center two-electron integrals by default. Much faster auxiliary basis algorithms (see Section 6.6 for an introduction), which are used for the correlation energy (not the reference SCF energy), can be enabled by specifying a valid string for AUX_BASIS. The example below illustrates a simple IP calculation.
$molecule 0 1 H F 1 1.0 $end $rem JOBTYPE opt CORRELATION gvb_ip BASIS cc-pVDZ AUX_BASIS rimp2-cc-pVDZ $end
If further improvement in the orbitals are needed, the GVB_SIP, GVB_DIP, OP, NP
and 2P models are also included.
J. Phys. Chem. A
(2010), 114, pp. 2930. The GVB_SIP model includes all the amplitudes of GVB_IP plus a set of quadratic amplitudes the represent the single ionization of a pair. The GVB_DIP model includes the GVB_SIP models amplitudes and the doubly ionized pairing amplitudes which are analogous to the correlation of the occupied electrons of the th pair exciting into the virtual orbitals of the th pair. These two models have the implementation limit of no analytic orbital gradient, meaning that a slow finite differences calculation must be performed to optimize their orbitals, or they must be computed using orbitals from a different method. The 2P model is the same as the GVB_DIP model, except it only allows the amplitudes to couple via integrals that span only two pairs. This allows for a fast implementation of it’s analytic orbital gradient and enables the optimization of it’s own orbitals. The OP method is like the 2P method except it removes the “direct”-like IP amplitudes and all of the same-spin amplitudes. The NP model is the GVB_IP model with the DIP amplitudes. This model is the one that works best with the symmetry breaking corrections that will be discussed later. All GVB methods except GVB_SIP and GVB_DIP have an analytic nuclear gradient implemented for both regular and RI four-center two-electron integrals.
There are often considerable challenges in converging the orbital optimization associated with these GVB-type calculations. The situation is somewhat analogous to SCF calculations but more severe because there are more orbital degrees of freedom that affect the energy (for instance, mixing occupied active orbitals amongst each other, mixing active virtual orbitals with each other, mixing core and active occupied, mixing active virtual and inactive virtual). Furthermore, the energy changes associated with many of these new orbital degrees of freedom are rather small and delicate. As a consequence, in cases where the correlations are strong, these GVB-type jobs often require many more iterations than the corresponding GDM calculations at the SCF level. This is a reflection of the correlation model itself. To deal with convergence issues, a number of REM values are available to customize the calculations, as listed below.
Another issue that a user of these methods should be aware of is the fact that
there is a multiple minimum challenge associated with GVB calculations. In SCF
calculations it is sometimes possible to converge to more than one set of
orbitals that satisfy the SCF equations at a given geometry. The same problem
can arise in GVB calculations, and based on our experience to date, the problem
in fact is more commonly encountered in GVB calculations than in SCF
calculations. A user may therefore want to (or have to!) tinker with the
initial guess used for the calculations. One way is to set GVB_RESTART
= TRUE (see above), to replace the default initial guess (the converged SCF
orbitals which are then localized). Another way is to change the localized
orbitals that are used in the initial guess, which is controlled by the
GVB_LOCAL variable, described below. Sometimes different localization
criteria, and thus different initial guesses, lead to different converged
solutions. Using the new amplitude regularization keywords enables some control
over the solution GVB optimizes.
J. Chem. Phys.
(2009), 130, pp. 184113. A calculation can be performed with amplitude regularization to find a desired solution, and then the calculation can be rerun with GVB_RESTART = TRUE and the regularization turned off to remove the energy penalty of regularization.
Other $rem variables relevant to GVB calculations are given below. It is possible to explicitly set the number of active electron pairs using the GVB_N_PAIRS variable. The default is to make all valence electrons active. Other reasonable choices are certainly possible. For instance all electron pairs could be active (. Or alternatively one could make only formal bonding electron pairs active (. Or in some cases, one might want only the most reactive electron pair to be active (1). Clearly making physically appropriate choices for this variable is essential for obtaining physically appropriate results!
It is known that symmetry breaking of the orbitals to favor localized solutions
over non-local solutions is an issue with GVB methods in general. A combined
coupled-cluster perturbation theory approach to solving symmetry breaking (SB)
using perturbation theory level double amplitudes that connect up to three
pairs has been examined in the literature,
J. Chem. Phys.
(2008), 128, pp. 024107. , 640 Mol. Phys.
(2008), 106, pp. 2309. and it seems to alleviate the SB problem to a large extent. It works in conjunction with the GVB_IP, NP, and 2P levels of correlation for both restricted and unrestricted wave functions (barring that there is no restricted implementation of the 2P model, but setting GVB_DO_ROHF to the same number as the number of pairs in the system is equivalent).
We have already mentioned a few issues associated with the GVB calculations:
the neglect of dynamic correlation [which can be remedied with PP(2)], the
convergence challenges and the multiple minimum issues. Another weakness of
these GVB methods is the occasional symmetry-breaking artifacts that are a
consequence of the limited number of retained pair correlation amplitudes. For
example, benzene in the PP approximation prefers D symmetry over
D by 3 kcal/mol (with a 2o distortion), while in IP, this difference
is reduced to 0.5 kcal/mol and less than 1o.
Chem. Phys. Lett.
(2000), 317, pp. 575. Likewise the allyl radical breaks symmetry in the unrestricted PP model, 89 J. Phys. Chem. A
(2005), 109, pp. 9183. although to a lesser extent than in restricted open shell HF. Another occasional weakness is the limitation to the perfect pairing active space, which is not necessarily appropriate for molecules with expanded valence shells, such as in some transition metal compounds (e.g. expansion from into ) or possibly hyper-valent molecules (expansion from into ). The singlet strongly orthogonal geminal method (see the next section) is capable of dealing with expanded valence shells and could be used for such cases. The perfect pairing active space is satisfactory for most organic and first row inorganic molecules.
To summarize, while these GVB methods are powerful and can yield much insight when used properly, they do have enough pitfalls for not to be considered true “black box” methods.