10. The pCCD module¶

Two-electron functions, called geminals, can be used to incorporate electron correlation effects into the many-particle wave function. PyBEST supports a special type of geminal-based wave function models, the antisymmetric product of 1-reference orbital geminals (pCCD), which is equivalent to pair-coupled cluster doubles. The pCCD wave function effectively parameterizes the doubly occupied configuration interaction wave function (DOCI), but requires only mean-field computational cost in contrast to the factorial scaling of traditional DOCI implementations [limacher2013]. Currently, the pCCD module is limited to closed-shell systems only.

10.1. The pCCD model¶

The pCCD wave function ansatz can be rewritten in terms of one-particle functions as a fully general pair-coupled-cluster wave function,

(10.1)¶\[\vert \textrm{pCCD}\rangle = \exp(\sum_{ia} c_i^a a^\dagger_a a^\dagger_{\bar{a}} a_{\bar{i}} a_i) \vert \Phi_0 \rangle\]

where \(a_p^{\dagger}\), \(a_{\bar{p}}^{\dagger}\), and \(a_p\), \(a_{\bar{p}}\) are the electron creation and annihilation operators and \(p\) and \(\bar{p}\) denote \(\alpha\) and \(\beta\) spins, respectively. \(\vert \Phi_0 \rangle\) is some independent-particle wave function (for instance, the Hartree-Fock determinant). Indices \(i\) and \(a\) correspond to virtual and occupied orbitals with respect to \(\vert \Phi_0 \rangle\), \(P\) and \(K\) denote the number of electron pairs (\(P = N/2\) with \(N\) being the total number of electrons) and orbitals, respectively. The geminal coefficient matrix (\(\mathbf{C}\)) of pCCD links the geminals with the underlying one-particle basis functions and has the following form,

(10.2)¶\[\begin{split}\mathbf{C} = \begin{pmatrix} 1 & 0 & \cdots & 0 & c_{1;P+1} & c_{1;P+2}&\cdots &c_{1;K}\\ 0 & 1 & \cdots & 0 & c_{2;P+1} & c_{2;P+2}&\cdots &c_{2;K}\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots &\ddots &\vdots\\ 0 & 0 & \cdots & 1 & c_{P;P+1} & c_{P;P+2}&\cdots & c_{P;K} \end{pmatrix}\end{split}\]

The exponential form of eq. (10.1) assures the size extensivity of the geminal wave function, however, to ensure the size consistency, one has to optimize the orbitals (see [boguslawski2014a] and [boguslawski2014b]). The simplest and most robust way is to use the variational orbital optimization (vOO-pCCD) method implemented in PyBEST (see pCCD with orbital optimization).

10.2. Currently supported features¶

If not mentioned otherwise, the pCCD module supports spin-restricted orbitals and the DenseLinalgFactory and CholeskyLinalgFactory. Specifically, the following features are provided:

Optimization of pCCD (eq. (10.1)) with (see pCCD with orbital optimization and [boguslawski2014a]) and without orbital optimization (see pCCD without orbital optimization) for a given Hamiltonian (see Getting started).

Variational orbital optimization and PS2c orbital optimization (see Summary of keyword arguments to choose the orbital optimizer and [boguslawski2014b])

Calculation of response one- and two-particle reduced density matrices (see Response density matrices)

Determination of pCCD natural orbitals and occupation numbers (see Response density matrices and Natural orbitals and occupation numbers)

Calculation of the exact orbital Hessian (see The exact orbital Hessian). Note that the orbital optimizer uses only a diagonal approximation to the orbital Hessian. Since the exact orbital Hessian is internally evaluated for the DenseLinalgFactory, PyBEST might require large amounts of memory to store the matrix representation of the orbital Hessian. Thus, it is recommended to evaluate the orbital Hessian only for small molecules or basis sets.

The pCCD wave function and its response density matrices can then be used for post-processing.

This version of PyBEST offers the following post-pCCD features:

A posteriori addition of dynamic electron correlation using perturbation theory (see The perturbation theory module and [boguslawski2017a]) or coupled cluster corrections (see The linearized coupled cluster module and [boguslawski2015b])

Analysis of orbital entanglement and correlations in the pCCD wave function using the orbital entanglement module (see Seniority zero wave functions)

10.3. Quick Guide: pCCD and orbital-optimized pCCD¶

If you use this part of the code, please cite [boguslawski2014a] and [boguslawski2014b].

