4.2. The AP1roG module¶

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

4.2.1. The AP1roG model¶

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

(1)\[\vert \textrm{AP1roG}\rangle = \exp(\sum_{ia} c_i^a a^\dagger_a a^\dagger_{\bar{a}} a_{\bar{i}} a_i) \vert \Psi_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 \Psi_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 \Psi_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 AP1roG links the geminals with the underlying one-particle basis functions and has the following form,

(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. (1) assures the size extensivity of the geminal wavefunction, however, in order 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-AP1roG) method implemented in HORTON (see AP1roG with orbital optimization).

4.2.2. Currently supported features¶

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

Optimization of AP1roG (eq. (1)) with (see AP1roG with orbital optimization) and without orbital optimization (see AP1roG 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)

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

Determination of AP1roG 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 Hessian. The exact orbital Hessian can only be evaluated in combination with the DenseLinalgFactory.

The AP1roG wave function and its response density matrices can then be used for post-processing. This version of HORTON offers:

A posteriori addition of dynamic electron correlation using the perturbation module (see The PTa module and The PTb module for documentation)

Analysis of orbital correlations in the AP1roG wave function using the orbital entanglement module (see Orbital entanglement and correlation for a seniority zero wavefunction for documentation)

Dump the Hamiltonian (collection of one- and two-electron integrals and the energy term due to core electrons and external potentials) in the AP1roG MO basis. The one- and two-electron integrals can be calculated for any pre-defined active space, that is, a selected number of active electrons and orbitals. The Hamiltonian is stored in the Molpro file format (see Dumping/Loading a Hamiltonian to/from a file for documentation)

4.2.3. Input structure¶

4.2.3.1. Getting started¶

To optimize an AP1roG wavefunction, the module requires a Hamiltonian and an initial guess for the orbitals (either an AO/MO coefficient matrix or an MO/MO coefficient matrix) as input arguments. HORTON provides different options for specifying the Hamiltonian and an orbital guess.

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:
1. In-house calculation of the quantum chemical Hamiltonian expressed in the AO basis (kinetic energy of the electrons, electron-nuclear attraction, electron-electron repulsion, and nuclear-nuclear repulsion). All terms are calculated separately in HORTON (see Computing (Hamiltonian) matrix elements for documentation). Note, however, that all one-electron terms have to be combined into one single operator term. This can be done in the following way
  ############################################################################### ## Calculate kinetic energy (kin) and nuclear attraction (na) term ############ ############################################################################### kin = obasis.compute_kinetic(lf) na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf) ############################################################################### ## Combine one-electron integrals to single Hamiltonian ####################### ############################################################################### one = kin.copy() one.iadd(na)
2. In-house calculation of model Hamiltonians. Supported model Hamiltonians are summarized in Defining a model Hamiltonian. If the model Hamiltonian contains separate one-electron contributions, they have to be combined to a single operator as shown under point 1.
3. External (one- and two-electron) integrals (in an orthonormal basis) and core energy can be read from file. The integral file must use the Molpro file format (see Dumping/Loading a Hamiltonian to/from a file for more details). To load a Hamiltonian from file, run
  one, two, coreenergy = load_fcidump(lf, filename='./FCIDUMP')
with arguments

lf: A linear algebra factory. Must be of type DenseLinalgFactory. Note that CholeskyLinalgFactory is not supported

and optional arguments

filename: (str) the filename of the fcidump file (default FCIDUMP)

The function load_fcidump has three return values; the one-electron integrals (one) stored as a TwoIndex object, the two-electron integrals (two) stored as a FourIndex object, and the core energy (coreenergy, float).
A set of initial guess orbitals can be either generated in HORTON (including the AO overlap matrix) or read from disk (see How to restart to use orbitals generated in HORTON as initial guess). Examples for initial guess orbitals are:
1. Restricted canonical Hartree-Fock orbitals (see How to use HORTON as a Hartree-Fock/DFT program)
2. Localized orbitals. HORTON supports Pipek-Mezey localization of canonical Hartree-Fock orbitals. See Localization of molecular orbitals for documentation.
3. If external integrals (expressed in an orthonormal basis) are used to define the Hamiltonian, the initial orbitals and the overlap matrix are the identity matrix and can be set as follows:
  orb = lf.create_expansion(nbasis) olp = lf.create_two_index(nbasis) olp.assign_diagonal(1.0) orb.assign(olp)
where nbasis is the number of basis function (total number of orbitals in the active space).

4.2.3.2. AP1roG with orbital optimization¶

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

4.2.3.2.1. How to set-up a calculation¶

After specifying a Hamiltonian and initial guess orbitals, you can create an instance of the RAp1rog class,

apr1og =  RAp1rog(lf, occ_model, npairs=None, nvirt=None)

with arguments

lf: A LinalgFactory instance

occ_model: (AufbauOccModel instance) an Aufbau occupation model

and optional arguments

npairs: (int) number of electron pairs. If not specified, the number of pairs equals the number of occupied orbitals in the AufbauOccModel

nvirt: (int) number of virtual orbitals. If not specified, the number of virtual orbitals is calculated from the total number of basis functions minus the number of electron pairs.

Note that no optional arguments need to be specified for the AP1roG model and the number of electron pairs and virtual orbitals is automatically determined using the AufbauOccModel. A restricted, orbital-optimized AP1roG calculation can be initiated with a function call,

energy, c, l = ap1rog(one, two, external, orb, olp, scf, **keywords)

with arguments

one: (TwoIndex instance) the one-electron integrals

two: (FourIndex or CholeskyLinalgFactory instance) the two-electron integrals

external: (float) energy contribution due to an external potential, e.g., nuclear-nuclear repulsion term, etc.

orb: (Expansion instance) the AO/MO or MO/MO coefficient matrix. It also contains information about orbital energies (not defined in the AP1roG model) and occupation numbers

olp: (TwoIndex instance) the AO overlap matrix or, in case of an orthonormal basis, the identity matrix

scf: (boolean) if True, orbitals are optimized

The keyword arguments are optional and contain optimization-specific options (like the number of orbital optimization steps, etc.) as well as orbital manipulation schemes (Givens rotations, swapping orbitals in reference determinant, etc.). Their default values are chosen to give reasonable performance and can be adjusted if convergence difficulties are encountered. All keyword arguments are summarized in the following section (Summary of keyword arguments).

The function call gives 3 return values,

energy: (float) the total AP1roG electronic energy (the external term included)

c: (TwoIndex instance) the geminal coefficient matrix (without the diagonal occupied sub-block, see The AP1roG model)

l: (TwoIndex instance) the Lagrange multipliers (can be used to calculated the response 1-RDM)

After the AP1roG calculation is finished (because AP1roG converged or the maximum number of iterations was reached), the orbitals and the overlap matrix are, by default, stored in a checkpoint file checkpoint.h5 and can be used for a subsequent restart. Note that the geminal coefficient matrix and Lagrange multipliers are not stored after the calculation is completed.

4.2.3.2.2. Summary of keyword arguments¶

indextrans:
(str) 4-index Transformation. 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 random (random numbers, default), const (1.0 scaled by factor)

factor: (float) a scaling factor for the initial guess of type type (default -0.1)

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. (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 wavefunction/Lagrange solver:

wfn: (str) wavefunction solver (default krylov)

lagrange: (str) Lagrange multiplier solver (default krylov)

Note that 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) 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 True). 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 wavefunction 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 that the reconstruction of the wavefunction 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 wavefunction amplitudes to file (default False)

