OCC: OrbitalOptimized CoupledCluster and Møller–Plesset Perturbation Theories¶
Code author: Ugur Bozkaya
Section author: Ugur Bozkaya
Module: Keywords, PSI Variables, OCC
Module: Keywords, PSI Variables, DFOCC
Introduction¶
Orbitaloptimized methods have several advantages over their nonoptimized counterparts. Once the orbitals are optimized, the wave function will obey the Hellmann–Feynman theorem for orbital rotation parameters. Therefore, there is no need for orbital response terms in the evaluation of analytic gradients. In other words, it is unnecessary to solve the first order coupledperturbed CC and manybody perturbation theory (MBPT) equations. Further, computation of oneelectron properties is easier because there are no response contributions to the particle density matrices (PDMs). Moreover, active space approximations can be readily incorporated into the CC methods [Krylov:2000:vod]. Additionally, orbitaloptimized coupledcluster avoids spurious secondorder poles in its response function, and its transition dipole moments are gauge invariant [Pedersen:1999:od].
Another advantage is that the orbitaloptimized methods do not suffer from artifactual symmetrybreaking instabilities [Crawford:1997:instability], [Sherrill:1998:od], [Bozkaya:2011:omp2], and [Bozkaya:2011:omp3]. Furthermore, Kurlancheek and HeadGordon [Kurlancek:2009] demonstrated that first order properties such as forces or dipole moments are discontinuous along nuclear coordinates when such a symmetry breaking occurs. They also observed that although the energy appears well behaved, the MP2 method can have natural occupation numbers greater than 2 or less than 0, hence may violate the Nrepresentability condition. They further discussed that the orbital response equations generally have a singularity problem at the unrestriction point where spinrestricted orbitals become unstable to unrestriction. This singularity yields to extremely large or small eigenvalues of the oneparticle density matrix (OPDM). These abnormal eigenvalues may lead to unphysical molecular properties such as vibrational frequencies. However, orbitaloptimized MP2 (also MP3) will solve this Nrepresentability problem by disregarding orbital response contribution of oneparticle density matrix.
Although the performance of coupledcluster singles and doubles (CCSD) and orbitaloptimized CCD (OD) is similar, the situation is different in the case of triples corrections, especially at stretched geometries [Bozkaya:2012:odtl]. Bozkaya and Schaefer demonstrated that orbitaloptimized coupled cluster based triple corrections, especially those of asymmetrics, provide significantly better potential energy curves than CCSD based triples corrections.
A lot of the functionality in OCC has been enabled with Density Fitting (DF) and Cholesky Decomposition (CD) techniques, which can greatly speed up calculations and reduce memory requirements for typically negligible losses in accuracy.
NOTE: As will be discussed later, all methods with orbitaloptimization functionality have nonorbital
optimized counterparts. Consequently, there arise two possible ways to call densityfitted MP2. In most
cases, users should prefer the DFMP2 code described in the DFMP2 section because it is
faster. If gradients are needed (like in a geometry optimization), then the procedures outlined hereafter
should be followed.
In general, choose the desired method, reference, and ERI type (e.g.,
set reference uhf
, set mp2_type df
, opt('mp2')
) and the most
efficient module will be selected automatically, according to
Crossmodule Redundancies.
Thus, there arise a few categories of method, each with corresponding input keywords:
Orbitaloptimized MP and CC methods with conventional integrals (OMP Methods OCC keywords)
Orbitaloptimized MP and CC methods with DF and CD integrals (OMP Methods DFOCC keywords)
Nonorbitaloptimized MP and CC methods with conventional integrals (MP/CC Methods OCC keywords)
Nonorbitaloptimized MP and CC methods with DF and CD integrals (MP/CC Methods DFOCC keywords)
Theory¶
What follows is a very basic description of orbitaloptimized Møller–Plesset perturbation theory as implemented in PSI4. We will follow our previous presentations ([Bozkaya:2011:omp2], [Bozkaya:2011:omp3], and [Bozkaya:2012:odtl])
The orbital variations may be expressed by means of an exponential unitary operator
where \(\hat{K}\) is the orbital rotation operator
The effect of the orbital rotations on the MO coefficients can be written as
where \({\bf C^{(0)}}\) is the initial MO coefficient matrix and \({\bf C({\bf \kappa})}\) is the new MO coefficient matrix as a function of \({\bf \kappa}\). Now, let us define a variational energy functional (Lagrangian) as a function of \({\bf \kappa}\)
OMP2
OMP3
OLCCD
where subscript c means only connected diagrams are allowed, and \(\hat{H}^{\kappa}\), \(\hat{f}_{N}^{\kappa}\), and \(\hat{W}_{N}^{\kappa}\) defined as
where \(\hat{f}_{N}\), and \(\hat{W}_{N}\) are the one and twoelectron components of normalordered Hamiltonian. Then, first and second derivatives of the energy with respect to the \({\bf \kappa}\) parameter at \({\bf \kappa} = 0\)
Then the energy can be expanded up to secondorder as follows
where \({\bf w}\) is the MO gradient vector, \({\bf \kappa}\) is the MO rotation vector, and \({\bf A}\) is the MO Hessian matrix. Therefore, minimizing the energy with respect to \({\bf \kappa}\) yields
This final equation corresponds to the usual NewtonRaphson step.
OREMP
The REMP hybrid perturbation theory is a constrained mixture of the Møller–Plesset perturbation theory and the Retaining the Excitation degree perturbation theory([Fink:2006:RE], [Behnle:2019:REMP]). The mixing ratio is determined by the parameter :math’:A:
Technically, the second order of RE corresponds to LCCD for RHF and UHF references. REMP2 and its orbitaloptimized variant OREMP2 are thus straightforward to implement in a (O)LCCD program by appropriate scaling of residual vector contributions and density matrices.
OMP3 [Bozkaya:2011:omp3] , [Bozkaya:2013:omp3], and [Bozkaya:2013:omp3grad]
OMP2.5 [Bozkaya:2011:omp3]
OLCCD [Bozkaya:2013:ocepa]
LCCD [Bozkaya:2013:ocepa]
OREMP2 [Behnle:2021:OREMP], and [Behnle:2022:OREMP]
Convergence Problems¶
For problematic openshell systems, we recommend to use the ROHF or DFT orbitals as an initial guess for orbitaloptimized methods. Both ROHF and
DFT orbitals may provide better initial guesses than UHF orbitals, hence convergence may be significantly speeded up with ROHF or DFT orbitals.
In order to use ROHF orbitals, simply set reference rohf
. For DFT orbitals, set reference uks
and set dft_functional b3lyp
. Of
course users can use any DFT functional available in PSI4.
Methods¶
The various orbitaloptimized methods supported by the OCC/DFOCC
modules in PSI4 are summarized in Table OCC OO
Methods and detailed in Table OCC
OO Capabilities. Note that while two
separate libraries OCC (conventional integrals CONV
) and DFOCC
(densityfitted DF
and Choleskydecomposed CD
) together provide
the methods described on this page, they are controlled through one
QC_MODULE value OCC
. Without set qc_module occ
,
these methods may default to implementations in other modules based on efficiency considerations.
name 
calls method 
OO 