We assume that you have performed a restricted Hartree-Fock calculation, whose results are stored in the IOData container hf_ (see The Self-Consistent Field Module). Furthermore, we will assume the following names for all PyBEST objects

lf

A LinalgFactory instance (see Preliminaries).

occ_model

An Aufbau occupation model of the AufbauOccModel class

kin

the kinetic energy integrals

ne

the nucleus-electron attraction integrals

eri

the two-electron repulsion integrals

10.3.1. How to: pCCD¶

The code snippet below shows how to perform a restricted pCCD calculation in PyBEST,

pccd = RpCCD(lf, occ_model)
pccd_ = pccd(kin, ne, eri, hf_)

The results are returned as an IOData container, while all results are saved to the pybest-results/checkpoint_pccd.h5 file. Specifically, the IOData container contains the following attributes

e_tot

The total pCCD energy (including all external terms)

e_corr

The pCCD correlation energy

e_ref

The energy of the reference determinant (here, the Hartree-Fock determinant)

t_p

The pair amplitudes

The pCCD module features additional options, which are discussed in greater detail below.

10.3.2. How to: orbital-optimized pCCD¶

The code snippet below shows how to perform an orbital-optimized pCCD calculation in PyBEST,

pccd = ROOpCCD(lf, occ_model)
pccd_ = pccd(kin, ne, eri, hf_)

The results are returned as an IOData container, while all results are saved to the pybest-resutls/checkpoint_pccd.h5 file. Similar to above, Similar to pCCD (without orbital optimization), the IOData container contains the following attributes

e_tot

The total pCCD energy (including all external terms)

e_corr

The pCCD correlation energy

e_ref

The energy of the reference determinant (here, the Hartree-Fock determinant)

t_p

The pair amplitudes

If the orbitals are optimized, the IOData container also includes the \(\lambda_p\) (electron-pair de-excitation) amplitudes, the optimized pCCD natural orbitals, and the response 1- and 2-particle reduced density matrices (RDMs).

orb_a

The pCCD natural orbitals

l_p

The pair de-excitation amplitudes

dm_1

The response 1-RDM

dm_2

The response 2-RDM

The OOpCCD module features additional options, which are discussed in greater detail below.

10.4. Detailed input structure of pCCD¶

10.4.1. Getting started¶

To optimize a pCCD wave function, the module requires some one- and two-electron integrals (that is, some Hamiltonian) and an initial guess for the orbitals as input arguments. PyBEST provides different options for specifying the Hamiltonian and an orbital guess.

Note

See Defining Basis Sets, Molecular Geometries, and Hamiltonians for how to define the molecular geometry, basis set, or the Hamiltonian.

The Hamiltonian is divided into three contributions: the one- and two-electron integrals as well as an external term (also referred to as core energy). Possible choices are (see below for examples):
1. In-house calculation of the quantum chemical Hamiltonian expressed in the AO basis (all terms are calculated separately in PyBEST; see Computing the matrix representation of the Hamiltonian for documentation).
2. In-house calculation of model Hamiltonians. Supported model Hamiltonians are summarized in Model Hamiltonians.
3. External (one- and two-electron) integrals (in an orthonormal basis) and core energy can be read from a file. The integral file must use the Molpro file format (see Dumping/Loading a Hamiltonian to/from a file for more details).
A set of initial guess orbitals can be either generated in PyBEST (including the AO overlap matrix) or read from disk. Examples for initial guess orbitals are:
1. Restricted canonical Hartree-Fock orbitals (see The Self-Consistent Field Module)
2. Localized orbitals (see Localization of molecular orbitals for documentation)

10.4.2. pCCD with orbital optimization¶

If you use this part of the code, please cite [boguslawski2014a] and [boguslawski2014b].

10.4.2.1. How to set up a calculation¶

After having specified a Hamiltonian and initial guess orbitals, a pCCD calculation with variational orbital optimization can be initiated as follows

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, hf_)

All arguments above have been introduced in Quick Guide: pCCD and orbital-optimized pCCD.

Note

The core energy (also referred to as external energy) is automatically determined from the IOData container hf_ and hence does not need to be passed explicitly.

The final results of the calculation are returned as an IOData container (here oopccd_)

Note

The number of electron pairs is automatically determined using the OccModel instance occ_model.

A set of keyword arguments is optional and contain optimization-specific options (like the number of orbital optimization steps, etc.). Their default values are chosen to give a reasonable performance and can be adjusted if convergence difficulties are encountered (see also Summary of keyword arguments).

