6.22 Variational Two-Electron Reduced-Density-Matrix Methods

6.22.2 Theory

(May 16, 2021)

The electronic energy is an exact functional of the 1-RDM and 2-RDM

E=12pqrsDrspq2(pr|qs)+pqDqp1hpq, (6.77)

where the 1-RDM (𝐃1) and 2-RDM are represented in a given spin-orbital basis indexed by p, q, r, and s. The one-hole RDM (𝐐1), two-hole RDM (𝐐2), particle-hole RDM (𝐆2), partial three-particle RDMs (𝐓𝟏 and 𝐓𝟐), and full three-particle RDMs (𝐃3, 𝐐3, 𝐄3, 𝐅3) are linear functions of 𝐃1 and 𝐃2. 321 Fosso-Tande J. et al.
J. Chem. Theory Comput.
(2016), 12, pp. 2260.
Link
Minimizing the electronic energy with respect to 𝐃2 while enforcing the linear relations among these RDMs, the contraction and spin constraints placed on 𝐃2, and the positive semidefinite property of all RDMs constitutes a semidefinite program (SDP). The current v2RDM implementation uses a boundary-point SDP (BPSDP) algorithm to solve the SDP. 891 Povh J., Rendl F., Wiegele A.
Computing
(2006), 78, pp. 277.
Link
, 706 Malick J. et al.
SIAM J. Optim.
(2009), 20, pp. 336.
Link
, 747 Mazziotti D. A.
Phys. Rev. Lett.
(2011), 106, pp. 083001.
Link

The primal formulation of the SDP is

minimize Eprimal =𝐜T𝐱 (6.78)
such that 𝐀𝐱 =𝐛
and M(𝐱) 0.

Here, 𝐱 represents the primal solution vector, the vector 𝐜 contains all information defining the quantum system (the one- and two-electron integrals), and the mapping M(𝐱) maps the primal solution onto the set of positive semidefinite RDMs:

M(𝐱)=(𝐃10000𝐐10000𝐃20000)0. (6.79)

Additional RDMs can be included in M(𝐱), depending on the choice of N-representability conditions applied. The action of the constraint matrix, 𝐀, on 𝐱 is a compact representation of the N-representability conditions. 𝐀 maintains the appropriate mappings between each block of M(𝐱) and enforces the appropriate spin and contraction conditions. Alternatively, one could consider the dual formulation of the semidefinite problem, expressed as

maximize Edual =𝐛T𝐲 (6.80)
such that 𝐳 =𝐜-𝐀T𝐲
and M(𝐳) 0

where 𝐲 and 𝐳 are the dual solutions, and M(𝐳) is constrained to be positive semidefinite.

The BPSDP algorithm involves an iterative two-step procedure:

  1. 1.

    Solve 𝐀𝐀T𝐲=𝐀(𝐜-𝐳)+τμ(𝐛-𝐀𝐱) for 𝐲 by conjugate gradient methods.

  2. 2.

    Update 𝐱 and 𝐳 by separating 𝐔=M(μ𝐱+𝐀T𝐲-𝐜) into its positive and negative components (by diagonalization). The updated primal and dual solutions 𝐱 and 𝐳 are given by M(𝐱)=𝐔(+)/μ and M(𝐳)=-𝐔(-).

Here, τ is a step-length parameter that lies in the interval [1.0,1.6] 747 Mazziotti D. A.
Phys. Rev. Lett.
(2011), 106, pp. 083001.
Link
. The penalty parameter μ controls how strictly the primal or dual constraints are enforced and is updated dynamically according to the protocol outlined in Ref. 747. The frequency with which μ is updated is controlled by the $rem keyword RDM_MU_UPDATE_FREQUENCY. The algorithm is considered converged when the primal error ||𝐀𝐱-𝐛||, the dual error ||𝐀T𝐲-𝐜+𝐳||, and the primal/dual energy gap |Eprimal-Edual| are sufficiently small. The convergence in the primal/dual errors and the primal/dual energy gap are controlled by the $rem keywords RDM_EPS_CONVERGENCE and RDM_E_CONVERGENCE, respectively. The BPSDP algorithm scales n6 for the DQG conditions and n9 for the T1, T2, and 3POS conditions where n is the number of active orbitals in the v2RDM computation. In v2RDM-CASSCF, the molecular orbitals are optimized after a chosen number of v2RDM iterations (Steps 1. and 2. above) indicated by the $rem keyword RDM_ORBOPT_FREQUENCY.