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