Adds better examples

examples/00-ccsd.py

@@ -1,15 +1,22 @@
+"""
+Example of a simple CCSD calculation.
+"""
+
 import numpy as np
 from pyscf import gto, scf
 
 from ebcc import EBCC
 
+# Define the molecule using PySCF
 mol = gto.Mole()
 mol.atom = "H 0 0 0; F 0 0 1.1"
 mol.basis = "cc-pvdz"
 mol.build()
 
+# Run a Hartree-Fock calculation using PySCF
 mf = scf.RHF(mol)
 mf.kernel()
 
+# Run a CCSD calculation using EBCC
 ccsd = EBCC(mf)
 ccsd.kernel()

examples/01-ebccsd.py

@@ -1,15 +1,21 @@
+"""
+Example of some electron-boson CCSD calculations.
+"""
+
 import numpy as np
 from pyscf import gto, scf
 
 from ebcc import EBCC
 
 np.random.seed(123)
 
+# Define the molecule using PySCF
 mol = gto.Mole()
 mol.atom = "H 0 0 0; F 0 0 1.1"
 mol.basis = "cc-pvdz"
 mol.build()
 
+# Run a Hartree-Fock calculation using PySCF
 mf = scf.RHF(mol)
 mf.kernel()
 
@@ -26,13 +32,13 @@
 #    v v v v
 # ____ _ _ _
 # CCSD-S-1-1: One-boson amplitudes and one-boson-one-fermion coupling