omp2 
orbitaloptimized secondorder MP perturbation theory 
E/G 
omp2.5 
orbitaloptimized average of MP2 and MP3 
E/G 
omp3 
orbitaloptimized thirdorder MP perturbation theory 
E/G 
oremp2 
orbitaloptimized secondorder REMP hybrid PT 
E/G 
olccd 
orbitaloptimized linear coupled cluster doubles 
E/G 
◻ ◻ name ↓ → ◻ ◻ 
◻ ◻ type[1] ↓ → 
QC_MODULE=OCC Capabilities 


Restricted (RHF) 
Unrestricted (UHF) 
Restricted Open (ROHF) 

CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 

A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 

omp2[4] 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 

omp2.5[4] 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 

omp3[4] 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 

oremp2[4] 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 

olccd[4] 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
name 
calls method 
Energy 
Gradient 

scsomp3 
SpinComponent Scaled OrbitalOptimized MP3 
RHF/UHF/ROHF/RKS/UKS 
— 
sosomp3 
SpinOpposite Scaled OrbitalOptimized MP3 
RHF/UHF/ROHF/RKS/UKS 
— 
scs(n)omp3 
A special version of SCSOMP3 for nucleobase interactions 
RHF/UHF/ROHF/RKS/UKS 
— 
scsomp3vdw 
A special version of SCSOMP3 (from ethene dimers) 
RHF/UHF/ROHF/RKS/UKS 
— 
sospiomp3 
A special version of SOSOMP3 for \(\pi\)systems 
RHF/UHF/ROHF/RKS/UKS 
— 
scsomp2 
SpinComponent Scaled OrbitalOptimized MP2 
RHF/UHF/ROHF/RKS/UKS 
— 
sosomp2 
SpinOpposite Scaled OrbitalOptimized MP2 
RHF/UHF/ROHF/RKS/UKS 
— 
scs(n)omp2 
A special version of SCSOMP2 for nucleobase interactions 
RHF/UHF/ROHF/RKS/UKS 
— 
scsomp2vdw 
A special version of SCSOMP2 (from ethene dimers) 
RHF/UHF/ROHF/RKS/UKS 
— 
sospiomp2 
A special version of SOSOMP2 for \(\pi\)systems 
RHF/UHF/ROHF/RKS/UKS 
— 
Basic OCC Keywords¶
E_CONVERGENCE¶
Convergence criterion for energy. See Table PostSCF Convergence for default convergence criteria for different calculation types.
Type: conv double
Default: 1e6
R_CONVERGENCE¶
Convergence criterion for amplitudes (residuals).
Type: conv double
Default: 1e5
RMS_MOGRAD_CONVERGENCE¶
Convergence criterion for RMS orbital gradient. If this keyword is not set by the user, OCC will estimate and use a value required to achieve the desired E_CONVERGENCE The listed default will be used for the default value of E_CONVERGENCE
Type: conv double
Default: 1e4
MAX_MOGRAD_CONVERGENCE¶
Convergence criterion for maximum orbital gradient. If this keyword is not set by the user, OCC will estimate and use a value required to achieve the desired E_CONVERGENCE The listed default will be used for the default value of E_CONVERGENCE
Type: conv double
Default: 1e4
MO_MAXITER¶
Maximum number of iterations to determine the orbitals
Type: integer
Default: 50
WFN_TYPE¶
Type of the wavefunction.
Type: string
Possible Values: OMP2, OMP3, OCEPA, OMP2.5, REMP, OREMP
Default: OMP2
ORB_OPT¶
Do optimize the orbitals?
Type: boolean
Default: true
Advanced OCC Keywords¶
OPT_METHOD¶
The optimization algorithm. Modified SteepestDescent (MSD) takes a NewtonRaphson (NR) step with a crude approximation to diagonal elements of the MO Hessian. The ORB_RESP option obtains the orbital rotation parameters with a crude approximation to all elements of the MO Hessian. Additionally, for both methods a DIIS extrapolation will be performed with the DO_DIIS = TRUE option.
Type: string
Possible Values: MSD, ORB_RESP
Default: MSD
ORTH_TYPE¶
The algorithm for orthogonalization of MOs
Type: string
Possible Values: GS, MGS
Default: MGS
MP2_OS_SCALE¶
Removed in 1.4. Will raise an error in 1.5.
Type: double
Default: 6.0
MP2_SS_SCALE¶
Removed in 1.4. Will raise an error in 1.5.
Type: double
Default: 1.0
MP2_SOS_SCALE¶
Removed in 1.4. Will raise an error in 1.5.
Type: double
Default: 1.3
MP2_SOS_SCALE2¶
Removed in 1.4. Will raise an error in 1.5.
Type: double
Default: 1.2
NAT_ORBS¶
Do compute natural orbitals?
Type: boolean
Default: false
OCC_ORBS_PRINT¶
Do print OCC orbital energies?
Type: boolean
Default: false
TPDM_ABCD_TYPE¶
How to take care of the TPDM VVVVblock. The COMPUTE option means it will be computed via an IC/OOC algorithm. The DIRECT option (default) means it will not be computed and stored, instead its contribution will be directly added to GeneralizedFock Matrix.
Type: string
Possible Values: DIRECT, COMPUTE
Default: DIRECT
DO_DIIS¶
Do apply DIIS extrapolation?
Type: boolean
Default: true
DO_LEVEL_SHIFT¶
Removed in 1.4. Will raise an error in 1.5.
Type: boolean
Default: true
Basic DFOCC Keywords¶
E_CONVERGENCE¶
Convergence criterion for energy. See Table PostSCF Convergence for default convergence criteria for different calculation types.
Type: conv double
Default: 1e6
R_CONVERGENCE¶
Convergence criterion for amplitudes (residuals).
Type: conv double
Default: 1e5
RMS_MOGRAD_CONVERGENCE¶
Convergence criterion for RMS orbital gradient. If this keyword is not set by the user, DFOCC will estimate and use a value required to achieve the desired E_CONVERGENCE The listed default will be used for the default value of E_CONVERGENCE
Type: conv double
Default: 1e4
MAX_MOGRAD_CONVERGENCE¶
Convergence criterion for maximum orbital gradient. If this keyword is not set by the user, DFOCC will estimate and use a value required to achieve the desired E_CONVERGENCE The listed default will be used for the default value of E_CONVERGENCE
Type: conv double
Default: 1e4
MO_MAXITER¶
Maximum number of iterations to determine the orbitals
Type: integer
Default: 100
ORB_OPT¶
Do optimize the orbitals?
Type: boolean
Default: true
Advanced DFOCC Keywords¶
OPT_METHOD¶
The orbital optimization algorithm. Presently quasiNewtonRaphson algorithm available with several Hessian * options.
Type: string
Possible Values: QNR
Default: QNR
HESS_TYPE¶
Type of the MO Hessian matrix
Type: string
Possible Values: APPROX_DIAG, APPROX_DIAG_EKT, APPROX_DIAG_HF, HF
Default: HF
MO_DIIS_NUM_VECS¶
Number of vectors used in orbital DIIS
Type: integer
Default: 6
ORTH_TYPE¶
The algorithm for orthogonalization of MOs
Type: string
Possible Values: GS, MGS
Default: MGS
DO_DIIS¶
Do apply DIIS extrapolation?
Type: boolean
Default: true
DO_LEVEL_SHIFT¶
Do apply level shifting?
Type: boolean
Default: true
Conventional (NonOO) CoupledCluster and Møller–Plesset Perturbation Theories¶
The various nonorbitaloptimized methods supported by the OCC/DFOCC
modules in PSI4 are summarized in Table OCC nonOO
Methods and detailed in Table OCC
nonOO Capabilities. Note that while two
separate libraries OCC (conventional integrals CONV
) and DFOCC
(densityfitted DF
and Choleskydecomposed CD
) together provide
the methods described on this page, they are controlled through one
QC_MODULE value OCC
. Without set qc_module occ
,
these methods may default to implementations in other modules based on efficiency considerations.
Starting in v1.4, MP2.5 and MP3 default to the densityfit algorithm. Set MP_TYPE to CONV
to get previous behavior.
Publications resulting from the use of the nonOO CC codes should cite the following publications:
MP2 [Bozkaya:2011:omp2], [Bozkaya:2013:omp2grad], and [Bozkaya:2014:dfomp2grad]
MP3 [Bozkaya:2011:omp3] , [Bozkaya:2013:omp3], [Bozkaya:2013:omp3grad], [Bozkaya:2016:dfomp3], and [Bozkaya:2018:dfomp3grad]
MP2.5 [Bozkaya:2011:omp3], [Bozkaya:2016:dfomp3], and [Bozkaya:2018:dfomp3grad]
CCSD(T) [Bozkaya:2017:dfccsdtgrad]
CCSD(AT) [Bozkaya:2016:dfccsdat]
name 
calls method 
plain 
FNO 

