PNO-Based Orbital Refinement using FCI

In this tutorial, we will demonstrate how to perform the orbital refinement using pair natural orbitals (PNOs) as initial guess orbitals for a ground state FCI calculation. The tutorial is structured as follows:

1. Define the molecule and parameters for the calculation
2. Get the initial pair natural orbitals
3. Calculate the initial integrals
4. Perform the FCI calculation and extract the one-body and two-body reduced density matrices (rdms)
5. Perform the orbital refinement using the 1-body and 2-body rdms obtained from the FCI calculation
6. Repeat the orbital refinement until convergence is reached

molecule_name = "h4"
n_electrons = 4
iterations = 10
box_size = 50.0
wavelet_order = 7
madness_thresh = 0.0001
econv = 1.0e-6

1

Number of electrons

2

Maximal number of iterations for orbital refinement

3

Size of the simulation box

4

Order of wavelet basis functions

5

Threshold for numerical precision of function representation

6

Energy convergence threshold for orbital refinement

Get pair natural orbitals (PNOs)

To get an initial guess for the FCI calculation, we use the pair natural orbitals (PNOs) generated with the FrayedEnds MadPNO class. First, the geometry of the H4 molecule is defined, and the numerical environment for the multiwavelet representation is set up using the MadWorld3D class. The PNOs are then computed and orthonormalized. Additionally, the nuclear repulsion energy and the nuclear potential are computed for later use in the FCI calculation.

import frayedends as fe

molecule = fe.MolecularGeometry(units="angstrom")
molecule.add_atom(0.0, 0.0, -1.5, "H")
molecule.add_atom(0.0, 0.0, -0.5, "H")
molecule.add_atom(0.0, 0.0, 0.5, "H")
molecule.add_atom(0.0, 0.0, 1.5, "H")
geom = "H 0.0 0.0 -1.5\nH 0.0 0.0 -0.5\nH 0.0 0.0 0.5\nH 0.0 0.0 1.5\n"

world = fe.MadWorld3D(L=box_size, k=wavelet_order, thresh=madness_thresh)

madpno = fe.MadPNO(world, geom, units="angstrom", n_orbitals=8)
orbs = madpno.get_orbitals()

integrals = fe.Integrals3D(world)
orbs = integrals.orthonormalize(orbitals=orbs)

nuc_repulsion = madpno.get_nuclear_repulsion()
Vnuc = madpno.get_nuclear_potential()

n_orbitals = len(orbs)

for i in range(len(orbs)):
    orbs[i].type = "active"

for i in range(len(orbs)):
    world.cube_plot(f"initial_orb{i}", orbs[i], molecule)

MADNESS runtime initialized with 9 threads in the pool and affinity OFF

1. Define a linear H4 molecule geometry with 1.0 Angstrom spacing between adjacent atoms
2. Setting up the numerical environment for the MRA calculations by creating a frayedends world object with the specified parameters
3. Compute the pair natural orbitals (PNOs) for the H4 molecule using the MadPNO class in FrayedEnds
4. Orthonormalize the PNOs using the orthonormalization method in the Integrals3D class
5. Compute the nuclear repulsion energy
6. Compute the nuclear potential for the H4 molecule

Orbital Refinement using FCI

In each iteration of the orbital refinement, the two-body, kinetic, and potential integrals are first computed using the current orbitals. With these integrals, the FCI calculation is performed to obtain the electronic energy and the FCI vector. From the FCI vector, the 1-body and 2-body reduced density matrices (rdms) are constructed and used for the orbital refinement. The orbital refinement loop continues until the energy converges within the specified threshold.

import numpy as np
from pyscf import fci

current = 0.0
for iteration in range(iterations):
    integrals = fe.Integrals3D(world)
    G = integrals.compute_two_body_integrals(orbs, ordering="chem").elems
    T = integrals.compute_kinetic_integrals(orbs)
    V = integrals.compute_potential_integrals(orbs, Vnuc)

    e, fcivec = fci.direct_spin0.kernel(T + V, G, n_orbitals, n_electrons)
    rdm1, rdm2 = fci.direct_spin0.make_rdm12(
        fcivec,
        n_orbitals,
        n_electrons,
    )
    rdm2 = np.swapaxes(rdm2, 1, 2)

    e_tot = e + nuc_repulsion

    print("iteration {} FCI electronic energy {:+2.8f}, total energy {:+2.8f}".format(iteration, e, e_tot))

    opti = fe.Optimization3D(world, Vnuc, nuc_repulsion)
    orbs = opti.get_orbitals(
        orbitals=orbs,
        rdm1=rdm1,
        rdm2=rdm2,
        opt_thresh=0.001,
        occ_thresh=0.001,
    )

    for i in range(n_orbitals):
        world.cube_plot(f"orb{i}", orbs[i], molecule)

    if np.isclose(e_tot, current, atol=econv, rtol=0.0):
        break
    current = e_tot

fe.cleanup(globals())

iteration 0 FCI electronic energy -4.53404243, total energy -2.24094119
iteration 1 FCI electronic energy -4.53661359, total energy -2.24351234
iteration 2 FCI electronic energy -4.53672886, total energy -2.24362762
iteration 3 FCI electronic energy -4.53673548, total energy -2.24363423
iteration 4 FCI electronic energy -4.53673586, total energy -2.24363461

1. Create an integrals object to compute the initial integrals
2. Compute the two-body, kinetic and potential integrals using the orbitals obtained from NWChem
3. Perform the FCI calculation using the integrals
4. Compute the 1-body and 2-body reduced density matrices (rdms) from the FCI vector
5. Compute the total energy by adding the nuclear repulsion energy to the electronic energy obtained from the FCI calculation
6. Perform the orbital refinement using the 1-body and 2-body rdms obtained from the FCI calculation
7. Loop terminates as soon as the energy changes less than econv in one iteration step

The orbital refinement process is repeated until convergence is reached. At the end of the orbital refinement, the ground state energy should be improved compared to the initial FCI calculation with the PNO orbitals as guess orbitals.

The initial PNOs and the refined orbitals after the orbital optimization can be visualized using the cube files and plotted interactively using py3Dmol.

import py3Dmol

initial_orb = open("initial_orb2.cube").read()
refined_orb = open("orb2.cube").read()

v = py3Dmol.view(width=800, height=400, viewergrid=(1, 2))

v.addVolumetricData(initial_orb, "cube", {"isoval": -0.001, "color": "red", "opacity": 0.75}, viewer=(0, 0))
v.addVolumetricData(initial_orb, "cube", {"isoval": 0.001, "color": "blue", "opacity": 0.75}, viewer=(0, 0))
v.addModel(initial_orb, "cube", viewer=(0, 0))
v.setStyle({"sphere": {}}, viewer=(0, 0))

v.addVolumetricData(refined_orb, "cube", {"isoval": -0.001, "color": "red", "opacity": 0.75}, viewer=(0, 1))
v.addVolumetricData(refined_orb, "cube", {"isoval": 0.001, "color": "blue", "opacity": 0.75}, viewer=(0, 1))
v.addModel(refined_orb, "cube", viewer=(0, 1))
v.setStyle({"sphere": {}}, viewer=(0, 1))

v.zoomTo()
v.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.