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:

 1import pyscf
 2import pyscf.gto
 3import pyscf.scf
 4import pyscf.cc
 5import vayesta
 6import vayesta.ewf
 9mol = pyscf.gto.Mole()
10mol.atom = """
11O  0.0000   0.0000   0.1173
12H  0.0000   0.7572  -0.4692
13H  0.0000  -0.7572  -0.4692
15mol.basis = 'cc-pVDZ'
16mol.output = 'pyscf.txt'
19# Hartree-Fock
20mf = pyscf.scf.RHF(mol)
23# Embedded CCSD
24emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6))
27# Reference full system CCSD:
28cc = pyscf.cc.CCSD(mf)
31print("E(HF)=            %+16.8f Ha" % mf.e_tot)
32print("E(CCSD)=          %+16.8f Ha" % cc.e_tot)
33print("E(Emb. CCSD)[WF]= %+16.8f Ha" % emb.e_tot)
34print("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:

which, in turn, is equivalent to

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

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):


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.

 1import numpy as np
 2import pyscf
 3import pyscf.pbc
 4import pyscf.pbc.scf
 5import pyscf.pbc.cc
 6import vayesta
 7import vayesta.ewf
10cell = pyscf.pbc.gto.Cell()
11a = 3.615
12cell.atom = 'B 0 0 0 ; N %f %f %f' % (a/4, a/4, a/4)
13cell.a = np.asarray([
14    [a/2, a/2, 0],
15    [0, a/2, a/2],
16    [a/2, 0, a/2]])
17cell.basis = 'sto-6g'
18cell.output = 'pyscf.out'
21# Hartree-Fock with k-points
22kmesh = [2,2,2]
23kpts = cell.make_kpts(kmesh)
24mf = pyscf.pbc.scf.KRHF(cell, kpts)
25mf = mf.density_fit(auxbasis='sto-6g')
28# Full system CCSD
29cc = pyscf.pbc.cc.KCCSD(mf)
32# Embedded calculation will automatically fold the k-point sampled mean-field to the supercell
33# solve_lambda=True is required, if density-matrix is needed!
34emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6), solver_options=dict(solve_lambda=True))
37print("E(HF)=             %+16.8f Ha" % mf.e_tot)
38print("E(Emb. CCSD)=      %+16.8f Ha" % emb.e_tot)
39print("E(CCSD)=           %+16.8f Ha" % cc.e_tot)
41# One-body density matrix in the supercell AO-basis
42dm1 = emb.make_rdm1(ao_basis=True)
43# Population analysis (q: charge, s: spin)
44# Possible options for local_orbitals: 'mulliken', 'lowdin', 'iao+pao', or custom N(AO) x N(AO) coefficient matrix
45# orbital_resolved=True is used to print orbital resolved (rather than only atom resolved) analysis
46emb.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.


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