In Q-Chem, the unrestricted and restricted GVB methods are implemented with a
resolution of the identity (RI) algorithm that makes them computationally very
efficient.^{937, 851} 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.^{535} 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 $i$th pair exciting into the virtual orbitals of
the $j$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.

GVB_ORB_MAX_ITER

Controls the number of orbital iterations allowed in GVB-CC calculations.
Some jobs, particularly unrestricted PP jobs can require 500–1000 iterations.

TYPE:

INTEGER

DEFAULT:

256

OPTIONS:

User-defined number of iterations.

RECOMMENDATION:

Default is typically adequate, but some jobs, particularly UPP jobs, can
require 500–1000 iterations if converged tightly.

GVB_ORB_CONV

The GVB-CC wave function is considered converged when the root-mean-square
orbital gradient and orbital step sizes are less than
${10}^{-\mathrm{GVB}\mathrm{\_}\mathrm{ORB}\mathrm{\_}\mathrm{CONV}}$. Adjust THRESH simultaneously.

TYPE:

INTEGER

DEFAULT:

5

OPTIONS:

$n$
User-defined

RECOMMENDATION:

Use 6 for PP(2) jobs or geometry optimizations.
Tighter convergence (*i.e.* 7 or higher) cannot always be reliably achieved.

GVB_ORB_SCALE

Scales the default orbital step size by $n$/1000.

TYPE:

INTEGER

DEFAULT:

1000
Corresponding to 100%

OPTIONS:

$n$
User-defined, 0–1000

RECOMMENDATION:

Default is usually fine, but for some stretched geometries it
can help with convergence to use smaller values.

GVB_AMP_SCALE

Scales the default orbital amplitude iteration step size by $n$/1000 for IP/RCC.
PP amplitude equations are solved analytically, so this parameter does not
affect PP.

TYPE:

INTEGER

DEFAULT:

1000
Corresponding to 100%

OPTIONS:

$n$
User-defined, 0–1000

RECOMMENDATION:

Default is usually fine, but in some highly-correlated
systems it can help with convergence to use smaller values.

GVB_RESTART

Restart a job from previously-converged GVB-CC orbitals.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE/FALSE

RECOMMENDATION:

Useful when trying to converge to the same GVB
solution at slightly different geometries, for example.

GVB_REGULARIZE

Coefficient for GVB_IP exchange type amplitude regularization
to improve the convergence of the amplitude equations especially
for spin-unrestricted amplitudes near dissociation. This is the
leading coefficient for an amplitude dampening term
$-(c/10000)({e}^{{t}_{ij}^{p}}-1)/({e}^{1}-1)$

TYPE:

INTEGER

DEFAULT:

0
For restricted
1
For unrestricted

OPTIONS:

$c$
User-defined

RECOMMENDATION:

Should be increased if unrestricted amplitudes do not converge or
converge slowly at dissociation. Set this to zero to remove
all dynamically-valued amplitude regularization.

GVB_POWER

Coefficient for GVB_IP exchange type amplitude regularization
to improve the convergence of the amplitude equations especially
for spin-unrestricted amplitudes near dissociation. This is the
leading coefficient for an amplitude dampening term included in
the energy denominator: -($c$/10000)$({e}^{{t}_{ij}^{p}}-1)/({e}^{1}-1)$

TYPE:

INTEGER

DEFAULT:

6

OPTIONS:

$p$
User-defined

RECOMMENDATION:

Should be decreased if unrestricted amplitudes do not converge or
converge slowly at dissociation, and should be kept even valued.

GVB_SHIFT

Value for a statically valued energy shift in the energy
denominator used to solve the coupled cluster amplitude equations, $n$/10000.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined

RECOMMENDATION:

Default is fine, can be used in lieu of the dynamically
valued amplitude regularization if it does not aid convergence.

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.^{534} 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.

GVB_LOCAL

Sets the localization scheme used in the initial guess wave function.

TYPE:

INTEGER

DEFAULT:

2
Pipek-Mezey orbitals

OPTIONS:

0
No Localization
1
Boys localized orbitals
2
Pipek-Mezey orbitals

RECOMMENDATION:

Different initial guesses can sometimes lead to different solutions.
It can be helpful to try both to ensure the global minimum has been found.

GVB_DO_SANO

Sets the scheme used in determining the active virtual orbitals
in a Unrestricted-in-Active Pairs GVB calculation.

TYPE:

INTEGER

DEFAULT:

2

OPTIONS:

0
No localization or Sano procedure
1
Only localizes the active virtual orbitals
2
Uses the Sano procedure

RECOMMENDATION:

Different initial guesses can sometimes lead to different solutions. Disabling
sometimes can aid in finding more non-local solutions for the orbitals.

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 (${n}_{\mathrm{active}}={n}_{\beta})$. Or
alternatively one could make only formal bonding electron pairs active
(${n}_{\mathrm{active}}={N}_{\mathrm{STO}-3\mathrm{G}}-{n}_{\alpha})$. Or in some cases, one might
want only the most reactive electron pair to be active (${n}_{\mathrm{active}}=$1).
Clearly making physically appropriate choices for this variable is essential
for obtaining physically appropriate results!

GVB_N_PAIRS