After the pCCD calculation is finished (because pCCD converged or the maximum number of iterations was reached), all final results are, by default, stored in a checkpoint file checkpoint_pccd.h5 of the pybest-results directory and can be used for a subsequent restart.

10.4.2.2. Defining a frozen core¶

The pCCD module in PyBEST supports frozen core orbitals. These orbitals are not optimized during the orbital optimization, neither are they included in the amplitude equations. The frozen (core) orbitals are thus by construction doubly occupied. The number of frozen orbitals can be defined when creating an instance of the ROOpCCD class,

###############################################################################
## Do OO-pCCD optimization ####################################################
## freeze the 1s orbitals of the oxygen atom ##################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model, ncore=1)
oopccd_ = oopccd(kin, ne, er, hf_)

Note

If ncore is set to some nonzero number, the orbitals are frozen with respect to their order in the orb_a attribute passed to ROOpCCD. In general, these are the energetically lowest-lying orbitals or those with the largest occupation numbers.

The number of virtual orbitals cannot be changed and hence PyBEST assumes all virtual orbitals to be active.

10.4.2.3. How to restart¶

10.4.2.3.1. PyBEST checkpoint files¶

To restart a pCCD calculation (for instance, using the orbitals from a different molecular geometry as initial guess or from a previous calculation using the same molecular geometry), you can use the restart keyword in the function call providing the filename (or its full path)

###############################################################################
## Restart an OO-pCCD calculation from default restart file ###################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart="pybest-results/checkpoint_pccd.h5")

pybest-results/checkpoint_pccd.h5 is the default checkpoint file generated py PyBEST if not specified otherwise by the user.

10.4.2.3.2. Perturbed restart orbitals¶

Sometimes, it might be advantageous to slightly perturb some of the orbitals used for restarts (for instance, to overcome local minima or to facilitate orbital localization). This can be achieved similar to the SCF module as explained in Restarting from perturbed orbitals.

###############################################################################
## Restart an OO-pCCD calculation from perturbed orbitals #####################
###############################################################################
# read in some orbitals
restart = IOData.from_file("pybest-results/checkpoint_pccd.h5")
# perturb them as you wish, eg, HOMO-LUMO 2x2 rotation (python indexing, starts
# with 0)
restart.orb_a.rotate_2orbitals(60 * deg, 4, 5)

# Now pass the restart container as argument
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart)

After some restart file pybest-results/checkpoint_pccd.h5 has been read in using the from_file() method of the IOData container class, you can perturb the orbitals using the swapping feature (see Swapping orbitals) or 2x2 rotations (see Rotating orbitals). The modified/updated/perturbed IOData container can then be simply passed as an argument to the ROOpCCD function call (as done above with the RHF container).

10.4.2.4. Response density matrices¶

PyBEST supports the calculation of the response 1- and 2-particle reduced density matrices (1-RDM and 2-RDM), \(\gamma_{pq}\) and \(\Gamma_{pqrs}\), respectively. Since pCCD is a product of natural geminals, the 1-RDM is diagonal and is calculated from

\[\gamma_p = \langle \Psi_0| (1+\hat{\Lambda}) a^\dagger_p a_p | \textrm{pCCD} \rangle,\]

where \(\hat{\Lambda}\) contains the de-excitation operator,

\[\hat{\Lambda} = \sum_{ia} \lambda_i^a (a^\dagger_i a^\dagger_{\bar{i}} a_{\bar{a}} a_a - c_i^a).\]

The most recent response 1-RDM (a OneIndex instance) of OO-pCCD is saved in the return value of pCCD (that is the IOData container). The 1-RDM is stored as the dm_1 attribute of the IOData container and can be accessed as follows

one_dm = oopccd_.dm_1

The response 2-RDM is defined as

\[\Gamma_{pqrs} = \langle \Psi_0| (1+\hat{\Lambda})a^\dagger_p a^\dagger_{q} a_{s} a_r| \textrm{pCCD} \rangle.\]

In PyBEST, only the non-zero elements of the response 2-RDM are calculated, which are \(\Gamma_{pqpq}=\Gamma_{p\bar{q}p\bar{q}}\) and \(\Gamma_{p\bar{p}q\bar{q}}\). Specifically, the non-zero elements \(\Gamma_{pqpq}\) and \(\Gamma_{ppqq}\) (where we have omitted the information about electron spin) are calculated separately and stored as TwoIndex objects. Note that \(\gamma_p=\Gamma_{p\bar{p}p\bar{p}}\). Similar to the 1-RDM, the most recent 2-RDM is stored in the IOData container as the attribute dm_2. Since most of elements of the 2-RDM are zero, only its non-zero blocks are calculated and stored in PyBEST. To distinguish between the \(\Gamma_{pqpq}=\Gamma_{p\bar{q}p\bar{q}}\) and \(\Gamma_{p\bar{p}q\bar{q}}\) block, the 2-RDM is a dictionary with attributes ppqq and pqpq. They can be accessed as follows

