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