1import sys 2sys.path.append('/usr/lib/python2.7/site-packages') 3import PyCheMPS2 4from pyscf import gto, scf, ao2mo, symm, fci 5import numpy as np 6import ctypes 7import time 8 9################## 10# Molecule 1 # 11################## 12mol1 = gto.Mole() 13mol1.atom = ''' 14 N 0.0000 0.0000 0.000000 15 N 0.0000 0.0000 1.090105 16 ''' 17mol1.basis = 'sto-3g' 18mol1.symmetry = 1 19mol1.charge = 0 20mol1.spin = 0 #2*S; multiplicity-1 21mol1.build() 22 23################## 24# Molecule 2 # 25################## 26mol2 = gto.Mole() 27mol2.atom = ''' 28 O 0.000000000 0.00000000 0.000000000 29 H 0.790689766 0.00000000 0.612217330 30 H -0.790689766 0.00000000 0.612217330 31 ''' 32mol2.basis = '6-31g' 33mol2.symmetry = 1 34mol2.charge = 0 35mol2.spin = 0 #2*S; multiplicity-1 36mol2.build() 37 38for geval in range(1): 39 40 if ( geval == 0 ): 41 mol = mol1 42 psi4group = 7 #d2h 43 print "####################" 44 print "# N2 in sto-3g #" 45 print "####################" 46 if ( geval == 1 ): 47 mol = mol2 48 psi4group = 5 #c2v 49 print "####################" 50 print "# H2O in 6-31g #" 51 print "####################" 52 mf = scf.RHF(mol) 53 mf.verbose = 3 54 mf.scf() 55 56 ############################# 57 # Build the Hamiltonian # 58 ############################# 59 L = mf.mol.nao_nr() 60 N = mf.mol.nelectron 61 CONST = mf.mol.energy_nuc() 62 OEImo = np.dot( np.dot( mf.mo_coeff.T, mf.mol.intor('cint1e_kin_sph') + mf.mol.intor('cint1e_nuc_sph') ), mf.mo_coeff ) 63 TEImo = ao2mo.outcore.full_iofree(mf.mol, mf.mo_coeff, compact=False).reshape(L,L,L,L) 64 orb2irrep = np.array(symm.label_orb_symm(mf.mol, mf.mol.irrep_id, mf.mol.symm_orb, mf.mo_coeff), dtype=ctypes.c_int) 65 for cnt in range( len( orb2irrep ) ): 66 orb2irrep[ cnt ] = orb2irrep[ cnt ] % 10 67 68 ####################################### 69 # PySCF FCI calculation and 4-RDM # 70 ####################################### 71 fci_solver = fci.FCI(mol, mf.mo_coeff) 72 fci_energy, fci_vector = fci_solver.kernel() 73 print "PySCF FCI energy =", fci_energy + CONST 74 pyscf_start = time.time() 75 dm1, dm2, dm3, dm4 = fci.rdm.make_dm1234('FCI4pdm_kern_spin0', fci_vector, fci_vector, L, N) 76 dm1, dm2, dm3, dm4 = fci.rdm.reorder_dm1234(dm1, dm2, dm3, dm4, inplace=True) 77 pyscf_trace = np.einsum( 'ii->', np.einsum( 'jjkl->kl', np.einsum( 'iijklm->jklm', np.einsum( 'iijklmno->jklmno', dm4 ) ) ) ) 78 pyscf_dm4 = np.array( dm4, copy=True ) 79 pyscf_dm4 = pyscf_dm4.swapaxes( 1, 2 ) # dm4[ 1, 5, 2, 6, 3, 7, 4, 8 ] -> dm4[ 1, 2, 5, 6, 3, 7, 4, 8 ] 80 pyscf_dm4 = pyscf_dm4.swapaxes( 2, 4 ) # dm4[ 1, 2, 5, 6, 3, 7, 4, 8 ] -> dm4[ 1, 2, 3, 6, 5, 7, 4, 8 ] 81 pyscf_dm4 = pyscf_dm4.swapaxes( 3, 6 ) # dm4[ 1, 2, 3, 6, 5, 7, 4, 8 ] -> dm4[ 1, 2, 3, 4, 5, 7, 6, 8 ] 82 pyscf_dm4 = pyscf_dm4.swapaxes( 5, 6 ) # dm4[ 1, 2, 3, 4, 5, 7, 6, 8 ] -> dm4[ 1, 2, 3, 4, 5, 6, 7, 8 ] 83 pyscf_end = time.time() 84 85 ############################# 86 # PyCheMPS2 Hamiltonian # 87 ############################# 88 Initializer = PyCheMPS2.PyInitialize() 89 Initializer.Init() 90 Ham = PyCheMPS2.PyHamiltonian( L, psi4group, orb2irrep, 'none' ) 91 Ham.setEconst( CONST ) 92 for cnt1 in range( L ): 93 for cnt2 in range( L ): 94 irrep12 = Ham.getOrbitalIrrep( cnt1 ) ^ Ham.getOrbitalIrrep( cnt2 ) 95 if ( irrep12 == 0 ): 96 Ham.setTmat(cnt1, cnt2, OEImo[cnt1, cnt2]) 97 for cnt3 in range( L ): 98 for cnt4 in range( L ): 99 irrep34 = Ham.getOrbitalIrrep( cnt3 ) ^ Ham.getOrbitalIrrep( cnt4 ) 100 if ( irrep12 == irrep34 ): 101 Ham.setVmat(cnt1, cnt2, cnt3, cnt4, TEImo[cnt1, cnt3, cnt2, cnt4]) #From chemist to physics notation 102 103 ########################################### 104 # PyCheMPS2 FCI calculation and 3-RDM # 105 ########################################### 106 Irrep = 0 107 Nel_up = ( N + mol.spin ) / 2 108 Nel_down = ( N - mol.spin ) / 2 109 work_mem = 25.0 110 verbosity = 0 111 theFCI = PyCheMPS2.PyFCI(Ham, Nel_up, Nel_down, Irrep, work_mem, verbosity) 112 GSvector = np.zeros([ theFCI.getVecLength() ], dtype=ctypes.c_double) 113 GSvector[ theFCI.LowestEnergyDeterminant() ] = 1.0 114 EnergyFCI = theFCI.GSDavidson(GSvector) 115 print "PyCheMPS2 FCI energy =", EnergyFCI 116 start2 = time.time() 117 FourRDM = np.zeros([ L**8 ], dtype=ctypes.c_double) 118 theFCI.Fill4RDM( GSvector, FourRDM ) 119 FourRDM = FourRDM.reshape(L,L,L,L,L,L,L,L) 120 print "PySCF Tr(\'dm4\') =", pyscf_trace 121 print "PyCheMPS2 Tr(4-RDM) =", np.einsum( 'ii->', np.einsum( 'jkjl->kl', np.einsum( 'ijkilm->jklm', np.einsum( 'ijkliqrs->jklqrs', FourRDM ) ) ) ) 122 end2 = time.time() 123 print "Time 4-RDM PySCF =", pyscf_end - pyscf_start, "seconds." 124 print "Time 4-RDM PyCheMPS2 =", end2 - start2, "seconds." 125 126 ############################## 127 # Compare the two 4-RDMs # 128 ############################## 129 for orb1 in range(L): 130 for orb2 in range(L): 131 for orb3 in range(L): 132 for orb4 in range(L): 133 for orb5 in range(L): 134 for orb6 in range(L): 135 for orb7 in range(L): 136 for orb8 in range(L): 137 temp = FourRDM[ orb1, orb2, orb3, orb4, orb5, orb6, orb7, orb8 ] - pyscf_dm4[ orb1, orb2, orb3, orb4, orb5, orb6, orb7, orb8 ] 138 if abs( temp ) > 1e-5: 139 print "4-RDM[",orb1,",",orb2,",",orb3,",",orb4,",",orb5,",",orb6,",",orb7,",",orb8,"] diff =", temp 140 print "RMS difference PySCF and PyCheMPS2 4-RDM =", np.linalg.norm( FourRDM - pyscf_dm4 ) 141 142 143 144