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.
1218
J. Chem. Phys.
(2002),
117,
pp. 9190.
Link
,
1117
J. Chem. Theory Comput.
(2006),
2,
pp. 300.
Link
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.
679
J. Phys. Chem. A
(2010),
114,
pp. 2930.
Link
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.
GVB_ORB_MAX_ITER
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
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
. Adjust THRESH simultaneously.
TYPE:
INTEGER
DEFAULT:
5
OPTIONS:
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
GVB_ORB_SCALE
Scales the default orbital step size by /1000.
TYPE:
INTEGER
DEFAULT:
1000
Corresponding to 100%
OPTIONS:
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
GVB_AMP_SCALE
Scales the default orbital amplitude iteration step size by /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:
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
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
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
TYPE:
INTEGER
DEFAULT:
0
For restricted
1
For unrestricted
OPTIONS:
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
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: -(/10000)
TYPE:
INTEGER
DEFAULT:
6
OPTIONS:
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
GVB_SHIFT
Value for a statically valued energy shift in the energy
denominator used to solve the coupled cluster amplitude equations, /10000.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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.
678
J. Chem. Phys.
(2009),
130,
pp. 184113.
Link
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
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
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 (. 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!
GVB_N_PAIRS
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:
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
GVB_PRINT
Controls the amount of information printed during a GVB-CC job.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
User-defined
RECOMMENDATION:
Should never need to go above 0 or 1.
GVB_TRUNC_OCC
GVB_TRUNC_OCC
Controls how many pairs’ occupied orbitals are truncated from the GVB active space.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
User-defined
RECOMMENDATION:
This allows for asymmetric GVB active spaces removing the
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 (i.e. NP and 2P)
benefit from this all other models this equivalent to just
reducing the total number of pairs.
GVB_TRUNC_VIR
GVB_TRUNC_VIR
Controls how many pairs’ virtual orbitals are truncated from the GVB active space.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
User-defined
RECOMMENDATION:
This allows for asymmetric GVB active spaces removing
the 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 (i.e. NP and 2P) benefit from this all
other models this equivalent to just reducing the total number of pairs.
GVB_REORDER_PAIRS
GVB_REORDER_PAIRS
Tells the code how many GVB pairs to switch around.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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
GVB_REORDER_1
Tells the code which two pairs to swap first.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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 1.
GVB_REORDER_2
GVB_REORDER_2
Tells the code which two pairs to swap second.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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 2.
GVB_REORDER_3
GVB_REORDER_3
Tells the code which two pairs to swap third.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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 3.
GVB_REORDER_4
GVB_REORDER_4
Tells the code which two pairs to swap fourth.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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 4.
GVB_REORDER_5
GVB_REORDER_5
Tells the code which two pairs to swap fifth.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
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 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,
676
J. Chem. Phys.
(2008),
128,
pp. 024107.
Link
,
677
Mol. Phys.
(2008),
106,
pp. 2309.
Link
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
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
GVB_SYMPEN
Sets the pre-factor for the amplitude regularization term for the SB amplitudes.
TYPE:
INTEGER
DEFAULT:
160
OPTIONS:
User-defined
RECOMMENDATION:
Sets the pre-factor for the amplitude regularization term for
the SB amplitudes: .
GVB_SYMSCA
GVB_SYMSCA
Sets the weight for the amplitude regularization term for the SB amplitudes.
TYPE:
INTEGER
DEFAULT:
125
OPTIONS:
User-defined
RECOMMENDATION:
Sets the weight for the amplitude regularization term for the SB amplitudes:
.
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.
1213
Chem. Phys. Lett.
(2000),
317,
pp. 575.
Link
Likewise the allyl radical breaks symmetry in the unrestricted PP
model,
93
J. Phys. Chem. A
(2005),
109,
pp. 9183.
Link
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.