1#!/usr/bin/env python 2# 3# Contributors: 4# Qiming Sun <osirpt.sun@gmail.com> 5# 6 7from functools import reduce 8import numpy 9import scipy.linalg 10from pyscf import scf 11from pyscf import gto 12from pyscf import mcscf 13from pyscf import dmrgscf 14from pyscf import mrpt 15# 16# Adjust mpi runtime schedular to execute the calculation with multi-processor 17# 18# NOTE DMRG-NEVPT2 requires about 10 GB memory per processor in this example 19# 20dmrgscf.settings.MPIPREFIX = 'mpirun -np 8' 21 22 23''' 24Triplet and quintet energy gap of Iron-Porphyrin molecule using DMRG-CASSCF 25and DMRG-NEVPT2 methods. DMRG is an approximate FCI solver. It can be used 26to handle large active space. This example is the next step to example 27018-dmet_cas_for_feporph.py 28''' 29 30# 31# Following 018-dmet_cas_for_feporph.py, we still use density matrix embedding 32# theory (DMET) to generate CASSCF initial guess. The active space includes 33# the Fe double d-shell, 4s shell, and the ligand N 2pz orbitals to describe 34# metal-ligand pi bond and pi backbond. 35# 36 37 38################################################## 39# 40# Define DMET active space 41# 42################################################## 43def dmet_cas(mc, dm, implst): 44 from pyscf import lo 45 nao = mc.mol.nao_nr() 46 ncore = mc.ncore 47 ncas = mc.ncas 48 nocc = ncore + ncas 49 nimp = len(implst) 50 nbath = ncas - nimp 51 corth = lo.orth.orth_ao(mol, method='meta_lowdin') 52 s = mol.intor_symmetric('cint1e_ovlp_sph') 53 cinv = numpy.dot(corth.T, s) 54 # 55 # Sum over spin-orbital DMs, then transform spin-free DM to orthogonal basis 56 # 57 dm = reduce(numpy.dot, (cinv, dm[0]+dm[1], cinv.T)) 58 59 # 60 # Decomposing DM to get impurity orbitals, doubly occupied core orbitals 61 # and entangled bath orbitals. Active space is consist of impurity plus 62 # truncated bath. 63 # 64 implst = numpy.asarray(implst) 65 notimp = numpy.asarray([i for i in range(nao) if i not in implst]) 66 occi, ui = scipy.linalg.eigh(-dm[implst][:,implst]) 67 occb, ub = scipy.linalg.eigh(-dm[notimp][:,notimp]) 68 bathorb = numpy.dot(corth[:,notimp], ub) 69 imporb = numpy.dot(corth[:,implst], ui) 70 mocore = bathorb[:,:ncore] 71 mocas = numpy.hstack((imporb, bathorb[:,ncore:ncore+nbath])) 72 moext = bathorb[:,ncore+nbath:] 73 74 # 75 # Restore core, active and external space to "canonical" form. Spatial 76 # symmetry is reserved in this canonicalization. 77 # 78 hf_orb = mc._scf.mo_coeff 79 fock = reduce(numpy.dot, (s, hf_orb*mc._scf.mo_energy, hf_orb.T, s)) 80 81 fockc = reduce(numpy.dot, (mocore.T, fock, mocore)) 82 e, u = scipy.linalg.eigh(fockc) 83 mocore = numpy.dot(mocore, u) 84 focka = reduce(numpy.dot, (mocas.T, fock, mocas)) 85 e, u = scipy.linalg.eigh(focka) 86 mocas = numpy.dot(mocas, u) 87 focke = reduce(numpy.dot, (moext.T, fock, moext)) 88 e, u = scipy.linalg.eigh(focke) 89 moext = numpy.dot(moext, u) 90 91 # 92 # Initial guess 93 # 94 mo_init = numpy.hstack((mocore, mocas, moext)) 95 return mo_init 96 97 98 99 100################################################## 101# 102# Quintet 103# 104################################################## 105 106mol = gto.Mole() 107mol.atom = [ 108 ['Fe', (0. , 0.0000 , 0.0000)], 109 ['N' , (1.9764 , 0.0000 , 0.0000)], 110 ['N' , (0.0000 , 1.9884 , 0.0000)], 111 ['N' , (-1.9764 , 0.0000 , 0.0000)], 112 ['N' , (0.0000 , -1.9884 , 0.0000)], 113 ['C' , (2.8182 , -1.0903 , 0.0000)], 114 ['C' , (2.8182 , 1.0903 , 0.0000)], 115 ['C' , (1.0918 , 2.8249 , 0.0000)], 116 ['C' , (-1.0918 , 2.8249 , 0.0000)], 117 ['C' , (-2.8182 , 1.0903 , 0.0000)], 118 ['C' , (-2.8182 , -1.0903 , 0.0000)], 119 ['C' , (-1.0918 , -2.8249 , 0.0000)], 120 ['C' , (1.0918 , -2.8249 , 0.0000)], 121 ['C' , (4.1961 , -0.6773 , 0.0000)], 122 ['C' , (4.1961 , 0.6773 , 0.0000)], 123 ['C' , (0.6825 , 4.1912 , 0.0000)], 124 ['C' , (-0.6825 , 4.1912 , 0.0000)], 125 ['C' , (-4.1961 , 0.6773 , 0.0000)], 126 ['C' , (-4.1961 , -0.6773 , 0.0000)], 127 ['C' , (-0.6825 , -4.1912 , 0.0000)], 128 ['C' , (0.6825 , -4.1912 , 0.0000)], 129 ['H' , (5.0441 , -1.3538 , 0.0000)], 130 ['H' , (5.0441 , 1.3538 , 0.0000)], 131 ['H' , (1.3558 , 5.0416 , 0.0000)], 132 ['H' , (-1.3558 , 5.0416 , 0.0000)], 133 ['H' , (-5.0441 , 1.3538 , 0.0000)], 134 ['H' , (-5.0441 , -1.3538 , 0.0000)], 135 ['H' , (-1.3558 , -5.0416 , 0.0000)], 136 ['H' , (1.3558 , -5.0416 , 0.0000)], 137 ['C' , (2.4150 , 2.4083 , 0.0000)], 138 ['C' , (-2.4150 , 2.4083 , 0.0000)], 139 ['C' , (-2.4150 , -2.4083 , 0.0000)], 140 ['C' , (2.4150 , -2.4083 , 0.0000)], 141 ['H' , (3.1855 , 3.1752 , 0.0000)], 142 ['H' , (-3.1855 , 3.1752 , 0.0000)], 143 ['H' , (-3.1855 , -3.1752 , 0.0000)], 144 ['H' , (3.1855 , -3.1752 , 0.0000)], 145] 146mol.basis = 'ccpvdz' 147mol.verbose = 4 148mol.output = 'fepor-dmrgscf.out' 149mol.spin = 4 150mol.symmetry = True 151mol.build() 152 153mf = scf.ROHF(mol) 154mf = scf.fast_newton(mf) 155 156# 157# CAS(16e, 20o) 158# 159# mcscf.approx_hessian approximates the orbital hessian. It does not affect 160# results. The N-2pz orbitals introduces more entanglement to environment. 161# 5 bath orbitals which have the strongest entanglement to impurity are 162# considered in active space. 163# 164mc = mcscf.approx_hessian(dmrgscf.dmrgci.DMRGSCF(mf, 20, 16)) 165# Function mol.search_ao_label returns the indices of the required AOs 166# It is equivalent to the following expression 167#idx = [i for i,s in enumerate(mol.ao_labels()) 168# if 'Fe 3d' in s or 'Fe 4d' in s or 'Fe 4s' in s or 'N 2pz' in s] 169idx = mol.search_ao_label(['Fe 3d', 'Fe 4d', 'Fe 4s', 'N 2pz']) 170mo = dmet_cas(mc, mf.make_rdm1(), idx) 171 172mc.fcisolver.wfnsym = 'Ag' 173mc.kernel(mo) 174#mc.analyze() 175e_q = mc.e_tot # -2244.90267106288 176cas_q = mc.mo_coeff[:,mc.ncore:mc.ncore+mc.ncas] 177 178# 179# call DMRG-NEVPT2 (about 2 days, 100 GB memory) 180# 181ept2_q = mrpt.NEVPT(mc).kernel() 182 183 184 185 186 187################################################## 188# 189# Triplet 190# 191################################################## 192 193mol.spin = 2 194mol.build(0, 0) 195 196mf = scf.ROHF(mol) 197mf = scf.fast_newton(mf) 198 199# 200# CAS(16e, 20o) 201# 202# Unlike CAS(8e, 11o) which is easily to draw 4s-character orbitals into the 203# active space, the larger active space, which includes 4s orbitals, does not 204# have such issue on MCSCF wfn. 205# 206mc = mcscf.approx_hessian(dmrgscf.dmrgci.DMRGSCF(mf, 20, 16)) 207idx = mol.search_ao_label(['Fe 3d', 'Fe 4d', 'Fe 4s', 'N 2pz']) 208mo = dmet_cas(mc, mf.make_rdm1(), idx3d) 209mc.fcisolver.wfnsym = 'B1g' 210mc.kernel(mo) 211mo = mc.mo_coeff 212#mc.analzye() 213e_t = mc.e_tot # -2244.88920313881 214cas_t = mc.mo_coeff[:,mc.ncore:mc.ncore+mc.ncas] 215 216# 217# call DMRG-NEVPT2 (about 2 days, 100 GB memory) 218# 219ept2_t = mrpt.NEVPT(mc).kernel() 220 221print('E(T) = %.15g E(Q) = %.15g gap = %.15g' % (e_t, e_q, e_t-e_q)) 222# E(T) = -2244.88920313881 E(Q) = -2244.90267106288 gap = 0.0134679240700279 223 224# The triplet and quintet active space are not perfectly overlaped 225s = reduce(numpy.dot, (cas_t.T, mol.intor('cint1e_ovlp_sph'), cas_q)) 226print('Active space overlpa <T|Q> ~ %f' % numpy.linalg.det(s)) 227 228print('NEVPT2: E(T) = %.15g E(Q) = %.15g' % (ept2_t, ept2_q)) 229# E(T) = -3.52155285166390 E(Q) = -3.46277436661638 230 231 232 233 234 235 236################################################## 237# 238# Output the active space orbitals to molden format 239# 240################################################## 241from pyscf import tools 242tools.molden.from_mo(mol, 'triplet-cas.molden', cas_t) 243tools.molden.from_mo(mol, 'quintet-cas.molden', cas_q) 244