Alternative to CC_REST_OCC and CC_REST_VIR for setting
active space size in GVB and valence coupled cluster methods.

TYPE:

INTEGER

DEFAULT:

PP active space (1 occ and 1 virt for each valence electron pair)

OPTIONS:

$n$
user-defined

RECOMMENDATION:

Use the default unless one wants to study a special active space. When using
small active spaces, it is important to ensure that the proper orbitals are
incorporated in the active space. If not, use the *$reorder_mo* feature to adjust
the SCF orbitals appropriately.

GVB_PRINT

Controls the amount of information printed during a GVB-CC job.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined

RECOMMENDATION:

Should never need to go above 0 or 1.

GVB_TRUNC_OCC

Controls how many pairs’ occupied orbitals are truncated from the GVB active space.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined

RECOMMENDATION:

This allows for asymmetric GVB active spaces removing the $n$
lowest energy occupied orbitals from the GVB active space
while leaving their paired virtual orbitals in the active space.
Only the models including the SIP and DIP amplitudes (ie NP and 2P)
benefit from this all other models this equivalent to just
reducing the total number of pairs.

GVB_TRUNC_VIR

Controls how many pairs’ virtual orbitals are truncated from the GVB active space.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined

RECOMMENDATION:

This allows for asymmetric GVB active spaces removing
the $n$ highest energy occupied orbitals from the GVB
active space while leaving their paired virtual orbitals
in the active space. Only the models including the SIP
and DIP amplitudes (ie NP and 2P) benefit from this all
other models this equivalent to just reducing the total number of pairs.

GVB_REORDER_PAIRS

Tells the code how many GVB pairs to switch around.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
$0\le n\le 5$

RECOMMENDATION:

This allows for the user to change the order the active
pairs are placed in after the orbitals are read in or are
guessed using localization and the Sano procedure. Up to 5
sequential pair swaps can be made, but it is best to leave this alone.

GVB_REORDER_1

Tells the code which two pairs to swap first.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined XXXYYY

RECOMMENDATION:

This is in the format of two 3-digit pair indices that
tell the code to swap pair XXX with YYY, for example swapping
pair 1 and 2 would get the input 001002. Must be specified
in GVB_REORDER_PAIRS $\ge $ 1.

GVB_REORDER_2

Tells the code which two pairs to swap second.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined XXXYYY

RECOMMENDATION:

This is in the format of two 3-digit pair indices that tell the code to swap
pair XXX with YYY, for example swapping pair 1 and 2 would get the input
001002. Must be specified in GVB_REORDER_PAIRS $\ge $ 2.

GVB_REORDER_3

Tells the code which two pairs to swap third.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined XXXYYY

RECOMMENDATION:

This is in the format of two 3-digit pair indices that tell the code to swap
pair XXX with YYY, for example swapping pair 1 and 2 would get the input
001002. Must be specified in GVB_REORDER_PAIRS $\ge $ 3.

GVB_REORDER_4

Tells the code which two pairs to swap fourth.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined XXXYYY

RECOMMENDATION:

This is in the format of two 3-digit pair indices that tell the code to swap
pair XXX with YYY, for example swapping pair 1 and 2 would get the input
001002. Must be specified in GVB_REORDER_PAIRS $\ge $ 4.

GVB_REORDER_5

Tells the code which two pairs to swap fifth.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$
User-defined XXXYYY

RECOMMENDATION:

This is in the format of two 3-digit pair indices that tell the code to swap
pair XXX with YYY, for example swapping pair 1 and 2 would get the input
001002. Must be specified in GVB_REORDER_PAIRS $\ge $ 5.

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,^{532, 533} 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).

GVB_SYMFIX

Should GVB use a symmetry breaking fix.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0
no symmetry breaking fix
1
symmetry breaking fix with virtual orbitals spanning the active space
2
symmetry breaking fix with virtual orbitals spanning the whole virtual space

RECOMMENDATION:

It is best to stick with type 1 to get a symmetry breaking correction
with the best results coming from CORRELATION=NP and
GVB_SYMFIX = 1.

GVB_SYMPEN

Sets the pre-factor for the amplitude regularization term for the SB amplitudes.

TYPE:

INTEGER

DEFAULT:

160

OPTIONS:

$\gamma $
User-defined

RECOMMENDATION:

Sets the pre-factor for the amplitude regularization term for
the SB amplitudes: $-(\gamma /1000)({e}^{(c*100)*{t}^{2}}-1)$.

GVB_SYMSCA

Sets the weight for the amplitude regularization term for the SB amplitudes.

TYPE:

INTEGER

DEFAULT:

125

OPTIONS:

$c$
User-defined

RECOMMENDATION:

Sets the weight for the amplitude regularization term for the SB amplitudes:
$-(\gamma /1000)({e}^{(c*100)*{t}^{2}}-1)$.

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${}_{3h}$ symmetry over
D${}_{6h}$ by 3 kcal/mol (with a 2o distortion), while in IP, this difference
is reduced to 0.5 kcal/mol and less than 1o.^{932}
Likewise the allyl radical breaks symmetry in the unrestricted PP
model,^{79} 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 $4s3d$ into $4s4p3d$) or possibly hyper-valent molecules (expansion from
$3s3p$ into $3s3p3d$). 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.