diff --git a/src/utils_periodic/Gen_Ezfio_from_integral_complex_3idx.sh b/src/utils_periodic/Gen_Ezfio_from_integral_complex_3idx.sh
new file mode 100755
index 00000000..10d42223
--- /dev/null
+++ b/src/utils_periodic/Gen_Ezfio_from_integral_complex_3idx.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+
+ezfio=$1
+# Create the integral
+echo 'Create Integral'
+
+echo 'Create EZFIO'
+read nel nmo natom <<< $(cat param) 
+read e_nucl <<< $(cat e_nuc)
+read nao <<< $(cat num_ao)
+read nkpts <<< $(cat num_kpts)
+read ndf <<< $(cat num_df)
+#./create_ezfio_complex_4idx.py $ezfio $nel $natom $nmo $e_nucl $nao $nkpts
+./create_ezfio_complex_3idx.py $ezfio $nel $natom $nmo $e_nucl $nao $nkpts $ndf
+#Handle the orbital consitensy check
+qp_edit -c $ezfio &> /dev/null
+cp $ezfio/{ao,mo}_basis/ao_md5 
+
+#Read the integral
+echo 'Read Integral'
+
+
+################################################
+##  using AO mono, 4-idx from pyscf           ##
+################################################
+qp_run import_integrals_ao_periodic $ezfio 
+
+
+################################################
+##  using AO mono, 3-idx, mo coef from pyscf  ##
+################################################
+
+#qp_run read_ao_mono_complex $ezfio 
+#qp_run read_kconserv $ezfio
+#qp_run read_ao_df_complex $ezfio
+#qp_run read_mo_coef_complex $ezfio    #start from converged pyscf MOs
+#
+#qp_run save_mo_df_to_disk $ezfio
+#qp_run save_mo_bielec_to_disk $ezfio
+
+#qp_run mo_from_ao_orth $ezfio        #use canonical orthonormalized AOs as initial MO guess
+#qp_run print_H_matrix_restart $ezfio > hmat.out
+
+
+###############################################################
+##  using AO mono, full 4-idx AO bielec, mo coef from pyscf  ##
+###############################################################
+
+#qp_run read_ao_mono_complex $ezfio 
+#qp_run read_kconserv $ezfio
+#qp_run read_ao_eri_chunk_complex $ezfio 
+#qp_run read_mo_coef_complex $ezfio    #start from converged pyscf MOs
+##qp_run mo_from_ao_orth $ezfio        #use canonical orthonormalized AOs as initial MO guess
+
+
+######################################################
+##  using MO mono, full 4-idx MO bielec from pyscf  ##
+######################################################
+
+#qp_run read_mo_mono_complex $ezfio 
+#qp_run read_kconserv $ezfio
+#qp_run read_mo_eri_chunk_complex $ezfio 
+
diff --git a/src/utils_periodic/MolPyscfToQPkpts.py b/src/utils_periodic/MolPyscfToQPkpts.py
new file mode 100644
index 00000000..1f72788e
--- /dev/null
+++ b/src/utils_periodic/MolPyscfToQPkpts.py
@@ -0,0 +1,881 @@
+import numpy as np
+from functools import reduce
+
+
+def memoize(f):
+    memo = {}
+    def helper(x):
+        if x not in memo:
+            memo[x] = f(x)
+        return memo[x]
+    return helper
+
+@memoize
+def idx2_tri(iijj):
+    '''
+    iijj should be a 2-tuple
+    return triangular compound index for (0-indexed counting)
+    '''
+    ij1=min(iijj)
+    ij2=max(iijj)
+    return ij1+(ij2*(ij2+1))//2
+#    return ij1+(ij2*(ij2-1))//2
+
+def pad(arr_in,outshape):
+    arr_out = np.zeros(outshape,dtype=np.complex128)
+    dataslice = tuple(slice(0,arr_in.shape[dim]) for dim in range(len(outshape)))
+    arr_out[dataslice] = arr_in
+    return arr_out
+
+def makesq(vlist,n1,n2):
+    '''
+    make hermitian matrices of size (n2 x n2) from from lower triangles
+    vlist is n1 lower triangles in flattened form
+    given: ([a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t],2,4)
+           output a 2x4x4 array, where each 4x4 is the square constructed from the lower triangle
+    [
+     [
+      [a  b* d* g*]
+      [b  c  e* h*]
+      [d  e  f  i*]
+      [g  h  i  j ]
+     ],
+     [
+      [k  l* n* q*]
+      [l  m  o* r*]
+      [n  o  p  s*]
+      [q  r  s  t ]
+     ]
+    ]
+    '''
+    out=np.zeros([n1,n2,n2],dtype=np.complex128)
+    n0 = vlist.shape[0]
+    lmask=np.tri(n2,dtype=bool)
+    for i in range(n0):
+        out[i][lmask] = vlist[i].conj()
+    out2=out.transpose([0,2,1])
+    for i in range(n0):
+        out2[i][lmask] = vlist[i]
+    return out2
+
+
+def makesq3(vlist,n2):
+    '''
+    make hermitian matrices of size (n2 x n2) from from lower triangles
+    vlist is n1 lower triangles in flattened form
+    given: ([a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t],2,4)
+           output a 2x4x4 array, where each 4x4 is the square constructed from the lower triangle
+    [
+     [
+      [a  b* d* g*]
+      [b  c  e* h*]
+      [d  e  f  i*]
+      [g  h  i  j ]
+     ],
+     [
+      [k  l* n* q*]
+      [l  m  o* r*]
+      [n  o  p  s*]
+      [q  r  s  t ]
+     ]
+    ]
+    '''
+    n0 = vlist.shape[0]
+    out=np.zeros([n0,n2,n2],dtype=np.complex128)
+    lmask=np.tri(n2,dtype=bool)
+    for i in range(n0):
+        out[i][lmask] = vlist[i].conj()
+    out2=out.transpose([0,2,1])
+    for i in range(n0):
+        out2[i][lmask] = vlist[i]
+    return out2
+
+def makesq2(vlist,n1,n2):
+    out=np.zeros([n1,n2,n2],dtype=np.complex128)
+    lmask=np.tri(n2,dtype=bool)
+    tmp=np.zeros([n2,n2],dtype=np.complex128)
+    tmp2=np.zeros([n2,n2],dtype=np.complex128)
+    for i in range(n1):
+        tmp[lmask] = vlist[i].conj()
+        tmp2=tmp.T
+        tmp2[lmask] = vlist[i]
+        out[i] = tmp2.copy()
+    return out
+
+
+def get_phase(cell, kpts, kmesh=None):
+    '''
+    The unitary transformation that transforms the supercell basis k-mesh
+    adapted basis.
+    '''
+    from pyscf.pbc import tools
+    from pyscf import lib
+
+    latt_vec = cell.lattice_vectors()
+    if kmesh is None:
+        # Guess kmesh
+        scaled_k = cell.get_scaled_kpts(kpts).round(8)
+        kmesh = (len(np.unique(scaled_k[:,0])),
+                 len(np.unique(scaled_k[:,1])),
+                 len(np.unique(scaled_k[:,2])))
+
+    R_rel_a = np.arange(kmesh[0])
+    R_rel_b = np.arange(kmesh[1])
+    R_rel_c = np.arange(kmesh[2])
+    R_vec_rel = lib.cartesian_prod((R_rel_a, R_rel_b, R_rel_c))
+    R_vec_abs = np.einsum('nu, uv -> nv', R_vec_rel, latt_vec)
+
+    NR = len(R_vec_abs)
+    phase = np.exp(1j*np.einsum('Ru, ku -> Rk', R_vec_abs, kpts))
+    phase /= np.sqrt(NR)  # normalization in supercell
+
+    # R_rel_mesh has to be construct exactly same to the Ts in super_cell function
+    scell = tools.super_cell(cell, kmesh)
+    return scell, phase
+
+def mo_k2gamma(cell, mo_energy, mo_coeff, kpts, kmesh=None):
+    '''
+    Transform MOs in Kpoints to the equivalents supercell
+    '''
+    from pyscf import lib
+    import scipy.linalg as la
+    scell, phase = get_phase(cell, kpts, kmesh)
+
+    E_g = np.hstack(mo_energy)
+    C_k = np.asarray(mo_coeff)
+    Nk, Nao, Nmo = C_k.shape
+    NR = phase.shape[0]
+
+    # Transform AO indices
+    C_gamma = np.einsum('Rk, kum -> Rukm', phase, C_k)
+    C_gamma = C_gamma.reshape(Nao*NR, Nk*Nmo)
+
+    E_sort_idx = np.argsort(E_g)
+    E_g = E_g[E_sort_idx]
+    C_gamma = C_gamma[:,E_sort_idx]
+    s = scell.pbc_intor('int1e_ovlp')
+    assert(abs(reduce(np.dot, (C_gamma.conj().T, s, C_gamma))
+               - np.eye(Nmo*Nk)).max() < 1e-7)
+
+    # Transform MO indices
+    E_k_degen = abs(E_g[1:] - E_g[:-1]).max() < 1e-5
+    if np.any(E_k_degen):
+        degen_mask = np.append(False, E_k_degen) | np.append(E_k_degen, False)
+        shift = min(E_g[degen_mask]) - .1
+        f = np.dot(C_gamma[:,degen_mask] * (E_g[degen_mask] - shift),
+                   C_gamma[:,degen_mask].conj().T)
+        assert(abs(f.imag).max() < 1e-5)
+
+        e, na_orb = la.eigh(f.real, s, type=2)
+        C_gamma[:,degen_mask] = na_orb[:, e>0]
+
+    if abs(C_gamma.imag).max() < 1e-7:
+        print('!Warning  Some complexe pollutions in MOs are present')
+
+    C_gamma = C_gamma.real
+    if  abs(reduce(np.dot, (C_gamma.conj().T, s, C_gamma)) - np.eye(Nmo*Nk)).max() < 1e-7:
+        print('!Warning  Some complexe pollutions in MOs are present')
+
+    s_k = cell.pbc_intor('int1e_ovlp', kpts=kpts)
+    # overlap between k-point unitcell and gamma-point supercell
+    s_k_g = np.einsum('kuv,Rk->kuRv', s_k, phase.conj()).reshape(Nk,Nao,NR*Nao)
+    # The unitary transformation from k-adapted orbitals to gamma-point orbitals
+    mo_phase = lib.einsum('kum,kuv,vi->kmi', C_k.conj(), s_k_g, C_gamma)
+
+    return mo_phase
+
+def qp2rename():
+    import shutil
+    qp2names={}
+    qp2names['mo_coef_complex'] = 'C.qp'
+    qp2names['bielec_ao_complex'] = 'W.qp'
+
+    qp2names['kinetic_ao_complex'] = 'T.qp'
+    qp2names['ne_ao_complex'] = 'V.qp'
+    qp2names['overlap_ao_complex'] = 'S.qp'
+
+
+    for old,new in qp2names.items():
+        shutil.move(old,new)
+    shutil.copy('e_nuc','E.qp')
+
+def pyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, 
+        print_ao_ints_bi=False, 
+        print_mo_ints_bi=False, 
+        print_ao_ints_df=True, 
+        print_mo_ints_df=False, 
+        print_ao_ints_mono=True, 
+        print_mo_ints_mono=False):
+    '''
+    kpts = List of kpoints coordinates. Cannot be null, for gamma is other script
+    kmesh = Mesh of kpoints (optional)
+    cas_idx = List of active MOs. If not specified all MOs are actives
+    int_threshold = The integral will be not printed in they are bellow that
+    '''
+  
+    from pyscf.pbc import ao2mo
+    from pyscf.pbc import tools
+    from pyscf.pbc.gto import ecp
+    import h5py
+    
+    mo_coef_threshold = int_threshold
+    ovlp_threshold = int_threshold
+    kin_threshold = int_threshold
+    ne_threshold = int_threshold
+    bielec_int_threshold = int_threshold
+
+    natom = len(cell.atom_coords())
+    print('n_atom per kpt',   natom)
+    print('num_elec per kpt', cell.nelectron)
+  
+    mo_coeff = mf.mo_coeff
+    # Mo_coeff actif
+    mo_k = np.array([c[:,cas_idx] for c in mo_coeff] if cas_idx is not None else mo_coeff)
+    e_k =  np.array([e[cas_idx] for e in mf.mo_energy] if cas_idx is not None else mf.mo_energy)
+  
+    Nk, nao, nmo = mo_k.shape
+    print("n Kpts", Nk)
+    print("n active Mos per kpt", nmo)
+    print("n AOs per kpt", nao)
+
+    naux = mf.with_df.auxcell.nao
+    print("n df fitting functions", naux)
+    with open('num_df','w') as f:
+        f.write(str(naux))
+  
+    # Write all the parameter need to creat a dummy EZFIO folder who will containt the integral after.
+    # More an implentation detail than a real thing
+    with open('param','w') as f:
+    # Note the use of nmo_tot
+        f.write(' '.join(map(str,(cell.nelectron*Nk, Nk*nmo, natom*Nk))))
+  
+    with open('num_ao','w') as f:
+        f.write(str(nao*Nk))
+    with open('num_kpts','w') as f:
+        f.write(str(Nk))
+    #                             _                             
+    # |\ |      _ |  _   _. ._   |_)  _  ._      |  _ o  _  ._  
+    # | \| |_| (_ | (/_ (_| |    | \ (/_ |_) |_| | _> | (_) | | 
+    #                                    |                      
+    
+    #Total energy shift due to Ewald probe charge = -1/2 * Nelec*madelung/cell.vol =
+    shift = tools.pbc.madelung(cell, kpts)*cell.nelectron * -.5 
+    e_nuc = (cell.energy_nuc() + shift)*Nk
+  
+    print('nucl_repul', e_nuc)
+    with open('e_nuc','w') as f:
+        f.write(str(e_nuc))
+  
+  
+  
+    #       __    __          _                                 
+    # |\/| |  |  |    _   _  |_  _ 
+    # |  | |__|  |__ (_) (/_ |  _> 
+    #                                               
+    with open('mo_coef_complex','w') as outfile:
+        c_kpts = np.reshape(mo_k,(Nk,nao,nmo))
+
+        for ik in range(Nk):
+            shift1=ik*nao+1
+            shift2=ik*nmo+1
+            for i in range(nao):
+                for j in range(nmo):
+                    cij = c_kpts[ik,i,j]
+                    if abs(cij) > mo_coef_threshold:
+                        outfile.write('%s %s %s %s\n' % (i+shift1, j+shift2, cij.real, cij.imag))
+    
+    # ___                                              
+    #  |  ._ _|_  _   _  ._ _. |  _   |\/|  _  ._   _  
+    # _|_ | | |_ (/_ (_| | (_| | _>   |  | (_) | | (_) 
+    #                 _|                              
+   
+    if mf.cell.pseudo:
+        v_kpts_ao = np.reshape(mf.with_df.get_pp(kpts=kpts),(Nk,nao,nao))
+    else:
+        v_kpts_ao = np.reshape(mf.with_df.get_nuc(kpts=kpts),(Nk,nao,nao))
+    if len(cell._ecpbas) > 0:
+        v_kpts_ao += np.reshape(ecp.ecp_int(cell, kpts),(Nk,nao,nao))
+
+    ne_ao = ('ne',v_kpts_ao,ne_threshold)
+    ovlp_ao = ('overlap',np.reshape(mf.get_ovlp(cell=cell,kpts=kpts),(Nk,nao,nao)),ovlp_threshold)
+    kin_ao = ('kinetic',np.reshape(cell.pbc_intor('int1e_kin',1,1,kpts=kpts),(Nk,nao,nao)),kin_threshold)
+
+    for name, intval_kpts_ao, thresh in (ne_ao, ovlp_ao, kin_ao):
+        if print_ao_ints_mono:
+            with open('%s_ao_complex' % name,'w') as outfile:
+                for ik in range(Nk):
+                    shift=ik*nao+1
+                    for i in range(nao):
+                        for j in range(i,nao):
+                            int_ij = intval_kpts_ao[ik,i,j]
+                            if abs(int_ij) > thresh:
+                                outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
+        if print_mo_ints_mono:
+            intval_kpts_mo = np.einsum('kim,kij,kjn->kmn',mo_k.conj(),intval_kpts_ao,mo_k)
+            with open('%s_mo_complex' % name,'w') as outfile:
+                for ik in range(Nk):
+                    shift=ik*nmo+1
+                    for i in range(nmo):
+                        for j in range(i,nmo):
+                            int_ij = intval_kpts_mo[ik,i,j]
+                            if abs(int_ij) > thresh:
+                                outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
+
+  
+    # ___                              _    
+    #  |  ._ _|_  _   _  ._ _. |  _   |_) o 
+    # _|_ | | |_ (/_ (_| | (_| | _>   |_) | 
+    #                 _|                    
+    #
+    kconserv = tools.get_kconserv(cell, kpts)
+
+    with open('kconserv_complex','w') as outfile:
+        for a in range(Nk):
+            for b in range(Nk):
+                for c in range(Nk):
+                    d = kconserv[a,b,c]
+                    outfile.write('%s %s %s %s\n' % (a+1,c+1,b+1,d+1))
+    
+
+    intfile=h5py.File(mf.with_df._cderi,'r')
+
+    j3c = intfile.get('j3c')
+    naosq = nao*nao
+    naotri = (nao*(nao+1))//2
+    j3ckeys = list(j3c.keys())
+    j3ckeys.sort(key=lambda strkey:int(strkey))
+
+    # in new(?) version of PySCF, there is an extra layer of groups before the datasets
+    # datasets used to be [/j3c/0,   /j3c/1,   /j3c/2,   ...]
+    # datasets now are    [/j3c/0/0, /j3c/1/0, /j3c/2/0, ...]
+    j3clist = [j3c.get(i+'/0') for i in j3ckeys]
+    if j3clist==[None]*len(j3clist):
+    # if using older version, stop before last level
+        j3clist = [j3c.get(i) for i in j3ckeys]
+
+    nkinvsq = 1./np.sqrt(Nk)
+
+    # dimensions are (kikj,iaux,jao,kao), where kikj is compound index of kpts i and j
+    # output dimensions should be reversed (nao, nao, naux, nkptpairs)
+    j3arr=np.array([(i.value.reshape([-1,nao,nao]) if (i.shape[1] == naosq) else makesq3(i.value,nao)) * nkinvsq for i in j3clist])
+
+    nkpt_pairs = j3arr.shape[0]
+
+    if print_ao_ints_df:
+        with open('df_ao_integral_array','w') as outfile:
+            pass
+        with open('df_ao_integral_array','a') as outfile:
+            for k,kpt_pair in enumerate(j3arr):
+                for iaux,dfbasfunc in enumerate(kpt_pair):
+                    for i,i0 in enumerate(dfbasfunc):
+                        for j,v in enumerate(i0):
+                            if (abs(v) > bielec_int_threshold):
+                                outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
+
+    if print_mo_ints_df:
+        kpair_list=[]
+        for i in range(Nk):
+            for j in range(Nk):
+                if(i>=j):
+                    kpair_list.append((i,j,idx2_tri((i,j))))
+        j3mo = np.array([np.einsum('mij,ik,jl->mkl',j3arr[kij],mo_k[ki].conj(),mo_k[kj]) for ki,kj,kij in kpair_list])
+        with open('df_mo_integral_array','w') as outfile:
+            pass
+        with open('df_mo_integral_array','a') as outfile:
+            for k,kpt_pair in enumerate(j3mo):
+                for iaux,dfbasfunc in enumerate(kpt_pair):
+                    for i,i0 in enumerate(dfbasfunc):
+                        for j,v in enumerate(i0):
+                            if (abs(v) > bielec_int_threshold):
+                                outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
+
+
+
+
+#    eri_4d_ao = np.zeros((Nk,nao,Nk,nao,Nk,nao,Nk,nao), dtype=np.complex)
+#    for d, kd in enumerate(kpts):
+#        for c, kc in enumerate(kpts):
+#            if c > d: break
+#            idx2_cd = idx2_tri(c,d)
+#            for b, kb in enumerate(kpts):
+#                if b > d: break
+#                a = kconserv[b,c,d]
+#                if idx2_tri(a,b) > idx2_cd: continue
+#                if ((c==d) and (a>b)): continue
+#                ka = kpts[a]
+#                v = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
+#                v *= 1./Nk
+#                eri_4d_ao[a,:,b,:,c,:,d] = v
+#    
+#    eri_4d_ao = eri_4d_ao.reshape([Nk*nao]*4)
+
+
+    if (print_ao_ints_bi or print_mo_ints_bi):
+        if print_ao_ints_bi:
+            with open('bielec_ao_complex','w') as outfile: 
+                pass
+        if print_mo_ints_bi:
+            with open('bielec_mo_complex','w') as outfile: 
+                pass
+        for d, kd in enumerate(kpts):
+            for c, kc in enumerate(kpts):
+                if c > d: break
+                idx2_cd = idx2_tri((c,d))
+                for b, kb in enumerate(kpts):
+                    if b > d: break
+                    a = kconserv[b,c,d]
+                    if idx2_tri((a,b)) > idx2_cd: continue
+                    if ((c==d) and (a>b)): continue
+                    ka = kpts[a]
+
+                    if print_ao_ints_bi:
+                        with open('bielec_ao_complex','a') as outfile:
+                            eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
+                            eri_4d_ao_kpt *= 1./Nk
+                            for l in range(nao):
+                                ll=l+d*nao
+                                for j in range(nao):
+                                    jj=j+c*nao
+                                    if jj>ll: break
+                                    idx2_jjll = idx2_tri((jj,ll))
+                                    for k in range(nao):
+                                        kk=k+b*nao
+                                        if kk>ll: break
+                                        for i in range(nao):
+                                            ii=i+a*nao
+                                            if idx2_tri((ii,kk)) > idx2_jjll: break
+                                            if ((jj==ll) and (ii>kk)): break
+                                            v=eri_4d_ao_kpt[i,k,j,l]
+                                            if (abs(v) > bielec_int_threshold):
+                                                outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
+            
+                    if print_mo_ints_bi:
+                        with open('bielec_mo_complex','a') as outfile:
+                            eri_4d_mo_kpt = mf.with_df.ao2mo([mo_k[a], mo_k[b], mo_k[c], mo_k[d]],
+                                                              [ka,kb,kc,kd],compact=False).reshape((nmo,)*4)
+                            eri_4d_mo_kpt *= 1./Nk
+                            for l in range(nmo):
+                                ll=l+d*nmo
+                                for j in range(nmo):
+                                    jj=j+c*nmo
+                                    if jj>ll: break
+                                    idx2_jjll = idx2_tri((jj,ll))
+                                    for k in range(nmo):
+                                        kk=k+b*nmo
+                                        if kk>ll: break
+                                        for i in range(nmo):
+                                            ii=i+a*nmo
+                                            if idx2_tri((ii,kk)) > idx2_jjll: break
+                                            if ((jj==ll) and (ii>kk)): break
+                                            v=eri_4d_mo_kpt[i,k,j,l]
+                                            if (abs(v) > bielec_int_threshold):
+                                                outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
+
+
+def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8, 
+        print_ao_ints_bi=False, 
+        print_mo_ints_bi=False, 
+        print_ao_ints_df=True, 
+        print_mo_ints_df=False, 
+        print_ao_ints_mono=True, 
+        print_mo_ints_mono=False):
+    '''
+    kpts = List of kpoints coordinates. Cannot be null, for gamma is other script
+    kmesh = Mesh of kpoints (optional)
+    cas_idx = List of active MOs. If not specified all MOs are actives
+    int_threshold = The integral will be not printed in they are bellow that
+    '''
+  
+    from pyscf.pbc import ao2mo
+    from pyscf.pbc import tools
+    from pyscf.pbc.gto import ecp
+    import h5py
+    import scipy
+
+    qph5=h5py.File('qpdat.h5')
+    qph5.create_group('nuclei')
+    qph5.create_group('electrons')
+    qph5.create_group('ao_basis')
+    qph5.create_group('mo_basis')
+
+    
+    mo_coef_threshold = int_threshold
+    ovlp_threshold = int_threshold
+    kin_threshold = int_threshold
+    ne_threshold = int_threshold
+    bielec_int_threshold = int_threshold
+
+    natom = cell.natm
+    nelec = cell.nelectron
+    print('n_atom per kpt',   natom)
+    print('num_elec per kpt', nelec)
+  
+    mo_coeff = mf.mo_coeff
+    # Mo_coeff actif
+    mo_k = np.array([c[:,cas_idx] for c in mo_coeff] if cas_idx is not None else mo_coeff)
+    e_k =  np.array([e[cas_idx] for e in mf.mo_energy] if cas_idx is not None else mf.mo_energy)
+  
+    Nk, nao, nmo = mo_k.shape
+    print("n Kpts", Nk)
+    print("n active Mos per kpt", nmo)
+    print("n AOs per kpt", nao)
+
+    naux = mf.with_df.auxcell.nao
+    print("n df fitting functions", naux)
+  
+    #in old version: param < nelec*Nk, nmo*Nk, natom*Nk
+    qph5['electrons'].attrs['elec_alpha_num']=nelec*Nk 
+    qph5['electrons'].attrs['elec_beta_num']=nelec*Nk
+    qph5['mo_basis'].attrs['mo_num']=Nk*nmo
+    qph5['ao_basis'].attrs['ao_num']=Nk*nao
+    qph5['nuclei'].attrs['nucl_num']=Nk*natom
+    qph5['nuclei'].attrs['kpt_num']=Nk
+    qph5['ao_basis'].attrs['df_num']=naux
+
+    #                             _                             
+    # |\ |      _ |  _   _. ._   |_)  _  ._      |  _ o  _  ._  
+    # | \| |_| (_ | (/_ (_| |    | \ (/_ |_) |_| | _> | (_) | | 
+    #                                    |                      
+    
+    #Total energy shift due to Ewald probe charge = -1/2 * Nelec*madelung/cell.vol =
+    shift = tools.pbc.madelung(cell, kpts)*cell.nelectron * -.5 
+    e_nuc = (cell.energy_nuc() + shift)*Nk
+  
+    print('nucl_repul', e_nuc)
+    qph5['nuclei'].attrs['nuclear_repulsion']=e_nuc
+  
+    #       __    __          _                                 
+    # |\/| |  |  |    _   _  |_  _ 
+    # |  | |__|  |__ (_) (/_ |  _> 
+    #                                               
+    mo_coef_blocked=scipy.linalg.block_diag(*mo_k)
+    qph5.create_dataset('mo_basis/mo_coef_real',data=mo_coef_blocked.real)
+    qph5.create_dataset('mo_basis/mo_coef_imag',data=mo_coef_blocked.imag)
+    qph5.create_dataset('mo_basis/mo_coef_kpts_real',data=mo_k.real)
+    qph5.create_dataset('mo_basis/mo_coef_kpts_imag',data=mo_k.imag)
+
+    # ___                                              
+    #  |  ._ _|_  _   _  ._ _. |  _   |\/|  _  ._   _  
+    # _|_ | | |_ (/_ (_| | (_| | _>   |  | (_) | | (_) 
+    #                 _|                              
+   
+    if mf.cell.pseudo:
+        ao_n_e = np.reshape(mf.with_df.get_pp(kpts=kpts),(Nk,nao,nao))
+    else:
+        v_kpts_ao = np.reshape(mf.with_df.get_nuc(kpts=kpts),(Nk,nao,nao))
+    if len(cell._ecpbas) > 0:
+        v_kpts_ao += np.reshape(ecp.ecp_int(cell, kpts),(Nk,nao,nao))
+      
+    ne_ao = ('ne',v_kpts_ao,ne_threshold)
+    ovlp_ao = ('overlap',np.reshape(mf.get_ovlp(cell=cell,kpts=kpts),(Nk,nao,nao)),ovlp_threshold)
+    kin_ao = ('kinetic',np.reshape(cell.pbc_intor('int1e_kin',1,1,kpts=kpts),(Nk,nao,nao)),kin_threshold)
+
+    for name, intval_kpts_ao, thresh in (ne_ao, ovlp_ao, kin_ao):
+        if print_ao_ints_mono:
+            with open('%s_ao_complex' % name,'w') as outfile:
+                for ik in range(Nk):
+                    shift=ik*nao+1
+                    for i in range(nao):
+                        for j in range(i,nao):
+                            int_ij = intval_kpts_ao[ik,i,j]
+                            if abs(int_ij) > thresh:
+                                outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
+        if print_mo_ints_mono:
+            intval_kpts_mo = np.einsum('kim,kij,kjn->kmn',mo_k.conj(),intval_kpts_ao,mo_k)
+            with open('%s_mo_complex' % name,'w') as outfile:
+                for ik in range(Nk):
+                    shift=ik*nmo+1
+                    for i in range(nmo):
+                        for j in range(i,nmo):
+                            int_ij = intval_kpts_mo[ik,i,j]
+                            if abs(int_ij) > thresh:
+                                outfile.write('%s %s %s %s\n' % (i+shift, j+shift, int_ij.real, int_ij.imag))
+
+  
+    # ___                              _    
+    #  |  ._ _|_  _   _  ._ _. |  _   |_) o 
+    # _|_ | | |_ (/_ (_| | (_| | _>   |_) | 
+    #                 _|                    
+    #
+    kconserv = tools.get_kconserv(cell, kpts)
+
+    with open('kconserv_complex','w') as outfile:
+        for a in range(Nk):
+            for b in range(Nk):
+                for c in range(Nk):
+                    d = kconserv[a,b,c]
+                    outfile.write('%s %s %s %s\n' % (a+1,c+1,b+1,d+1))
+    
+
+    intfile=h5py.File(mf.with_df._cderi,'r')
+
+    j3c = intfile.get('j3c')
+    naosq = nao*nao
+    naotri = (nao*(nao+1))//2
+    j3ckeys = list(j3c.keys())
+    j3ckeys.sort(key=lambda strkey:int(strkey))
+
+    # in new(?) version of PySCF, there is an extra layer of groups before the datasets
+    # datasets used to be [/j3c/0,   /j3c/1,   /j3c/2,   ...]
+    # datasets now are    [/j3c/0/0, /j3c/1/0, /j3c/2/0, ...]
+    j3clist = [j3c.get(i+'/0') for i in j3ckeys]
+    if j3clist==[None]*len(j3clist):
+    # if using older version, stop before last level
+        j3clist = [j3c.get(i) for i in j3ckeys]
+
+    nkinvsq = 1./np.sqrt(Nk)
+
+    # dimensions are (kikj,iaux,jao,kao), where kikj is compound index of kpts i and j
+    # output dimensions should be reversed (nao, nao, naux, nkptpairs)
+    j3arr=np.array([(i.value.reshape([-1,nao,nao]) if (i.shape[1] == naosq) else makesq3(i.value,nao)) * nkinvsq for i in j3clist])
+
+    nkpt_pairs = j3arr.shape[0]
+
+    if print_ao_ints_df:
+        with open('df_ao_integral_array','w') as outfile:
+            pass
+        with open('df_ao_integral_array','a') as outfile:
+            for k,kpt_pair in enumerate(j3arr):
+                for iaux,dfbasfunc in enumerate(kpt_pair):
+                    for i,i0 in enumerate(dfbasfunc):
+                        for j,v in enumerate(i0):
+                            if (abs(v) > bielec_int_threshold):
+                                outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
+
+    if print_mo_ints_df:
+        kpair_list=[]
+        for i in range(Nk):
+            for j in range(Nk):
+                if(i>=j):
+                    kpair_list.append((i,j,idx2_tri((i,j))))
+        j3mo = np.array([np.einsum('mij,ik,jl->mkl',j3arr[kij],mo_k[ki].conj(),mo_k[kj]) for ki,kj,kij in kpair_list])
+        with open('df_mo_integral_array','w') as outfile:
+            pass
+        with open('df_mo_integral_array','a') as outfile:
+            for k,kpt_pair in enumerate(j3mo):
+                for iaux,dfbasfunc in enumerate(kpt_pair):
+                    for i,i0 in enumerate(dfbasfunc):
+                        for j,v in enumerate(i0):
+                            if (abs(v) > bielec_int_threshold):
+                                outfile.write('%s %s %s %s %s %s\n' % (i+1,j+1,iaux+1,k+1,v.real,v.imag))
+
+
+
+#    eri_4d_ao = np.zeros((Nk,nao,Nk,nao,Nk,nao,Nk,nao), dtype=np.complex)
+#    for d, kd in enumerate(kpts):
+#        for c, kc in enumerate(kpts):
+#            if c > d: break
+#            idx2_cd = idx2_tri(c,d)
+#            for b, kb in enumerate(kpts):
+#                if b > d: break
+#                a = kconserv[b,c,d]
+#                if idx2_tri(a,b) > idx2_cd: continue
+#                if ((c==d) and (a>b)): continue
+#                ka = kpts[a]
+#                v = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
+#                v *= 1./Nk
+#                eri_4d_ao[a,:,b,:,c,:,d] = v
+#    
+#    eri_4d_ao = eri_4d_ao.reshape([Nk*nao]*4)
+
+
+    if (print_ao_ints_bi or print_mo_ints_bi):
+        if print_ao_ints_bi:
+            with open('bielec_ao_complex','w') as outfile: 
+                pass
+        if print_mo_ints_bi:
+            with open('bielec_mo_complex','w') as outfile: 
+                pass
+        for d, kd in enumerate(kpts):
+            for c, kc in enumerate(kpts):
+                if c > d: break
+                idx2_cd = idx2_tri((c,d))
+                for b, kb in enumerate(kpts):
+                    if b > d: break
+                    a = kconserv[b,c,d]
+                    if idx2_tri((a,b)) > idx2_cd: continue
+                    if ((c==d) and (a>b)): continue
+                    ka = kpts[a]
+
+                    if print_ao_ints_bi:
+                        with open('bielec_ao_complex','a') as outfile:
+                            eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
+                            eri_4d_ao_kpt *= 1./Nk
+                            for l in range(nao):
+                                ll=l+d*nao
+                                for j in range(nao):
+                                    jj=j+c*nao
+                                    if jj>ll: break
+                                    idx2_jjll = idx2_tri((jj,ll))
+                                    for k in range(nao):
+                                        kk=k+b*nao
+                                        if kk>ll: break
+                                        for i in range(nao):
+                                            ii=i+a*nao
+                                            if idx2_tri((ii,kk)) > idx2_jjll: break
+                                            if ((jj==ll) and (ii>kk)): break
+                                            v=eri_4d_ao_kpt[i,k,j,l]
+                                            if (abs(v) > bielec_int_threshold):
+                                                outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
+            
+                    if print_mo_ints_bi:
+                        with open('bielec_mo_complex','a') as outfile:
+                            eri_4d_mo_kpt = mf.with_df.ao2mo([mo_k[a], mo_k[b], mo_k[c], mo_k[d]],
+                                                              [ka,kb,kc,kd],compact=False).reshape((nmo,)*4)
+                            eri_4d_mo_kpt *= 1./Nk
+                            for l in range(nmo):
+                                ll=l+d*nmo
+                                for j in range(nmo):
+                                    jj=j+c*nmo
+                                    if jj>ll: break
+                                    idx2_jjll = idx2_tri((jj,ll))
+                                    for k in range(nmo):
+                                        kk=k+b*nmo
+                                        if kk>ll: break
+                                        for i in range(nmo):
+                                            ii=i+a*nmo
+                                            if idx2_tri((ii,kk)) > idx2_jjll: break
+                                            if ((jj==ll) and (ii>kk)): break
+                                            v=eri_4d_mo_kpt[i,k,j,l]
+                                            if (abs(v) > bielec_int_threshold):
+                                                outfile.write('%s %s %s %s %s %s\n' % (ii+1,jj+1,kk+1,ll+1,v.real,v.imag))
+
+    
+  
+#def testpyscf2QP(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8):
+#    '''
+#    kpts = List of kpoints coordinates. Cannot be null, for gamma is other script
+#    kmesh = Mesh of kpoints (optional)
+#    cas_idx = List of active MOs. If not specified all MOs are actives
+#    int_threshold = The integral will be not printed in they are bellow that
+#    '''
+#
+#    from pyscf.pbc import ao2mo
+#    from pyscf.pbc import tools
+#    from pyscf.pbc.gto import ecp
+#
+#    mo_coef_threshold = int_threshold
+#    ovlp_threshold = int_threshold
+#    kin_threshold = int_threshold
+#    ne_threshold = int_threshold
+#    bielec_int_threshold = int_threshold
+#
+#    natom = len(cell.atom_coords())
+#    print('n_atom per kpt',   natom)
+#    print('num_elec per kpt', cell.nelectron)
+#
+#    mo_coeff = mf.mo_coeff
+#    # Mo_coeff actif
+#    mo_k = np.array([c[:,cas_idx] for c in mo_coeff] if cas_idx is not None else mo_coeff)
+#    e_k =  np.array([e[cas_idx] for e in mf.mo_energy] if cas_idx is not None else mf.mo_energy)
+#
+#    Nk, nao, nmo = mo_k.shape
+#    print("n Kpts", Nk)
+#    print("n active Mos per kpt", nmo)
+#    print("n AOs per kpt", nao)
+#
+#    naux = mf.with_df.get_naoaux()
+#    print("n df fitting functions", naux)
+#
+#    #                             _
+#    # |\ |      _ |  _   _. ._   |_)  _  ._      |  _ o  _  ._
+#    # | \| |_| (_ | (/_ (_| |    | \ (/_ |_) |_| | _> | (_) | |
+#    #                                    |
+#
+#    #Total energy shift due to Ewald probe charge = -1/2 * Nelec*madelung/cell.vol =
+#    shift = tools.pbc.madelung(cell, kpts)*cell.nelectron * -.5
+#    e_nuc = (cell.energy_nuc() + shift)*Nk
+#
+#    print('nucl_repul', e_nuc)
+#
+#
+#    # ___
+#    #  |  ._ _|_  _   _  ._ _. |  _   |\/|  _  ._   _  
+#    # _|_ | | |_ (/_ (_| | (_| | _>   |  | (_) | | (_) 
+#    #                 _|                              
+#   
+#    if mf.cell.pseudo:
+#        v_kpts_ao = np.reshape(mf.with_df.get_pp(kpts=kpts),(Nk,nao,nao))
+#    else:
+#        v_kpts_ao = np.reshape(mf.with_df.get_nuc(kpts=kpts),(Nk,nao,nao))
+#    if len(cell._ecpbas) > 0:
+#        v_kpts_ao += np.reshape(ecp.ecp_int(cell, kpts),(Nk,nao,nao))
+#
+#    ne_ao = ('ne',v_kpts_ao,ne_threshold)
+#    ovlp_ao = ('overlap',np.reshape(mf.get_ovlp(cell=cell,kpts=kpts),(Nk,nao,nao)),ovlp_threshold)
+#    kin_ao = ('kinetic',np.reshape(cell.pbc_intor('int1e_kin',1,1,kpts=kpts),(Nk,nao,nao)),kin_threshold)
+#
+#  
+#    # ___                              _    
+#    #  |  ._ _|_  _   _  ._ _. |  _   |_) o 
+#    # _|_ | | |_ (/_ (_| | (_| | _>   |_) | 
+#    #                 _|                    
+#    #
+#    kconserv = tools.get_kconserv(cell, kpts)
+#
+#    
+#    import h5py
+#
+#    intfile=h5py.File(mf.with_df._cderi,'r')
+#
+#    j3c = intfile.get('j3c')
+#    naosq = nao*nao
+#    naotri = (nao*(nao+1))//2
+#    j3keys = list(j3c.keys())
+#    j3keys.sort(key=lambda x:int(x))
+#    j3clist = [j3c.get(i) for i in j3keys]
+#    nkinvsq = 1./np.sqrt(Nk)
+#    
+#    # dimensions are (kikj,iaux,jao,kao), where kikj is compound index of kpts i and j
+#    # output dimensions should be reversed (nao, nao, naux, nkptpairs)
+#    j3arr=np.array([(pad(i.value.reshape([-1,nao,nao]),[naux,nao,nao]) if (i.shape[1] == naosq) else makesq(i.value,naux,nao)) * nkinvsq for i in j3clist])
+#
+#    nkpt_pairs = j3arr.shape[0]
+#
+#    kpair_list=[]
+#    for i in range(Nk):
+#        for j in range(Nk):
+#            if(i>=j):
+#                kpair_list.append((i,j,idx2_tri((i,j))))
+#    j3mo = np.array([np.einsum('mij,ik,jl->mkl',j3arr[kij,:,:,:],mo_k[ki,:,:].conj(),mo_k[kj,:,:]) for ki,kj,kij in kpair_list])
+#
+#
+#
+#    eri_mo = np.zeros(4*[nmo*Nk],dtype=np.complex128)
+#    eri_ao = np.zeros(4*[nao*Nk],dtype=np.complex128)
+#
+#    for d, kd in enumerate(kpts):
+#        for c, kc in enumerate(kpts):
+#            for b, kb in enumerate(kpts):
+#                a = kconserv[b,c,d]
+#                ka = kpts[a]
+#                eri_4d_ao_kpt = mf.with_df.get_ao_eri(kpts=[ka,kb,kc,kd],compact=False).reshape((nao,)*4)
+#                eri_4d_ao_kpt *= 1./Nk
+#                for l in range(nao):
+#                    ll=l+d*nao
+#                    for j in range(nao):
+#                        jj=j+c*nao
+#                        for k in range(nao):
+#                            kk=k+b*nao
+#                            for i in range(nao):
+#                                ii=i+a*nao
+#                                v=eri_4d_ao_kpt[i,k,j,l]
+#                                eri_ao[ii,kk,jj,ll]=v
+#        
+#                eri_4d_mo_kpt = mf.with_df.ao2mo([mo_k[a], mo_k[b], mo_k[c], mo_k[d]],
+#                                                  [ka,kb,kc,kd],compact=False).reshape((nmo,)*4)
+#                eri_4d_mo_kpt *= 1./Nk
+#                for l in range(nmo):
+#                    ll=l+d*nmo
+#                    for j in range(nmo):
+#                        jj=j+c*nmo
+#                        for k in range(nmo):
+#                            kk=k+b*nmo
+#                            for i in range(nmo):
+#                                ii=i+a*nmo
+#                                v=eri_4d_mo_kpt[i,k,j,l]
+#                                eri_mo[ii,kk,jj,ll]=v
+#        
+#    return (mo_k,j3arr,j3mo,eri_ao,eri_mo,kpair_list)
+    
+  
diff --git a/src/utils_periodic/create_ezfio_complex_3idx.py b/src/utils_periodic/create_ezfio_complex_3idx.py
new file mode 100755
index 00000000..d0bacf81
--- /dev/null
+++ b/src/utils_periodic/create_ezfio_complex_3idx.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+from ezfio import ezfio
+import h5py
+
+import sys
+filename = sys.argv[1]
+h5filename = sys.argv[2]
+#num_elec, nucl_num, mo_num = map(int,sys.argv[2:5])
+
+#nuclear_repulsion = float(sys.argv[5])
+#ao_num = int(sys.argv[6])
+#n_kpts = int(sys.argv[7])
+#n_aux = int(sys.argv[8])
+ezfio.set_file(filename)
+qph5=h5py.File(h5filename.'r')
+
+ezfio.electrons_elec_alpha_num = qph5['electrons'].attrs['elec_alpha_num']
+ezfio.electrons_elec_beta_num = qph5['electrons'].attrs['elec_beta_num']
+
+nucl_num = qph5['nuclei'].attrs['nucl_num']
+kpt_num = qph5['nuclei'].attrs['kpt_num']
+#df_num = qph5['???'].attrs['df_num']
+
+mo_num = qph5['mo_basis'].attrs['mo_num']
+ezfio.set_mo_basis_mo_num(mo_num)
+
+#ao_num = mo_num
+#Important !
+import math
+nelec_per_kpt = num_elec // n_kpts
+nelec_alpha_per_kpt = int(math.ceil(nelec_per_kpt / 2.))
+nelec_beta_per_kpt = int(math.floor(nelec_per_kpt / 2.))
+
+ezfio.electrons_elec_alpha_num = int(nelec_alpha_per_kpt * n_kpts)
+ezfio.electrons_elec_beta_num = int(nelec_beta_per_kpt * n_kpts)
+
+#ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
+#ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
+
+#ezfio.set_utils_num_kpts(n_kpts)
+#ezfio.set_integrals_bielec_df_num(n_aux)
+
+#Important
+ezfio.set_nuclei_nucl_num(nucl_num)
+ezfio.set_nuclei_nucl_charge([0.]*nucl_num)
+ezfio.set_nuclei_nucl_coord( [ [0.], [0.], [0.] ]*nucl_num )
+ezfio.set_nuclei_nucl_label( ['He'] * nucl_num )
+
+ezfio.set_nuclei_io_nuclear_repulsion('Read')
+ezfio.set_nuclei_nuclear_repulsion(nuclear_repulsion)
+
+# Ao num
+#ao_num = mo_num
+ezfio.set_ao_basis_ao_basis("Dummy one. We read MO")
+ezfio.set_ao_basis_ao_num(ao_num)
+ezfio.set_ao_basis_ao_nucl([1]*ao_num) #Maybe put a realy incorrect stuff
+
+#Just need one
+ao_prim_num_max  = 5
+
+d = [ [0] *ao_prim_num_max]*ao_num
+ezfio.set_ao_basis_ao_prim_num([ao_prim_num_max]*ao_num)
+ezfio.set_ao_basis_ao_power(d)
+ezfio.set_ao_basis_ao_coef(d)
+ezfio.set_ao_basis_ao_expo(d)
+
+#Dummy one
+ao_md5 = '3b8b464dfc95f282129bde3efef3c502'
+ezfio.set_ao_basis_ao_md5(ao_md5)
+ezfio.set_mo_basis_ao_md5(ao_md5)
+
+
+ezfio.set_mo_basis_mo_num(mo_num)
+c_mo = [[1 if i==j else 0 for i in range(mo_num)] for j in range(ao_num)]
+ezfio.set_mo_basis_mo_coef([ [0]*mo_num] * ao_num)
+#ezfio.set_mo_basis_mo_coef_real(c_mo)
+
+ezfio.set_nuclei_is_periodic(True)
diff --git a/src/utils_periodic/qp2-pbc-diff.txt b/src/utils_periodic/qp2-pbc-diff.txt
index 9bde53fc..d1953c8b 100644
--- a/src/utils_periodic/qp2-pbc-diff.txt
+++ b/src/utils_periodic/qp2-pbc-diff.txt
@@ -73,16 +73,24 @@ later:
 
 
 NOTES:
-  number of unique 4-tuples with 8-fold symmetry is a8(n)=n*(n+1)*(n^2+n+2)/8
-  number of unique 4-tuples with 4-fold symmetry is a4(n)=n^2*(n^2+3)/4
-  a8 is number of unique real 2e ints with n mos
-  a4 is number of unique* complex 2e ints with n mos (where p+i*q and p-i*q are counted as one, not two)
-  a4(n) = a8(n) + a8(n-1)
-  
-  we can already generate the list of <ij|kl> with unique values for the 8-fold case
-  the set of these for 4-fold symmetry is the union of the 8-fold set for n and the 8-fold set for n-1 with a simple transformation
-  <ij|kl>_{4,n} = <ij|kl>_{8,n} + <(k+1)j|i(l+1)>_{8,n-1}
+  3-index integrals
+    <ij|kl> = \sum_\mu (ik|\mu)(jl|\mu)
+    store (ik|\mu) for I<=K
+    if i>k, take conjugate transpose in first two dimensions
 
+    df_mo(:,:,mu,kjkl) = C(:,:,kj)^\dagger.df_ao(:,:,mu,kjkl).C(:,:,kl)
+
+  2e int compound indexing
+    number of unique 4-tuples with 8-fold symmetry is a8(n)=n*(n+1)*(n^2+n+2)/8
+    number of unique 4-tuples with 4-fold symmetry is a4(n)=n^2*(n^2+3)/4
+    a8 is number of unique real 2e ints with n mos
+    a4 is number of unique* complex 2e ints with n mos (where p+i*q and p-i*q are counted as one, not two)
+    a4(n) = a8(n) + a8(n-1)
+    
+    we can already generate the list of <ij|kl> with unique values for the 8-fold case
+    the set of these for 4-fold symmetry is the union of the 8-fold set for n and the 8-fold set for n-1 with a simple transformation
+    <ij|kl>_{4,n} = <ij|kl>_{8,n} + <(k+1)j|i(l+1)>_{8,n-1}
+