amplitudesfilename:

(str) file name (default ap1rog_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 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 Newton step. Used in backtracking and None method (default 1.00)

c1: (float) 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, orbitals (orb.hdf5) and overlap (olp.hdf5) are written to disk (default 1)

levelshift:
(float) level shift of Hessian (default 1e-8). 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 wavefunction after each orbital optimization step. Works only if orbitaloptimizer is set to variational (default True)

swapa:
(2-dim np.array) swap orbitals. Each row in swapa contains 2 orbital indices to be swapped (default np.array([[]]))

givensrot:
(2-dim np.array) rotate two orbitals using a Givens rotation. Each row in givensrot contains 2 orbital indices and the rotation angle in deg (default np.array([[]])). Orbitals are rotated sequentially according to the rows in givensrot. If a sequence of Givens rotation is performed, note that the indices in givensrot refer to the already rotated orbital basis. If givensrot is combined with swapa, all orbitals are swapped prior to any Givens rotation

orbitaloptimizer:

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

4.2.3.2.3. How to restart¶

To restart an AP1roG calculation (for instance, using the orbitals from a different molecular geometry as initial guess or from a previous calculation using the same molecular geometry), the molecular orbitals (of the previous wavefunction run) need to be read from disk,

old = IOData.from_file('checkpoint.h5')

If this checkpoint file is used as an initial guess for a new geometry, the orbitals must be re-orthogonalized w.r.t. the new basis:

# It is assume that the following two variables available:
#   olp: the overlap matrix of the new basis (geometry)
#   orb: an instance of DenseExpansion to which the new orbitals will be
#        written
project_orbitals_ortho(old.olp, olp, old.exp_alpha, orb)

Then, generate an instance of the RAp1rog class and perform a function call:

apr1og =  RAp1rog(lf, occ_model)
energy, c, l = ap1rog(one, two, external, orb, olp, True)

Note that all optional arguments have been omitted and orb was set to True.

4.2.3.2.4. Response density matrices¶

HORTON 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 AP1roG 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{AP1roG} \rangle,\]

