Wave Function Based Embedding (EWF)
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 7 8 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 14""" 15mol.basis = "cc-pVDZ" 16mol.output = "pyscf.txt" 17mol.build() 18 19# Hartree-Fock 20mf = pyscf.scf.RHF(mol) 21mf.kernel() 22 23# Embedded CCSD 24emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6)) 25emb.kernel() 26 27# Reference full system CCSD: 28cc = pyscf.cc.CCSD(mf) 29cc.kernel() 30 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).
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):
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 8 9 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([[a / 2, a / 2, 0], [0, a / 2, a / 2], [a / 2, 0, a / 2]]) 14cell.basis = "sto-6g" 15cell.output = "pyscf.out" 16cell.build() 17 18# Hartree-Fock with k-points 19kmesh = [2, 2, 2] 20kpts = cell.make_kpts(kmesh) 21mf = pyscf.pbc.scf.KRHF(cell, kpts) 22mf = mf.density_fit(auxbasis="sto-6g") 23mf.kernel() 24 25# Full system CCSD 26cc = pyscf.pbc.cc.KCCSD(mf) 27cc.kernel() 28 29# Embedded calculation will automatically fold the k-point sampled mean-field to the supercell 30# solve_lambda=True is required, if density-matrix is needed! 31emb = vayesta.ewf.EWF(mf, bath_options=dict(threshold=1e-6), solver_options=dict(solve_lambda=True)) 32emb.kernel() 33 34print("E(HF)= %+16.8f Ha" % mf.e_tot) 35print("E(Emb. CCSD)= %+16.8f Ha" % emb.e_tot) 36print("E(CCSD)= %+16.8f Ha" % cc.e_tot) 37 38# One-body density matrix in the supercell AO-basis 39dm1 = emb.make_rdm1(ao_basis=True) 40# Population analysis (q: charge, s: spin) 41# Possible options for local_orbitals: 'mulliken', 'lowdin', 'iao+pao', or custom N(AO) x N(AO) coefficient matrix 42# orbital_resolved=True is used to print orbital resolved (rather than only atom resolved) analysis 43emb.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
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
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 ...