Wave Function Based Embedding (EWF)

This introduces the EWF class to perform a more generalized wave function based quantum embedding that improves on DMET for ab initio systems, as presented in Phys. Rev. X 12, 011046 [1].

Water Molecule

An embedded wave function calculation of a simple water molecule can be performed with the following code:

import pyscf
import pyscf.gto
import pyscf.scf
import pyscf.cc
import vayesta
import vayesta.ewf


mol = pyscf.gto.Mole()
mol.atom = """
O  0.0000   0.0000   0.1173
H  0.0000   0.7572  -0.4692
H  0.0000  -0.7572  -0.4692
"""
mol.basis = "cc-pVDZ"
mol.output = "pyscf.txt"
mol.build()

# Hartree-Fock
mf = pyscf.scf.RHF(mol)
mf.kernel()

# Embedded CCSD
emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6))
emb.kernel()

# Reference full system CCSD:
cc = pyscf.cc.CCSD(mf)
cc.kernel()

print("E(HF)=            %+16.8f Ha" % mf.e_tot)
print("E(CCSD)=          %+16.8f Ha" % cc.e_tot)
print("E(Emb. CCSD)[WF]= %+16.8f Ha" % emb.e_tot)
print("E(Emb. CCSD)[DM]= %+16.8f Ha" % emb.get_dm_energy())

There are two main differences in setup compared to the DMET embedding class in the definition of the EWF class (lines 24–25).

• The keyword argument bath_options=dict(threshold=1e-6) defines the threshold for the MP2 bath natural orbitals (the variable \(\eta\) in Ref. 1) [2].

• No fragmentation scheme is explicitly provided, which by default results in the system being fully fragmented into simple atomic fragments as defined by intrinsic atomic orbitals (IAOs). This is equivalent to adding the following lines before calling the embedding kernel method:

with emb.iao_fragmentation() as f:
    f.add_all_atomic_fragments()

which, in turn, is equivalent to

with emb.iao_fragmentation() as f:
    for atom in range(mol.natm):
        f.add_atomic_fragment(atom)

Other fragmentations schemes can still be used, but have to be specified manually.

Cubic Boron Nitride (cBN)

In this example we calculate cubic boron nitride (Zinc Blende structure):

Note

The basis set, auxiliary (density-fitting) basis set, and k-point sampling in this example are much too small for accurate results and only chosen for demonstration.

import numpy as np
import pyscf
import pyscf.pbc
import pyscf.pbc.scf
import pyscf.pbc.cc
import vayesta
import vayesta.ewf


cell = pyscf.pbc.gto.Cell()
a = 3.615
cell.atom = "B 0 0 0 ; N %f %f %f" % (a / 4, a / 4, a / 4)
cell.a = np.asarray([[a / 2, a / 2, 0], [0, a / 2, a / 2], [a / 2, 0, a / 2]])
cell.basis = "sto-6g"
cell.output = "pyscf.out"
cell.build()

# Hartree-Fock with k-points
kmesh = [2, 2, 2]
kpts = cell.make_kpts(kmesh)
mf = pyscf.pbc.scf.KRHF(cell, kpts)
mf = mf.density_fit(auxbasis="sto-6g")
mf.kernel()

# Full system CCSD
cc = pyscf.pbc.cc.KCCSD(mf)
cc.kernel()

# Embedded calculation will automatically fold the k-point sampled mean-field to the supercell
# solve_lambda=True is required, if density-matrix is needed!
emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6), solver_options=dict(solve_lambda=True))
emb.kernel()

print("E(HF)=             %+16.8f Ha" % mf.e_tot)
print("E(Emb. CCSD)=      %+16.8f Ha" % emb.e_tot)
print("E(CCSD)=           %+16.8f Ha" % cc.e_tot)

# One-body density matrix in the supercell AO-basis
dm1 = emb.make_rdm1(ao_basis=True)
# Population analysis (q: charge, s: spin)
# Possible options for local_orbitals: 'mulliken', 'lowdin', 'iao+pao', or custom N(AO) x N(AO) coefficient matrix
# orbital_resolved=True is used to print orbital resolved (rather than only atom resolved) analysis
emb.pop_analysis(dm1, local_orbitals="iao+pao", orbital_resolved=True)

In line 34 the setup of the embedding class is performed in the same way as for molecules. Vayesta will detect if the mean field object mf has k-points defined. If these are found, then the k-point sampled mean-field will automatically be folded to the \(\Gamma\)-point of the equivalent (in this case \(2\times2\times2\)) Born–von Karman supercell.

Note

Only Monkhorst-pack k-point meshes which include the \(\Gamma\)-point are currently supported.

Note that instantiating the embedding class with a k-point sampled mean-field object will automatically add the translational symmetry to the symmetry group stored in emb.symmetry. This assumes that the user will only define fragments within the original primitive unit cell, which are then copied throughout the supercell using the translational symmetry (and this symmetry will be exploited to avoid the cost of solving embedded fragments which are symmetrically equivalent to others). For calculations of fragments across multiple primitive cells or where the primitive cell has not been explicitly enlarged to encompass the full fragment space, the translational symmetry should be removed by calling emb.symmetry.clear_translations() or overwritten via emb.symmetry.set_translations(nimages), as demonstrated in for the 1D Hubbard model here.

Performing the embedding in the supercell allows for optimal utilization of the locality of electron correlation, as the embedding problems are only restricted to have the periodicity of the supercell, rather than the k-point sampled primitive cell. Properties, such as density-matrix calculated in line 42, will recover the full, primitive cell symmetry, since they are obtained from a summation over all symmetry equivalent fragments in the supercell. This is confirmed by the population analysis, which shows that the boron atom 2 has the same population than boron atom 0, despite being part of a different primitive cell within the supercell:

Population analysis
-------------------
    0 B:       q=  0.17325874  s=  0.00000000
         0 0 B 1s          =  1.98971008
         1 0 B 2s          =  0.76417671
         2 0 B 2px         =  0.69095149
         3 0 B 2py         =  0.69095149
         4 0 B 2pz         =  0.69095149
    1 N:       q= -0.17325874  s=  0.00000000
         5 1 N 1s          =  1.99053993
         6 1 N 2s          =  1.17403392
         7 1 N 2px         =  1.33622830
         8 1 N 2py         =  1.33622830
         9 1 N 2pz         =  1.33622830
    2 B:       q=  0.17325874  s=  0.00000000
        10 2 B 1s          =  1.98971008
        11 2 B 2s          =  0.76417671
        12 2 B 2px         =  0.69095149
        13 2 B 2py         =  0.69095149
        14 2 B 2pz         =  0.69095149
    3 N:       q= -0.17325874  s=  0.00000000
        15 3 N 1s          =  1.99053993
        16 3 N 2s          =  1.17403392
        17 3 N 2px         =  1.33622830
        18 3 N 2py         =  1.33622830
        19 3 N 2pz         =  1.33622830
    ...

[1]

Max Nusspickel and George H. Booth, Phys. Rev. X 12, 011046 (2022).

[2]

The definition of \(\eta\) in Vayesta differs from Ref. 1 by a factor of \(1/2\). In Ref. 1 the occupation numbers of the spin traced density matrix (between 0 and 2) are compared against \(\eta\), whereas in Vayesta the eigenvalues of the spinned density-matrix (between 0 and 1) are used. The latter definition allows for a more natural comparison between spin-restricted and unrestricted calculations with the same value of \(\eta\).