where \(\hat{\Lambda}\) contains the deexcitation operator,

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

The response 1-RDM (a OneIndex instance) can be calculated in HORTON as follows

one_dm = lf.create_one_index()
ap1rog.compute_1dm(one_dm, c, l, factor=2.0, response=True)

where ap1rog is an instance of the RAp1rog class and (see also How to set-up a calculation to get c and l)

one_dm:	(`OneIndex` instance) output argument that contains the response 1-RDM in the _array attribute
c:	(`TwoIndex` instance) the geminal coefficient matrix (without the diagonal occupied sub-block, see The AP1roG model)
l:	(`TwoIndex` instance) the Lagrange multipliers

and optional arguments

factor:	(float) a scaling factor for the 1-RDM. If `factor=2.0`, the spin-summed 1-RDM is calculated, as \(\gamma_{p}=\gamma_{\bar{p}}\) (default `1.0`)
response:	(boolean) if `True`, the response 1-RDM is calculated (default `True`)

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{AP1roG} \rangle.\]

In HORTON, 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}}\).

twoppqq = lf.create_two_index()
twopqpq = lf.create_two_index()
###############################################################################
## Gamma_ppqq #################################################################
###############################################################################
ap1rog.compute_2dm(twoppqq, one_dm, c, l, 'ppqq', response=True)
###############################################################################
## Gamma_pqpq #################################################################
###############################################################################
ap1rog.compute_2dm(twopqpq, one_dm, c, l, 'pqpq', response=True)

with arguments (see again How to set-up a calculation to get c and l)

twopqpq/twoppqq:

(TwoIndex instance) output argument that contains the response 2-RDM

one_dm: (OneIndex instance) the response 1-RDM

c: (TwoIndex instance) the geminal coefficient matrix (without the diagonal occupied sub-block, see The AP1roG model)

l: (TwoIndex instance) the Lagrange multipliers

and optional arguments

response: (boolean) if True, the response 1-RDM is calculated (default True)

Note that, in HORTON, \(\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}\).

4.2.3.2.5. Natural orbitals and occupation numbers¶

If AP1roG converges, the final orbitals are the AP1roG natural orbitals and are stored in orb (see AP1roG with orbital optimization how to obtain orb). The natural orbitals can be exported to the molden file format (see Data file formats (input and output)) and visualized using, for instance, Jmol or VESTA.

The natural occupation numbers, the eigenvalues of the response 1-RDM (see Response density matrices for how to calculate response RDMs) are stored in the occupations attribute (a 1-dim np.array) of orb and can be directly accessed after an AP1roG calculation using

orb.occupations

4.2.3.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 AP1roG calculation. Note that this feature is only available for the DenseLinalgFactory. The CholeskyLinalgFactory does not allow for 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, the one-electron (one) and two-electron integrals (two) need to be transformed into the AP1roG MO basis first,

onemo, twomo = transform_integrals(one, two, indextrans, orb)

with arguments

one: (TwoIndex instance) the one-electron integrals

two: (FourIndex instance) the two-electron integrals

indextrans: (str) the 4-index transformation, either tensordot (preferred) or einsum

orb: (Expansion instance) the AO/MO coefficient matrix