two_dm = oopccd_.dm_2
# Access the ppqq block
two_dm_ppqq = two_dm['ppqq']
# Access the pqpq block
two_dm_pqpq = two_dm['pqpq']

Note

In PyBEST, we have \(\Gamma_{p\bar{p}q\bar{q}} = 0 \, \forall \, p=q \in \textrm{occupied}\) and \(\Gamma_{p\bar{q}p\bar{q}} = 0 \, \forall \, p=q \in \textrm{virtual}\).

10.4.2.5. Natural orbitals and occupation numbers¶

If pCCD converges, the final orbitals are the pCCD natural orbitals and are stored in the orb_a attribute of the IOData container. The natural orbitals can be exported to the molden file format (see Generating molden files) and visualized using, for instance, Jmol or VESTA.

The natural occupation numbers, the eigenvalues of the response 1-RDM (see Response density matrices) are stored in the occupations attribute (a 1-dim np.array) of orb_a and can be directly accessed after a pCCD calculation using

oopccd_.orb_a.occupations

10.4.2.6. The exact orbital Hessian¶

Although the orbital optimizer uses a diagonal approximation to the exact orbital Hessian, the exact orbital Hessian can be evaluated after an pCCD calculation. Note that this feature uses DenseLinalgFactory. If CholeskyLinalgFactory is chosen, the corresponding 2-electron integrals are transformed back to a dense object prior the calculation of the exact orbital Hessian. Thus, this feature is limited by the memory bottleneck of the 4-index transformation of DenseLinalgFactory (see also indextrans in Summary of keyword arguments). To calculate the exact orbital Hessian, all one-electron and two-electron integrals (two) need to be transformed into the pCCD MO basis first (see also The AO/MO Transformation of the Hamiltonian),

# transform integrals for restricted orbitals oopccd_.orb_a
ti_ = transform_integrals(kin, ne, er, oopccd_)

# transformed one-electron integrals: attribute 'one' (list)
(one,) = ti_.one   # or: one = ti_.one[0]
# transformed two-electron integrals: attribute 'two' (list)
(two,) = ti_.two   # or: two = ti_.two[0]

where kin and na (er) are the one- (two-)electron integrals expressed in the AO basis. In the above example, the molecular orbitals (here, the optimized pCCD orbitals) are passed as an IOData container oopccd_. PyBEST automatically assigns all arguments. Thus their order does not matter.

This step can be skipped if the one- and two-electron integrals are already expressed in the (optimized) MO basis. The transformed one- and two-electron integrals (first element in each list) are passed as function arguments to the get_exact_hessian() method of ROOpCCD, which returns a 2-dim np.array with elements \(H_{pq,rs} = H_{p,q,r,s}\),

# calculate the exact orbital hessian of OOpCCD
hessian = oopccd.get_exact_hessian(one, two)

where oopccd is an instance of ROOpCCD (see above pCCD with orbital optimization). The exact orbital Hessian can be diagonalized using, for instance, the np.linalg.eigvalsh routine,

# diagonalize the exact orbital Hessian
eigv = np.linalg.eigvalsh(hessian)

If frozen core orbitals have been defined for the pCCD optimization, the one- and two-electron integrals need to be transformed using the split_core_active() function (see also Defining an Active Space Hamiltonian for more details)

# Define an active space with ncore=1 frozen core orbitals and all virtual
# active orbitals from a pCCD wave function stored in oopccd_
cas = split_core_active(kin, ne, er, oopccd_, ncore=1)

# all one-electron integrals
one = cas.one
# all two-electron integrals
two = cas.two
# the core energy (not required for the orbital Hessian)
e_core = cas.e_core

The exact orbital Hessian can then be calculated and diagonalized using the transformed one- and two-electron integrals.

10.4.2.7. Summary of keyword arguments¶

The ROOpCCD module supports various keyword arguments that allow us to steer the optimization of the pCCD wave function, including orbital optimization. In the following, all supported keyword arguments are listed together with their default values. Please note, that for most cases the default values should be sufficient to reach convergence. For convergence difficulties see Troubleshooting in OOpCCD calculations.