-ccsd = EBCC(
+ccsd_s_1_1 = EBCC(
     mf,
     ansatz="CCSD-S-1-1",
     omega=omega,
     g=g,
 )
-ccsd.kernel()
+ccsd_s_1_1.kernel()
 
 #    ,--------- Fermionic ansatz
 #    |  ,------ Bosonic excitation amplitudes
@@ -41,13 +47,13 @@
 #    v  v v v
 # ____ __ _ _
 # CCSD-SD-1-1: Two-boson amplitudes and one-boson-one-fermion coupling
-ccsd = EBCC(
+ccsd_sd_1_1 = EBCC(
     mf,
     ansatz="CCSD-SD-1-1",
     omega=omega,
     g=g,
 )
-ccsd.kernel()
+ccsd_sd_1_1.kernel()
 
 #    ,--------- Fermionic ansatz
 #    |  ,------ Bosonic excitation amplitudes
@@ -56,10 +62,10 @@
 #    v  v v v
 # ____ __ _ _
 # CCSD-SD-1-2: Two-boson amplitudes and two-boson-one-fermion coupling
-ccsd = EBCC(
+ccsd_sd_1_2 = EBCC(
     mf,
     ansatz="CCSD-SD-1-2",
     omega=omega,
     g=g,
 )
-ccsd.kernel()
+ccsd_sd_1_2.kernel()

examples/02-eom_uccsd.py

@@ -1,18 +1,27 @@
+"""
+Example of a CCSD calculation using a UHF reference and a subsequent
+EOM-CCSD calculation for the ionization potential.
+"""
+
 import numpy as np
 from pyscf import gto, scf
 
 from ebcc import UEBCC
 
+# Define the molecule using PySCF
 mol = gto.Mole()
 mol.atom = "H 0 0 0; F 0 0 1.1"
 mol.basis = "cc-pvdz"
 mol.build()
 
+# Run a UHF calculation using PySCF
 mf = scf.UHF(mol)
 mf.kernel()
 
-ccsd = UEBCC(mf)
+# Run a UCCSD calculation
+ccsd = UEBCC(mf, ansatz="CCSD")
 ccsd.kernel()
 
+# Run an EOM-CCSD calculation
 eom = ccsd.ip_eom()
 eom.kernel()

examples/03-cc2.py

@@ -1,16 +1,22 @@
+"""
+Example of a CC2 calculation.
+"""
+
 import numpy as np
 from pyscf import gto, scf
 
 from ebcc import EBCC
 
+# Define the molecule using PySCF
 mol = gto.Mole()
 mol.atom = "H 0 0 0; F 0 0 1.1"
 mol.basis = "cc-pvdz"
 mol.build()
 
+# Run a RHF calculation using PySCF
 mf = scf.RHF(mol)
 mf.kernel()
 
-ccsd = EBCC(mf, ansatz="CC2")
-ccsd.kernel()
-
+# Run a CC2 calculation
+cc2 = EBCC(mf, ansatz="CC2")
+cc2.kernel()

examples/04-ccsdt_active_space.py

@@ -1,25 +1,35 @@
+"""
+Example of a CCSDt' calculation with T3 amplitudes in an active
+space.
+"""
+
 import numpy as np
 from pyscf import gto, scf
 
 from ebcc import REBCC, Space
 
+# Define the molecule using PySCF
 mol = gto.Mole()
 mol.atom = "H 0 0 0; F 0 0 1.1"
 mol.basis = "cc-pvdz"
 mol.build()
 
+# Run a RHF calculation using PySCF
 mf = scf.RHF(mol)
 mf.kernel()
 
-frozen = np.zeros_like(mf.mo_occ, dtype=bool)
-active = np.zeros_like(mf.mo_occ, dtype=bool)
+# Define the occupied, frozen, and active spaces
+occupied = mf.mo_occ > 0
+frozen = np.zeros_like(occupied)
+active = np.zeros_like(occupied)
 active[mol.nelectron // 2 - 1] = True  # HOMO
 active[mol.nelectron // 2] = True      # LUMO
 space = Space(
-    mf.mo_occ > 0,
+    occupied,
     frozen,
     active,
 )
 
+# Run a CCSDt' calculation
 ccsdt = REBCC(mf, ansatz="CCSDt'", space=space)
 ccsdt.kernel()

examples/05-bccd.py

@@ -0,0 +1,25 @@
+"""
+Example of a Brueckner orbital calculation using a CCSD reference.
+"""
+
+import numpy as np
+from pyscf import gto, scf
+
+from ebcc import EBCC
+
+# Define the molecule using PySCF
+mol = gto.Mole()
+mol.atom = "H 0 0 0; F 0 0 1.1"
+mol.basis = "cc-pvdz"
+mol.build()
+
+# Run a RHF calculation using PySCF
+mf = scf.RHF(mol)
+mf.kernel()
+
+# Run a CCSD calculation
+ccsd = EBCC(mf, ansatz="CCSD")
+ccsd.kernel()
+
+# Run a Brueckner orbital calculation using the CCSD reference
+ccsd.brueckner(e_tol=1e-6, t_tol=1e-5)

examples/06-restart.py

@@ -0,0 +1,37 @@
+"""
+Example of saving and restarting an EBCC calculation.
+"""
+
+import os
+import numpy as np
+from pyscf import gto, scf
+
+from ebcc import REBCC
+
+# Define the molecule using PySCF
+mol = gto.Mole()
+mol.atom = "H 0 0 0; F 0 0 1.1"
+mol.basis = "cc-pvdz"
+mol.build()
+
+# Run a RHF calculation using PySCF
+mf = scf.RHF(mol)
+mf.kernel()
+
+# Run a CC3 calculation that does not converge
+cc3 = REBCC(mf, ansatz="CC3")
+cc3.options.max_iter = 5
+cc3.kernel()
+
+# Save the calculation to a file
+cc3.write("restart.h5")
+
+# Load the calculation from the file
+cc3 = REBCC.read("restart.h5")
+
+# Run the calculation again, but this time with a higher max_iter
+cc3.options.max_iter = 20
+cc3.kernel()
+
+# Delete the file
+os.remove("restart.h5")

examples/07-rdms.py

@@ -0,0 +1,42 @@
+"""
+Example obtaining RDMs from EBCC.
+"""
+
+import numpy as np
+from pyscf import gto, scf, lib
+
+from ebcc import GEBCC
+
+# Define the molecule using PySCF
+mol = gto.Mole()
+mol.atom = "H 0 0 0; F 0 0 1.1"
+mol.basis = "cc-pvdz"
+mol.build()
+
+# Run a RHF calculation using PySCF
+mf = scf.RHF(mol)
+mf.kernel()
+
+# Convert the RHF object to a GHF object (not necessary for density
+# matrices, just simplifies the Hamiltonian)
+mf = mf.to_ghf()
+
+# Run a CCSD calculation
+ccsd = GEBCC(mf, ansatz="CCSD")
+ccsd.kernel()
+
+# If the Λ amplitudes are not solved, EBCC will use the approximation
+# Λ = T* and warn the user.
+ccsd.solve_lambda()
+
+# Fermionic RDMs
+dm1 = ccsd.make_rdm1_f()
+dm2 = ccsd.make_rdm2_f()
+
+# Compare the energies
+h1 = np.linalg.multi_dot((mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff))
+h2 = ccsd.get_eris().array
+e_rdm = lib.einsum("pq,qp->", h1, dm1)
+e_rdm += lib.einsum("pqrs,pqrs->", h2, dm2) * 0.5
+e_rdm += mol.energy_nuc()
+assert np.allclose(e_rdm, ccsd.e_tot)

examples/08-mp2.py

@@ -0,0 +1,27 @@
+"""
+Example of a simple MP2 calculation.
+"""
+
+import numpy as np
+from pyscf import gto, scf
+
+from ebcc import REBCC
+
+# Define the molecule using PySCF
+mol = gto.Mole()
+mol.atom = "H 0 0 0; F 0 0 1.1"
+mol.basis = "cc-pvdz"
+mol.build()
+
+# Run a Hartree-Fock calculation using PySCF
+mf = scf.RHF(mol)
+mf.kernel()
+
+# The CC solver can be used for Moller-Plesset perturbation theory,
+# and EBCC will detect that convergence isn't necessary via the
+# `ansatz.is_one_shot` attribute.
+mp2 = REBCC(mf, ansatz="MP2")
+mp2.kernel()
+
+# Note that this is not the most efficient way to compute MP energies,
+# as EBCC will still generate bare amplitudes stored in memory.

examples/09-spin_conversion.py

@@ -0,0 +1,36 @@
+"""
+Example of converting restricted and unrestricted EBCC calculations
+to spin-orbital (GHF) EBCC calculations.
+"""
+
+import numpy as np
+from pyscf import gto, scf
+
+from ebcc import REBCC, UEBCC, GEBCC
+
+# Define the molecule using PySCF
+mol = gto.Mole()
+mol.atom = "H 0 0 0; F 0 0 1.1"
+mol.basis = "cc-pvdz"
+mol.build()
+
+# Run a RHF calculation using PySCF
+rhf = scf.RHF(mol)
+rhf.kernel()
+
+# Run a REBCC calculation
+rcc = REBCC(rhf, ansatz="QCISD")
+rcc.kernel()
+
+# Convert to unrestricted and run kernel - unless the UEBCC solution
+# breaks some symmetry this should converge immediately to the same
+# solution as the REBCC calculation.
+uebcc_from_rebcc = UEBCC.from_rebcc(rcc)
+uebcc_from_rebcc.kernel()
+
+# Conversion of REBCC to GEBCC goes via a UEBCC intermediate, here
+# we just convert the UEBCC object we just created. Once again, in
+# the absence of symmetry breaking this should converge immediately
+# to the same solution as the REBCC calculation.
+gebcc_from_uebcc = GEBCC.from_uebcc(uebcc_from_rebcc)
+gebcc_from_uebcc.kernel()