and return values

onemo: (list of TwoIndex instances, one element for each spin combination) the transformed one-electron integrals

twomo: (list of FourIndex instances, one element for each spin combination) the transformed two-electron integrals

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 attribute function of RAp1rog which returns a 2-dim np.array with elements \(H_{pq,rs} = H_{p,q,r,s}\),

hessian = ap1rog.get_exact_hessian(onemo[0], twomo[0])

where ap1rog is an instance of RAp1rog (see AP1roG with orbital optimization). The exact orbital Hessian can be diagonalized using, for instance, the np.linalg.eigvalsh routine,

eigv = np.linalg.eigvalsh(hessian)

4.2.3.3. AP1roG without orbital optimization¶

If you use this part of the module, please cite [limacher2013]

4.2.3.3.1. How to set-up a calculation¶

There are two different ways of performing a restricted AP1roG calculation without orbital optimization. Similar to the orbital-optimized counterpart, we have to create an instance of the RAp1rog class (using the same abbreviations for arguments as in AP1roG with orbital optimization and excluding all optional arguments):

apr1og =  RAp1rog(lf, occ_model)

The geminal coefficients can be optimized as follows

Use a function call and set orb=False:

energy, c = ap1rog(one, two, external, orb, olp, False, **keywords)
where we have used the same abbreviations for function arguments as in AP1roG with orbital optimization. The keyword arguments contain optimisation-specific options and are summarized in Summary of keyword arguments.

Note that the function call gives only 2 return values (the Lagrange multipliers are not calculated):

energy: (float) the total AP1roG electronic energy (the external term included)

c: (TwoIndex instance) the geminal coefficient matrix (without the diagonal occupied sub-block, see The AP1roG model)

In contrast to the orbital-optimized code, the orbitals and the overlap matrix are not stored to disk after the AP1roG calculation is finished.

Use the orbital-optimized version of AP1roG (see AP1roG with orbital optimization), but set orbiter to 0, which suppresses an orbital-rotation step:
energy, c, l = ap1rog(one, two, external, orb, olp, True, **{ 'maxiter': {'orbiter': 0} })
The function call gives 3 return values (total energy energy, geminal coefficients c, and Lagrange multipliers l). The Lagrange multipliers can be used for post-processing.

4.2.3.3.2. Summary of keyword arguments¶

indextrans:
(str) 4-index Transformation. Choice between tensordot (default) and einsum. tensordot is faster than einsum, requires, however, 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 containing:

type: (str) guess type. One of random (random numbers, default), const (constant numbers)