restart

(str) if provided, a restart from a given checkpoint file is performed (see How to restart)

indextrans

(str) 4-index Transformation. The choice between tensordot (default) and einsum. tensordot is faster than einsum, but requires more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively.

warning

(boolean) if True, (scipy) solver-specific warnings are printed (default False)

guess

(dictionary) initial guess specifications:

type

(str) guess type. One of mp2 (default), random (random numbers), const (1.0 scaled by factor)

factor

(float) a scaling factor for the initial guess of type type (default -0.1, applies only to random and const)

geminal

(1-dim np.array) external guess for geminal coefficients (default None). If provided, type and factor are ignored. The elements of the geminal matrix of eq. (10.2) have to be indexed in C-like order. Note that the identity block is not required. The size of the 1-dim np.array is thus equal to the number of unknowns, that is, \(n_{\rm pairs}*n_{\rm virtuals}\).

lagrange

(1-dim np.array) external guess for Lagrange multipliers (default None). If provided, type and factor are ignored. The elements have to be indexed in C-like order. The size of the 1-dim np.array is equal to the number of unknowns, that is, \(n_{\rm pairs}*n_{\rm virtuals}\).

solver

(dictionary) scipy wave function/Lagrange solver:

wfn

(str) wave function solver (default krylov)

lagrange

(str) Lagrange multiplier solver (default krylov)

Note

The exact Jacobian of wfn and lagrange is not supported. Thus, scipy solvers that need the exact Jacobian cannot be used. See scipy root-solvers for more details.

maxiter

(dictionary) a maximum number of iterations:

wfniter

(int) maximum number of iterations for the wfn/lagrange solver (default 200)

orbiter

(int) maximum number of orbital optimization steps (default 100)

thresh

(dictionary) optimization thresholds:

wfn

(float) optimization threshold for geminal coefficients and Lagrange multipliers (default 1e-12)

energy

(float) convergence threshold for energy (default 1e-8)

gradientnorm

(float) convergence threshold for norm of orbital gradient (default 1e-4)

gradientmax

(float) threshold for maximum absolute value of the orbital gradient (default 5e-5)

printoptions

(dictionary) print level:

geminal

(boolean) if True, geminal matrix is printed (default False). Note that the identity block is omitted.

ci

(float) threshold for CI coefficients (requires evaluation of a permanent). All coefficients (for a given excitation order) larger than ci are printed (default 0.01)

excitationlevel

(int) a number of excited pairs w.r.t. the reference determinant for which the wave function amplitudes are reconstructed (default 1). At most, the coefficients corresponding to hextuply excited Slater determinants w.r.t the reference determinant can be calculated.

Note

The reconstruction of the wave function amplitudes requires evaluating a permanent which is in general slow (in addition to the factorial number of determinants in the active space).

dumpci

(dictionary) dump Slater determinants and corresponding CI coefficients to file:

amplitudestofile

(boolean) write wave function amplitudes to file (default False)

amplitudesfilename

(str) filename (default pccd_amplitudes.dat)

stepsearch

(dictionary) optimizes an orbital rotation step:

method

(str) step search method used. One of trust-region (default), None, backtracking

optimizer

(str) optimizes step to a boundary of trust radius in trust-region. One of pcg (preconditioned conjugate gradient), dogleg (Powell’s single dogleg step), ddl (Powell’s double-dogleg step) (default ddl)

alpha

(float) scaling factor for the Newton step. Used in backtracking and None method (default 1.00)

c1

(float) a parameter used in the Armijo condition of backtracking (default 1e-4)

minalpha

(float) minimum step length used in backracking (default 1e-6). If step length falls below minalpha, the backtracking line search is terminated and the most recent step is accepted

maxiterouter

(int) maximum number of iterations to optimize orbital rotation step (default 10)

maxiterinner

(int) maximum number of optimization steps in each step search (used only in pcg, default 500)

maxeta

(float) upper bound for estimated vs. actual change in trust-region (default 0.75)

mineta

(float) lower bound for estimated vs. actual change in trust-region (default 0.25)

upscale

(float) scaling factor to increase trust radius in trust-region (default 2.0)

downscale

(float) scaling factor to decrease trust radius in trust-region (default 0.25)

trustradius

(float) initial trust radius (default 0.75)

maxtrustradius

(float) maximum trust radius (default 0.75)

threshold

(float) trust-region optimization threshold, only used in pcg (default 1e-8)

checkpoint