mp2 
secondorder MP perturbation theory 
E/G 
n/a 
mp2.5 
average of MP2 and MP3 
E/G 

mp3 
thirdorder MP perturbation theory 
E/G 

remp2 
secondorder retainingtheexcitationdegree MP hybrid PT 
E 

lccd 
linear coupled cluster doubles 
E/G 

ccd 
coupled cluster doubles 
E/G 

ccsd 
coupled cluster singles and doubles 
E/G 

ccsd(t) 
coupled cluster singles and doubles with perturbative triples 
E/G 

accsd(t) 
CCSD with asymmetric perturbative triples 
E 
◻ ◻ name ↓ → ◻ ◻ 
◻ ◻ type[5] ↓ → 
QC_MODULE=OCC Capabilities 


Restricted (RHF) 
Unrestricted (UHF) 
Restricted Open (ROHF) 

CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 
CV 
DF 
CD 

A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 
A 
F 

mp2 
✓̲ 
✓̲ 
✓ 
✓ 
✓̲ 
✓̲ 
✓̲ 
✓ 
✓ 
✓̲ 
✓̲ 
✓ 
✓ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓ 
✓ 
✓̲ 
✓̲ 

mp2.5 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 

mp3 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 

remp2 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 

lccd 
✓ 
✓ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 
✓̳ 
✓̳ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̳ 
✓̲ 
✓̲ 

ccd 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 

ccsd 
✓ 
✓ 
✓ 
✓ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
✓̲ 

ccsd(t) 
✓ 
✓ 
✓ 
✓ 
✓̲ 
✓̲ 

accsd(t)[8] 
✓̲ 
✓̲ 
✓̲ 
✓̲ 
Algorithm type selection keyword below. Values to the right: conventional CV
, densityfitted DF
, and Choleskydecomposed CD
.
Active orbital values to the right: allelectron A
and frozencore F
.
Finite difference gradients are not marked explicitly by “∷”, but the capability can be gleaned from the energy availability.
aCCSD(T) also known as CCSD(aT), LambdaCCSD(T), and CCSD(T)_L