factor: (float) a scaling factor for the initial guess of type type (default -0.1)

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. (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:
wfn/Lagrange solver (dictionary) containing:

wfn: (str) wavefunction solver (default krylov)

Note that 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 containing:

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

thresh:
(dictionary) optimization thresholds containing:

wfn: (float) optimization threshold for geminal coefficients (default 1e-12)

printoptions:
print level; dictionary containing:

geminal: (boolean) if True, geminal matrix is printed (default True). 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 wavefunction 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 that the reconstruction of the wavefunction 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 wavefunction amplitudes to file (default False)

amplitudesfilename:

(str) file name (default ap1rog_amplitudes.dat)

swapa:
(2-dim np.array) swap orbitals. Each row in swapa contains 2 orbital indices to be swapped (default np.array([[]]))

4.2.4. Troubleshooting in AP1roG-SCF 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 HORTON does not support orbital optimization with the exact Hessian nor Hessian updates.
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.
How to perturb the orbitals:

The initial guess orbitals can be perturbed using a sequence of Givens rotations (see also Summary of keyword arguments),
```
'givensrot': np.array([[orbindex1a, orbindex1b, angle1],[orbindex2a, orbindex2b, angle2],...])
```
where orbitals with indices orbindex1a and orbindex1b are rotated by angle angle1, etc. Givens rotations between orbital pairs can be used if, for instance, the orbital optimizer converges to a saddle point.

4.2.5. Example Python scripts¶

Several complete examples can be found in the directory data/examples/ap1rog. Three of these are discussed in the following subsections.

4.2.5.1. The water molecule (a minimum input example)¶

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

data/examples/ap1rog/water_minimal.py¶

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from HORTON's test data directory
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
two = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(two, 'hartree'),
    RExchangeTerm(two, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, c, l = ap1rog(one, two, external['nn'], orb, olp, True)

4.2.5.2. The water molecule (with all default keyword arguments)¶

This is the same example as above, but all keyword arguments are mentioned explicitly using their default values.

data/examples/ap1rog/water_default.py¶

import numpy as np

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from HORTON's test data directory
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
two = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(two, 'hartree'),
    RExchangeTerm(two, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, c, l = ap1rog(one, two, external['nn'], orb, olp, True, **{
        'indextrans': 'tensordot',
        'warning': False,
        'checkpoint': 1,
        'levelshift': 1e-8,
        'absolute': False,
        'givensrot': np.array([[]]),
        'swapa': np.array([[]]),
        'sort': True,
        'guess': {'type': 'random', 'factor': -0.1, 'geminal': None, 'lagrange': None},
        'solver': {'wfn': 'krylov', 'lagrange': 'krylov'},
        'maxiter': {'wfniter': 200, 'orbiter': 100},
        'dumpci': {'amplitudestofile': False, 'amplitudesfilename': './ap1rog_amplitudes.dat'},
        'thresh': {'wfn':  1e-12, 'energy': 1e-8, 'gradientnorm': 1e-4, 'gradientmax': 5e-5},
        'printoptions': {'geminal': True, 'ci': 0.01, 'excitationlevel': 1},
        'stepsearch': {'method': 'trust-region', 'alpha': 1.0, 'c1': 0.0001, 'minalpha': 1e-6, 'maxiterouter': 10, 'maxiterinner': 500, 'maxeta': 0.75, 'mineta': 0.25, 'upscale': 2.0, 'downscale': 0.25, 'trustradius': 0.75, 'maxtrustradius': 0.75, 'threshold': 1e-8, 'optimizer': 'ddl'},
        'orbitaloptimizer': 'variational'
    })

4.2.5.3. AP1roG with external integrals¶

This is a basic example on how to perform an orbital-optimized AP1roG calculation using one- and two-electron integrals from an external file. The number of doubly-occupied orbitals is 5, while the total number of basis functions is 28. See Defining a model Hamiltonian.

data/examples/ap1rog/external_hamiltonian_n2_dense.py¶

from horton import *
from horton.test.common import numpy_seed

# Read Hamiltonian from file 'FCIDUMP'
# ------------------------------------
# The required FCIDUMP file can be generated with the script
# data/examples/hf_dft/rhf_n2_dense.py
mol = IOData.from_file('n2.FCIDUMP')

# Define Occupation model, expansion coefficients and overlap
# -----------------------------------------------------------
nocc = 5
occ_model = AufbauOccModel(nocc)
orb = mol.lf.create_expansion()
olp = mol.lf.create_two_index()
olp.assign_diagonal(1.0)
orb.assign(olp)

# Do OO-AP1roG optimization
# -------------------------
ap1rog = RAp1rog(mol.lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, c, l = ap1rog(mol.one_mo, mol.two_mo, mol.core_energy, orb, olp, True)

4.2.5.4. AP1roG using model Hamiltonians¶

This is a basic example on how to perform an orbital-optimized AP1roG calculation using 1-D Hubbard model Hamiltonian. The number of doubly-occupied sites is 3, the total number of sites is 6. The t parameter is set to -1, the U parameter is set to 2, and periodic boundary conditions are employed.

data/examples/ap1rog/hubbard.py¶

import numpy as np

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(6)
occ_model = AufbauOccModel(3)
modelham = Hubbard(pbc=True)
orb = lf.create_expansion(6)
olp = modelham.compute_overlap(lf)
###############################################################################
# t-param, t = -1 #############################################################
###############################################################################
kin = modelham.compute_kinetic(lf, -1)
###############################################################################
# U-param, U = 2 ##############################################################
###############################################################################
two = modelham.compute_er(lf, 2)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, orb)
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(two, 'hartree'),
    RExchangeTerm(two, 'x_hf'),
]
ham = REffHam(terms)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver()
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, c, l = ap1rog(kin, two, 0, orb, olp, True)

Note that for the Hubbard model, the external potential has to be set to 0,

energy, c, l = ap1rog(kin, two, 0, orb, olp, True)