(int) frequency of checkpointing. If checkpoint > 0, a checkpoint file is written every int iteration (default 1)

checkpoint_fn

(str) name of checkpoint file (default checkpoint_pccd.h5)

levelshift

(float) a level shift of Hessian (default 1e-8). The absolute value of elements of the orbital Hessian smaller than levelshift are shifted by levelshift

absolute

(boolean), if True, the absolute value of the orbital Hessian is taken (default False)

sort

(boolean), if True, orbitals are sorted according to their natural occupation numbers. This requires re-solving for the wave function after each orbital optimization step. Works only if orbitaloptimizer is set to variational (default True)

orbitaloptimizer

(str) switch between variational orbital optimization (variational) and PS2c orbital optimization (ps2c) (default variational)

10.4.3. pCCD without orbital optimization¶

If you use this part of the code, please cite [boguslawski2014a] and [boguslawski2014b].

10.4.3.1. How to set-up a calculation¶

If you only want to optimize the pCCD amplitudes (that is, to skip the orbital optimization step), use the RpCCD class instead of the ROOpCCD class. The remaining syntax is equivalent to the orbital optimized pCCD module,

pccd =  RpCCD(lf, occ_model)
pccd_ = pccd(kin, ne, eri, hf_)

The function call again returns an IOData container.

Note

In the RpCCD module, the \(\Lambda\) equations are not solved. Thus, no response density matrices are available.

If response density matrices are required, the ROOpCCD class has to be used. To skip the orbital optimization procedure, the maxiter keyword has to be adjusted as follows (see also Summary of keyword arguments)

pccd =  ROOpCCD(lf, occ_model)
pccd_ = pccd(kin, ne, eri, hf_, maxiter={'orbiter': 0})

This will force PyBEST to solve the \(\Lambda\) equations and calculate by default the 1- and 2-RDMs.

10.4.3.2. Defining a frozen core¶

Similar to the orbital-optimized version, frozen core orbitals are also supported. These orbitals are not included in the amplitude equations and are thus doubly occupied. The number of frozen orbitals can be defined when creating an instance of the RpCCD class,

###############################################################################
## Do pCCD calculations #######################################################
## freeze the 1s orbitals #####################################################
###############################################################################
pccd = RpCCD(lf, occ_model, ncore=1)
pccd_ = pccd(kin, ne, er, hf_)

Note

Similar to the ncore feature of ROOpCCD, the orbitals are frozen with respect to their order in the orb_a attribute passed to RpCCD. In general, these are the energetically lowest-lying orbitals or those with the largest occupation numbers.

The number of virtual orbitals cannot be changed and hence PyBEST assumes all virtual orbitals to be active.

10.4.3.3. Summary of keyword arguments¶

Similar to the ROOpCCD module, RpCCD also supports various keyword arguments that allow us to steer the optimization of the pCCD wave function. In the following, all supported keyword arguments are listed together with their default values.

indextrans

(str) 4-index Transformation. The choice between tensordot (default) and einsum. tensordot is faster than einsum, but requires more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively.

warning

(boolean) if True, (scipy) solver-specific warnings are printed (default False)

guess

(dictionary) initial guess specifications:

type

(str) guess type. One of mp2 (default), random (random numbers), const (1.0 scaled by factor)

factor

(float) a scaling factor for the initial guess of type type (default -0.1, applies only to random and const)

geminal

(1-dim np.array) external guess for geminal coefficients (default None). If provided, type and factor are ignored. The elements of the geminal matrix of eq. (10.2) have to be indexed in C-like order. Note that the identity block is not required. The size of the 1-dim np.array is thus equal to the number of unknowns, that is, \(n_{\rm pairs}*n_{\rm virtuals}\).

solver

(dictionary) scipy wave function solver:

wfn

(str) wave function solver (default krylov)

Note

The exact Jacobian of wfn is not supported. Thus, scipy solvers that need the exact Jacobian cannot be used. See scipy root-solvers for more details.

maxiter

(dictionary) maximum number of iterations:

wfniter

(int) maximum number of iterations for the wfn solver (default 200)

thresh

(dictionary) optimization thresholds:

wfn

(float) optimization threshold for geminal coefficients and Lagrange multipliers (default 1e-12)

printoptions

(dictionary) print level:

geminal

(boolean) if True, geminal matrix is printed (default False). Note that the identity block is omitted.

ci

(float) threshold for CI coefficients (requires evaluation of a permanent). All coefficients (for a given excitation order) larger than ci are printed (default 0.01)

excitationlevel

(int) number of excited pairs w.r.t. the reference determinant for which the wave function amplitudes are reconstructed (default 1). At most, the coefficients corresponding to hextuply excited Slater determinants w.r.t the reference determinant can be calculated.

Note

The reconstruction of the wave function amplitudes requires evaluating a permanent which is in general slow (in addition to the factorial number of determinants in the active space).

dumpci

(dictionary) dump Slater determinants and corresponding CI coefficients to file:

amplitudestofile

(boolean) write wave function amplitudes to file (default False)

amplitudesfilename

(str) filename (default pccd_amplitudes.dat)

10.5. Troubleshooting in OOpCCD calculations¶

How to change the number of orbital optimization steps:

To increase the number of iterations in the orbital optimization, adjust the keyword maxiter (see Summary of keyword arguments):
```
maxiter={'orbiter': int}
```
where int is the desired number of iterations
The energy oscillates during orbital optimization:

The occupied-virtual separation breaks down and the reference determinant cannot be optimized. In some cases, fixing the reference determinant might accelerate convergence. However, the final solution might not be reasonable if the optimized geminal coefficient matrix contains elements that are significantly larger than 1.0 in absolute value. To fix the reference determinant, the sort keyword has to be set to False (see Summary of keyword arguments):
```
sort=False
```
The orbital optimization converges very, very slowly: Usually, the orbital optimization converges fast around the equilibrium. For stretched distances (in the vicinity of dissociation, etc.) convergence can be very slow, especially if the final solution results in symmetry-broken orbitals. In such cases, the diagonal approximation to the Hessian is not optimal. However, the current version of PyBEST does not support orbital optimization with the exact Hessian nor Hessian updates. You can either try to start with localized guess orbitals, or some perturbed orbital guess (2x2 rotations of, for instance, s and p orbitals), or simply wait until the optimizer converges.
How to scan a potential energy surface:

To accelerate convergence, restart from adjacent points on the potential energy surface (see How to restart). Using Hartree-Fock orbitals as initial guess might result in convergence difficulties and optimization problems of the reference determinant.

Note

Different guess orbitals might result in different natural orbitals optimized within the OOpCCD module. This is especially true for stretched bond distances, where (partial) localization might occur rather quickly. Please check the final pCCD natural orbitals by, for instance, dumping them to Molden files (see Generating molden files).
How to perturb the orbitals:

The initial guess orbitals can be perturbed using a sequence of Givens rotations (see also Rotating orbitals). Givens rotations between orbital pairs can be used if, for instance, the orbital optimizer converges to a saddle point.

10.6. Example Python scripts¶

Several complete examples can be found in the directory data/examples/pccd.

10.6.1. The water molecule¶

This is a basic example on how to perform an orbital-optimized pCCD calculation in PyBEST. This script performs an orbital-optimized pCCD calculation on the water molecule using the cc-pVDZ basis set and canonical RHF orbitals as initial orbitals.

Listing 10.1 data/examples/pccd/water_oopccd_cc-pvdz.py¶

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals.rpccd import ROOpCCD
from pybest.io.iodata import IOData
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF
from pybest.units import deg


###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
fn_xyz = context.get_fn("test/water_2.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
er = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_ = hf(kin, ne, er, external, olp, orb_a)

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, hf_)

###############################################################################
## Restart an OO-pCCD calculation from default restart file ###################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart="pybest-results/checkpoint_pccd.h5")

###############################################################################
## Restart an OO-pCCD calculation from perturbed orbitals #####################
###############################################################################
# read in some orbitals
restart = IOData.from_file("pybest-results/checkpoint_pccd.h5")
# perturb them as you wish, eg, HOMO-LUMO 2x2 rotation (python indexing, starts
# with 0)
restart.orb_a.rotate_2orbitals(60 * deg, 4, 5)

# Now pass the restart container as argument
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart)

10.6.2. The water molecule including restart for the same geometry¶

This is a basic example on how to perform an orbital-optimized pCCD calculation in PyBEST. This script performs an orbital-optimized pCCD calculation on the water molecule using the cc-pVDZ basis set and canonical RHF orbitals as initial orbitals. Afterwards, the calculation is restarted from the default checkpoint file generated by the previous iteration.

Listing 10.2 data/examples/pccd/water_oopccd_cc-pvdz.py¶

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals.rpccd import ROOpCCD
from pybest.io.iodata import IOData
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF
from pybest.units import deg


###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
fn_xyz = context.get_fn("test/water_2.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
er = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_ = hf(kin, ne, er, external, olp, orb_a)

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, hf_)

###############################################################################
## Restart an OO-pCCD calculation from default restart file ###################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart="pybest-results/checkpoint_pccd.h5")

###############################################################################
## Restart an OO-pCCD calculation from perturbed orbitals #####################
###############################################################################
# read in some orbitals
restart = IOData.from_file("pybest-results/checkpoint_pccd.h5")
# perturb them as you wish, eg, HOMO-LUMO 2x2 rotation (python indexing, starts
# with 0)
restart.orb_a.rotate_2orbitals(60 * deg, 4, 5)

# Now pass the restart container as argument
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, ne, er, restart)

The next example shows how to perform a simple restart from some checkpoint file without any RHF calculation.

Listing 10.3 data/examples/pccd/water_oopccd_restart_cc-pvdz.py¶

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
)
from pybest.geminals.rpccd import ROOpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel


# You need to execute `water_oopccd_cc-pvdz.py` first, before running
# this example script.

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
fn_xyz = context.get_fn("test/water.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
na = compute_nuclear(obasis)
er = compute_eri(obasis)

###############################################################################
## Restart an OO-pCCD calculation from default restart file ###################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, na, er, restart="pybest-results/checkpoint_pccd.h5")

Note

This will only work, if the restart is performed for the same molecular geometry.

10.6.3. The water molecule including restart for a different geometry¶

This is a basic example on how to perform an orbital-optimized pCCD calculation in PyBEST, where the orbitals from a different molecular geometry are taken as guess orbitals.

Listing 10.4 data/examples/pccd/water_oopccd_restart_perturbed_cc-pvdz.py¶

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals.rpccd import ROOpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.scf import AufbauOccModel


###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
fn_xyz = context.get_fn("test/water_2.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
na = compute_nuclear(obasis)
er = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Restart an OO-pCCD calculation from default restart file ###################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(
    kin, na, er, olp, orb_a, e_core=external, restart="pybest-results/checkpoint_pccd.h5"
)

Note

If the molecular structure is perturbed, we need to explicitly pass the external energy and the overlap matrix (as well as an empty set of orbitals). If not provided, the restart will result in wrong energies.

10.6.4. pCCD with external integrals (FCIDUMP)¶

This is a basic example on how to perform an orbital-optimized pCCD calculation using one- and two-electron integrals from an external file. The number of doubly-occupied orbitals is 5 and has to be specified, while the total number of basis functions is 28 (it is set automatically).

Listing 10.5 data/examples/pccd/water_oopccd_fcidump.py¶

from pybest import context
from pybest.geminals.rpccd import ROOpCCD
from pybest.io.iodata import IOData
from pybest.scf import AufbauOccModel


# Read Hamiltonian from some FCIDUMP file
# ---------------------------------------
fcidump = context.get_fn("test/water.FCIDUMP")
mol = IOData.from_file(fcidump)

# Define occupation model
# ------------------------
occ_model = AufbauOccModel(5)

# Do OOpCCD optimization
# ----------------------
oopccd = ROOpCCD(mol.lf, occ_model)
oopccd_ = oopccd(mol.one, mol.two, mol)

10.6.5. pCCD using model Hamiltonians: Hubbard¶

This is a basic example on how to perform an orbital-optimized pCCD calculation using the half-filled 1-D Hubbard model Hamiltonian with 6 sites. The number of doubly-occupied sites is thus 3. The t parameter is set to -1, the U parameter is set to 1. Periodic boundary conditions are employed.

Listing 10.6 data/examples/pccd/hubbard_oopccd.py¶

from pybest.geminals.rpccd import ROOpCCD
from pybest.linalg.dense import DenseLinalgFactory
from pybest.modelhamiltonians.physmodham import Hubbard
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF


###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(6)
occ_model = AufbauOccModel(3)
orb_a = lf.create_orbital()
###############################################################################
## Initialize an instance of the Hubbard class ################################
###############################################################################
modelham = Hubbard(lf, pbc=True)
olp = modelham.compute_overlap()
###############################################################################
# t-param, t = -1 #############################################################
###############################################################################
kin = modelham.compute_one_body(-1)
###############################################################################
# U-param, U = 2 ##############################################################
###############################################################################
two = modelham.compute_two_body(1)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_ = hf(kin, two, 0.0, orb_a, olp)

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
oopccd = ROOpCCD(lf, occ_model)
oopccd_ = oopccd(kin, two, hf_)