diff --git a/bin/qp_convert_h5_to_ezfio b/bin/qp_convert_h5_to_ezfio
index 5fb5554b..5c8f8ad3 100755
--- a/bin/qp_convert_h5_to_ezfio
+++ b/bin/qp_convert_h5_to_ezfio
@@ -21,7 +21,93 @@ from docopt import docopt
 import gzip
 #fname = sys.argv[1]
 #qph5name = sys.argv[2]
+def idx2_tri(i,j):
+    """
+    for 0-indexed counting
+    """
+    p = max(i,j)
+    q = min(i,j)
+    return q + (p*(p+1))//2
+def idx2_tri_1(i,j):
+    """
+    for 1-indexed counting
+    """
+    p = max(i,j)
+    q = min(i,j)
+    return q + (p*(p-1))//2
 
+def idx4_cplx_1(i,j,k,l):
+    """
+    original function from qp2 (fortran counting)
+    """
+    p = idx2_tri_1(i,k)
+    q = idx2_tri_1(j,l)
+    i1 = idx2_tri_1(p,q)
+    return (i1,p,q)
+
+def ao_idx_map_sign(i,j,k,l):
+    """
+    qp2 indexing
+    """
+    idx,ik,jl = idx4_cplx_1(i,j,k,l)
+    ij = idx2_tri_1(i,j)
+    kl = idx2_tri_1(k,l)
+    idx = 2*idx - 1
+
+    if (ij==kl):
+        sign = 0.0
+        use_map1 = False
+    else:
+        if ik==jl:
+            if i<k:
+                sign = 1.0
+                use_map1 = True
+            else:
+                sign = -1.0
+                use_map1 = True
+        elif i==k:
+            if j<l:
+                sign = 1.0
+                use_map1 = True
+            else:
+                sign = -1.0
+                use_map1 = True
+        elif j==l:
+            if i<k:
+                sign = 1.0
+                use_map1 = True
+            else:
+                sign = -1.0
+                use_map1 = True
+        elif ((i<k) == (j<l)):
+            if i<k:
+                sign = 1.0
+                use_map1 = True
+            else:
+                sign = -1.0
+                use_map1 = True
+        else:
+            if ((j<l) == (ik<jl)):
+                sign = 1.0
+                use_map1 = False
+            else:
+                sign = -1.0
+                use_map1 = False
+    return (idx, use_map1, sign)
+def get_ao_int_cplx(i,j,k,l,map1,map2):
+    idx,use_m1,sgn = ao_idx_map_sign(i,j,k,l)
+    if use_m1:
+        tmp_re = map1[idx]
+        tmp_im = map1[idx+1]
+        tmp_im *= sgn
+    else:
+        tmp_re = map2[idx]
+        if sgn != 0.0:
+            tmp_im = map2[idx+1]
+            tmp_im *= sgn
+        else:
+            tmp_im = 0.0
+    return tmp_re + 1j*tmp_im
 def kconserv_p_from_qkk2_mk(qkk2,mk):
     nk, nk2 = qkk2.shape
     assert(nk == nk2)
@@ -30,6 +116,7 @@ def kconserv_p_from_qkk2_mk(qkk2,mk):
         for j in range(nk):
             for k in range(nk):
                 kcon_p[i,j,k] = qkk2[mk[j],qkk2[k,i]]
+                assert(qkk2[mk[j],qkk2[k,i]] == qkk2[qkk2[j,k],i])
     return kcon_p
 
 
@@ -201,6 +288,11 @@ def convert_mol(filename,qph5path):
 def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
     import json
     from scipy.linalg import block_diag
+    dump_fci, dump_cd = (False, False)
+    #dump_fci, dump_cd = (True, False)
+    #dump_fci, dump_cd = (False, True)
+    #dump_fci2 = True
+    dump_fci2 = False
 
     ezfio.set_file(filename)
     ezfio.set_nuclei_is_complex(True)
@@ -385,10 +477,10 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
             isparse += 1
         else:
             kpt_sparse_map[i] = -kpt_sparse_map[minusk[i]-1]
-    ezfio.set_nuclei_kconserv(kconserv)
+    ezfio.set_nuclei_kconserv(kconserv.transpose(2,1,0))
     ezfio.set_nuclei_io_kconserv('Read')
     ezfio.set_nuclei_minusk(minusk)
-    ezfio.set_nuclei_qktok2(QKTok2)
+    ezfio.set_nuclei_qktok2(np.transpose(QKTok2)) #transposed for correct col-major ordering
     ezfio.set_nuclei_kpt_sparse_map(kpt_sparse_map)
     ezfio.set_nuclei_unique_kpt_num(unique_kpt_num)
     # kpt_sparse_map
@@ -445,7 +537,263 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
             #ao_chol_two_e_ints = np.vstack(L_list)
             #ao_chol_two_e_ints = ao_chol_two_e_ints.transpose()
             #TODO: check dims/reshape/transpose
-            ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all)
+            #ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all.transpose((5,2,3,4,1,0)))
+            ezfio.set_ao_two_e_ints_chol_ao_integrals_complex(L_all.transpose((0,1,4,3,2,5)))
+
+            def fortformat(x0):
+                x = f'{abs(x0):25.14E}'.strip()
+                xsign = '-' if (x0<0) else ''
+                e = x.find('E')
+                return xsign + f'0.{x[0]}{x[2:e]}{x[e:e+2]}{abs(int(x[e+1:])*1+1):02d}'
+            if dump_cd:
+                #for qi in range(unique_kpt_num):
+                #    for ki in range(kpt_num):
+                #        for i in range(ao_num_per_kpt):
+                #            for j in range(ao_num_per_kpt):
+                #                for ci in range(nchol_per_kpt_max):
+                #                    vr = L_all[qi,ki,i,j,ci,0]
+                #                    vi = L_all[qi,ki,i,j,ci,1]
+                #                    print(f'{qi:6d} {ki:6d} {i:6d} {j:6d} {ci:6d} {vr:25.15E} {vi:25.15E}')
+                Lnew = L_all.transpose((5,2,3,4,1,0))
+                #for qi in range(unique_kpt_num):
+                #    for ki in range(kpt_num):
+                #        for ci in range(nchol_per_kpt_max):
+                #            for j in range(ao_num_per_kpt):
+                #                for i in range(ao_num_per_kpt):
+                #                    vr = Lnew[0,i,j,ci,ki,qi]
+                #                    vi = Lnew[1,i,j,ci,ki,qi]
+                #                    print(f'{qi:6d} {ki:6d} {ci:6d} {j:6d} {i:6d} {vr:25.15E} {vi:25.15E}')
+                for qi in range(unique_kpt_num):
+                    for ki in range(kpt_num):
+                        for ci in range(nchol_per_kpt_max):
+                            for i in range(ao_num_per_kpt):
+                                for j in range(ao_num_per_kpt):
+                                    #vr = Lnew[0,i,j,ci,ki,qi]
+                                    #vi = Lnew[1,i,j,ci,ki,qi]
+                                    vr = L_all[qi,ki,i,j,ci,0]
+                                    vi = L_all[qi,ki,i,j,ci,1]
+                                    print(f'{qi+1:6d} {ki+1:6d} {ci+1:6d} {i+1:6d} {j+1:6d} {fortformat(vr):>25s} {fortformat(vi):>25s}')
+
+            if dump_fci:
+                Wfull = np.zeros((ao_num_tot, ao_num_tot, ao_num_tot, ao_num_tot), dtype=np.complex128)
+                for Qi in range(kpt_num):
+                    Qloc = abs(kpt_sparse_map[Qi])-1
+                    Qneg = (kpt_sparse_map[Qi] < 0)
+                    LQ00 = L_all[Qloc]
+                    #LQ0a = LQ00.view(dtype=np.complex128)
+                    #print(f'LQ0a.shape {LQ0a.shape}')
+                    #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1))
+                    #print(f'LQ0a1.shape {LQ0a1.shape}')
+                    LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1]
+                    #print(f'LQ0.shape {LQ0.shape}')
+                    #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}')
+
+                    
+
+                    for kp in range(kpt_num):
+                        kr = QKTok2[Qi,kp]-1
+                        for ks in range(kpt_num):
+                            kq = QKTok2[Qi,ks]-1
+                            if Qneg:
+                                A = LQ0[kr].transpose((1,0,2)).conj()
+                                B = LQ0[kq]
+                                W = np.einsum('prn,sqn->pqrs',A,B)
+                            else:
+                                A = LQ0[kp]
+                                B = LQ0[ks].transpose((1,0,2)).conj()
+                                W = np.einsum('rpn,qsn->pqrs',A,B)
+                            p0 = kp*ao_num_per_kpt
+                            r0 = kr*ao_num_per_kpt
+                            q0 = kq*ao_num_per_kpt
+                            s0 = ks*ao_num_per_kpt
+                            for ip in range(ao_num_per_kpt):
+                                for iq in range(ao_num_per_kpt):
+                                    for ir in range(ao_num_per_kpt):
+                                        for i_s in range(ao_num_per_kpt):
+                                            v = W[ip,iq,ir,i_s]
+                                            print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}')
+                            Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy()
+                H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128)
+                for Qi in range(kpt_num):
+                    hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()]
+                    hi = hi0[:,:,0] + 1j*hi0[:,:,1]
+                    H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy()
+                mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num
+                E1=0
+                E2j=0
+                E2k=0
+                for imo, iocc in enumerate(mo_occ):
+                    if iocc:
+                        E1 += 2*H1[imo,imo]
+                        for jmo, jocc in enumerate(mo_occ):
+                            if jocc:
+                                E2j += 2* Wfull[imo,jmo,imo,jmo]
+                                E2k -= Wfull[imo,jmo,jmo,imo]
+                print(f'E1  = {E1:25.15E}')
+                print(f'E2j = {E2j:25.15E}')
+                print(f'E2k = {E2k:25.15E}')
+                print(f'E2  = {E2j+E2k:25.15E}')
+
+
+            if dump_fci2:
+                ao_map1 = {}
+                ao_map2 = {}
+                ao_map1_idx = {}
+                ao_map2_idx = {}
+                Wdict = {}
+                for kQ in range(kpt_num):
+                    for kl in range(kpt_num):
+                        kj = QKTok2[kQ,kl]-1
+                        if (kj>kl):
+                            continue
+                        kjkl2 = idx2_tri(kj,kl)
+                        Qneg = (kpt_sparse_map[kQ] < 0)
+                        Qloc = abs(kpt_sparse_map[kQ]) - 1
+                        if not Qneg:
+                            ints_jl0 = L_all[Qloc,kl,:,:,:,:]
+                            ints_jl = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1]
+                        else:
+                            ints_jl0 = L_all[Qloc,kj,:,:,:,:]
+                            ints_jl1 = ints_jl0[:,:,:,0]+1j*ints_jl0[:,:,:,1]
+                            ints_jl = ints_jl1.transpose((1,0,2)).conj()
+                        for kk in range(kl+1):
+                            ki = QKTok2[minusk[kk]-1,kQ]-1 #TODO: check
+                            if ki != kconserv[kl,kk,kj]-1:
+                                print(ki,kconserv[kl,kk,kj],kl,kk,kj)
+                                assert( ki == kconserv[kl,kk,kj]-1) #TODO: check
+                            if (ki > kl):
+                                continue
+                            kikk2 = idx2_tri(ki,kk)
+                            if not Qneg:
+                                ints_ik0 = L_all[Qloc,ki,:,:,:]
+                                ints_ik = ints_ik0[:,:,:,0] + 1j*ints_ik0[:,:,:,1]
+                            else:
+                                ints_ik0 = L_all[Qloc,kk,:,:,:]
+                                ints_ik1 = ints_ik0[:,:,:,0]+1j*ints_ik0[:,:,:,1]
+                                ints_ik = ints_ik1.transpose((1,0,2)).conj()
+                            ints_jl_flat = ints_jl.reshape((ao_num_per_kpt**2, -1))
+                            ints_ik_flat = ints_ik.reshape((ao_num_per_kpt**2, -1))
+
+                            ints_ikjl = np.einsum('an,bn->ab',ints_ik_flat,ints_jl_flat).reshape((ao_num_per_kpt,)*4)
+
+                            for il in range(ao_num_per_kpt):
+                                l = il + kl*ao_num_per_kpt
+                                for ij in range(ao_num_per_kpt):
+                                    j = ij + kj*ao_num_per_kpt
+                                    if (j>l):
+                                        break
+                                    jl2 = idx2_tri(j,l)
+                                    for ik in range(ao_num_per_kpt):
+                                        k = ik + kk*ao_num_per_kpt
+                                        if (k>l):
+                                            break
+                                        for ii in range(ao_num_per_kpt):
+                                            i = ii + ki*ao_num_per_kpt
+                                            if ((j==l) and (i>k)):
+                                                break
+                                            ik2 = idx2_tri(i,k)
+                                            if (ik2 > jl2):
+                                                break
+                                            integral = ints_ikjl[ii,ik,ij,il]
+                                            if abs(integral) < 1E-15:
+                                                continue
+                                            idx_tmp,use_map1,sign = ao_idx_map_sign(i+1,j+1,k+1,l+1)
+                                            tmp_re = integral.real
+                                            tmp_im = integral.imag
+                                            Wdict[i,j,k,l] = tmp_re + 1j*tmp_im
+                                            if use_map1:
+                                                if idx_tmp in ao_map1:
+                                                    print(idx_tmp,1)
+                                                    raise
+                                                ao_map1[idx_tmp] = tmp_re
+                                                ao_map1_idx[idx_tmp] = (i,j,k,l,'re')
+                                                if sign != 0.0:
+                                                    ao_map1[idx_tmp+1] = tmp_im*sign
+                                                    ao_map1_idx[idx_tmp+1] = (i,j,k,l,'im')
+                                            else:
+                                                if idx_tmp in ao_map2:
+                                                    print(idx_tmp,2)
+                                                    raise
+                                                ao_map2[idx_tmp] = tmp_re
+                                                ao_map2_idx[idx_tmp] = (i,j,k,l,'re')
+                                                if sign != 0.0:
+                                                    ao_map2[idx_tmp+1] = tmp_im*sign
+                                                    ao_map2_idx[idx_tmp+1] = (i,j,k,l,'im')
+#                for idx in ao_map1:
+#                    i,j,k,l,ax = ao_map1_idx[idx]
+#                    print(f'1,{idx},{i},{j},{k},{l},{ax},{ao_map1[idx]}')
+#                for idx in ao_map2:
+#                    i,j,k,l,ax = ao_map2_idx[idx]
+#                    print(f'2,{idx},{i},{j},{k},{l},{ax},{ao_map2[idx]}')
+                for idx in Wdict:
+                    i,j,k,l = idx
+                    v = Wdict[idx]
+                    print(f'{i+1:6d} {j+1:6d} {k+1:6d} {l+1:6d} {fortformat(v.real):>25s} {fortformat(v.imag):>25s}')
+
+
+                                            
+
+
+                #for Qi in range(kpt_num):
+                #    Qloc = abs(kpt_sparse_map[Qi])-1
+                #    Qneg = (kpt_sparse_map[Qi] < 0)
+                #    LQ00 = L_all[Qloc]
+                #    #LQ0a = LQ00.view(dtype=np.complex128)
+                #    #print(f'LQ0a.shape {LQ0a.shape}')
+                #    #LQ0a1 = LQ0a.reshape((kpt_num,ao_num_per_kpt,ao_num_per_kpt,-1))
+                #    #print(f'LQ0a1.shape {LQ0a1.shape}')
+                #    LQ0 = LQ00[:,:,:,:,0] + 1j*LQ00[:,:,:,:,1]
+                #    #print(f'LQ0.shape {LQ0.shape}')
+                #    #print(f'abdiff {np.abs(LQ0a1 - LQ0).max()}')
+
+                #    
+
+                #    for kp in range(kpt_num):
+                #        kr = QKTok2[Qi,kp]-1
+                #        for ks in range(kpt_num):
+                #            kq = QKTok2[Qi,ks]-1
+                #            if Qneg:
+                #                A = LQ0[kr].transpose((1,0,2)).conj()
+                #                B = LQ0[kq]
+                #                W = np.einsum('prn,sqn->pqrs',A,B)
+                #            else:
+                #                A = LQ0[kp]
+                #                B = LQ0[ks].transpose((1,0,2)).conj()
+                #                W = np.einsum('rpn,qsn->pqrs',A,B)
+                #            p0 = kp*ao_num_per_kpt
+                #            r0 = kr*ao_num_per_kpt
+                #            q0 = kq*ao_num_per_kpt
+                #            s0 = ks*ao_num_per_kpt
+                #            for ip in range(ao_num_per_kpt):
+                #                for iq in range(ao_num_per_kpt):
+                #                    for ir in range(ao_num_per_kpt):
+                #                        for i_s in range(ao_num_per_kpt):
+                #                            v = W[ip,iq,ir,i_s]
+                #                            print(f'{p0+ip:5d} {q0+iq:5d} {r0+ir:5d} {s0+i_s:5d} {v.real:25.15E} {v.imag:25.15E}')
+                #            Wfull[p0:p0+ao_num_per_kpt,q0:q0+ao_num_per_kpt,r0:r0+ao_num_per_kpt,s0:s0+ao_num_per_kpt] = W.copy()
+                #H1 = np.zeros((ao_num_tot, ao_num_tot), dtype=np.complex128)
+                #for Qi in range(kpt_num):
+                #    hi0 = qph5[f'Hamiltonian/H1_kp{Qi}'][()]
+                #    hi = hi0[:,:,0] + 1j*hi0[:,:,1]
+                #    H1[Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt,Qi*ao_num_per_kpt:(Qi+1)*ao_num_per_kpt] = hi.copy()
+                #mo_occ = ([1,]* elec_alpha_num_per_kpt + [0,] * (ao_num_per_kpt - elec_alpha_num_per_kpt)) * kpt_num
+                #E1=0
+                #E2j=0
+                #E2k=0
+                #for imo, iocc in enumerate(mo_occ):
+                #    if iocc:
+                #        E1 += 2*H1[imo,imo]
+                #        for jmo, jocc in enumerate(mo_occ):
+                #            if jocc:
+                #                E2j += 2* Wfull[imo,jmo,imo,jmo]
+                #                E2k -= Wfull[imo,jmo,jmo,imo]
+                #print(f'E1  = {E1:25.15E}')
+                #print(f'E2j = {E2j:25.15E}')
+                #print(f'E2k = {E2k:25.15E}')
+                #print(f'E2  = {E2j+E2k:25.15E}')
+
+
+                            
 
             #(2,ao_basis.ao_num_per_kpt,ao_basis.ao_num_per_kpt,ao_two_e_ints.chol_num_max,nuclei.kpt_num,nuclei.unique_kpt_num#)
             """
@@ -991,6 +1339,12 @@ def convert_cplx(filename,qph5path):
 if __name__ == '__main__':
     ARGUMENTS = docopt(__doc__)
 
+    #for i in range(1,6):
+    #    for j in range(1,6):
+    #        for k in range(1,6):
+    #            for l in range(1,6):
+    #                idx,usem1,sgn = ao_idx_map_sign(i,j,k,l)
+    #                print(f'{i:4d} {j:4d} {k:4d} {l:4d}      {str(usem1)[0]:s}  {idx:6d}     {sgn:5.1f}')
     FILE = get_full_path(ARGUMENTS['FILE'])
     qmcpack = True
     rmg = False
diff --git a/src/ao_two_e_ints/cd_ao_ints.irp.f b/src/ao_two_e_ints/cd_ao_ints.irp.f
index 19305d42..e1888bf9 100644
--- a/src/ao_two_e_ints/cd_ao_ints.irp.f
+++ b/src/ao_two_e_ints/cd_ao_ints.irp.f
@@ -92,11 +92,13 @@ subroutine ao_map_fill_from_chol
       !TODO: verify the kj, kl as 4th index in expressions below
       if (kpt_sparse_map(kQ) > 0) then
         ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))
+        !ints_jl = dconjg(chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ)))
       else
         do i_ao=1,ao_num_per_kpt
           do j_ao=1,ao_num_per_kpt
             do i_cd=1,chol_num_max
               ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)))
+              !ints_jl(i_ao,j_ao,i_cd) = chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ))
             enddo
           enddo
         enddo
diff --git a/src/mo_two_e_ints/cd_mo_ints.irp.f b/src/mo_two_e_ints/cd_mo_ints.irp.f
index 19f6ccb3..1afc2729 100644
--- a/src/mo_two_e_ints/cd_mo_ints.irp.f
+++ b/src/mo_two_e_ints/cd_mo_ints.irp.f
@@ -81,7 +81,8 @@ subroutine mo_map_fill_from_chol_dot
         do i_mo=1,mo_num_per_kpt
           do j_mo=1,mo_num_per_kpt
             do i_cd=1,chol_num(kQ)
-              ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx)
+              !ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx)
+              ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(i_mo,j_mo,i_cd,kl,Q_idx))
             enddo
           enddo
         enddo
@@ -89,7 +90,8 @@ subroutine mo_map_fill_from_chol_dot
         do i_mo=1,mo_num_per_kpt
           do j_mo=1,mo_num_per_kpt
             do i_cd=1,chol_num(kQ)
-              ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx))
+              !ints_jl(i_cd,i_mo,j_mo) = dconjg(chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx))
+              ints_jl(i_cd,i_mo,j_mo) = chol_mo_integrals_complex(j_mo,i_mo,i_cd,kj,-Q_idx)
             enddo
           enddo
         enddo
@@ -268,7 +270,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k)
   complex*16,allocatable :: coef_i(:,:), coef_k(:,:), ints_ik(:,:), ints_tmp(:,:)
   double precision :: wall_1,wall_2,cpu_1,cpu_2
 
-  print*,'providing 3-index MO integrals from 3-index AO integrals'
+  print*,'providing 3-index CD MO integrals from 3-index CD AO integrals'
   
   cd_mo = 0.d0
 
@@ -315,7 +317,7 @@ subroutine chol_mo_from_chol_ao(cd_mo,cd_ao,n_mo,n_ao,n_cd,n_k,n_unique_k)
           )
   call wall_time(wall_2)
   call cpu_time(cpu_2)
-  print*,' 3-idx MO provided'
+  print*,' 3-idx CD MO provided'
   print*,'  cpu  time:',cpu_2-cpu_1,'s'
   print*,'  wall time:',wall_2-wall_1,'s  ( x ',(cpu_2-cpu_1)/(wall_2-wall_1),')'
 
diff --git a/src/utils_complex/dump_cd_ints.irp.f b/src/utils_complex/dump_cd_ints.irp.f
new file mode 100644
index 00000000..40adb2b4
--- /dev/null
+++ b/src/utils_complex/dump_cd_ints.irp.f
@@ -0,0 +1,29 @@
+program dump_cd_ksym
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+
+  integer ::q,k,n,i,j
+  double precision :: vr, vi
+  complex*16 :: v
+  print*,"chol_ao_integrals_complex q,k,n,i,j"
+  provide chol_ao_integrals_complex
+  do q = 1, unique_kpt_num
+    do k = 1, kpt_num
+      do n = 1, chol_num_max
+        do i = 1, ao_num_per_kpt
+          do j = 1, ao_num_per_kpt
+            v = chol_ao_integrals_complex(i,j,n,k,q)
+            vr = dble(v)
+            vi = dimag(v)
+            print '(5(I6,X),2(E25.15,X))', q, k, n, i, j, vr, vi
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end
diff --git a/src/utils_complex/dump_cd_ksym.irp.f b/src/utils_complex/dump_cd_ksym.irp.f
new file mode 100644
index 00000000..d5f9cfdb
--- /dev/null
+++ b/src/utils_complex/dump_cd_ksym.irp.f
@@ -0,0 +1,47 @@
+program dump_cd_ksym
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+
+  integer ::i,j,k,l
+  integer(key_kind) :: idx
+  logical :: use_map1
+  double precision :: sign
+  do i=1,5
+    do j=1,5
+      do k=1,5
+        do l=1,5
+          call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx,sign)
+          print'(4(I4,X),(L6),(I8),(F10.1))',i,j,k,l,use_map1,idx,sign
+        enddo
+      enddo
+    enddo
+  enddo
+
+  provide qktok2 minusk kconserv
+  print*,'minusk'
+  do i=1,kpt_num
+    j = minusk(i)
+    print'(2(I4))',i,j
+  enddo
+  print*,'qktok2'
+  do i=1,kpt_num
+    do j=1,kpt_num
+      k = qktok2(i,j)
+      print'(3(I4))',i,j,k
+    enddo
+  enddo
+  print*,'kconserv'
+  do i=1,kpt_num
+    do j=1,kpt_num
+      do k=1,kpt_num
+        l = kconserv(i,j,k)
+        print'(4(I4))',i,j,k,l
+      enddo
+    enddo
+  enddo
+
+end
diff --git a/src/utils_complex/test_cd_ksym.irp.f b/src/utils_complex/test_cd_ksym.irp.f
new file mode 100644
index 00000000..be1433be
--- /dev/null
+++ b/src/utils_complex/test_cd_ksym.irp.f
@@ -0,0 +1,234 @@
+program test_cd_ksym
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+
+  !integer ::i,j,k,l
+
+  provide qktok2 minusk kconserv
+  !print*,'minusk'
+  !do i=1,kpt_num
+  !  j = minusk(i)
+  !  print'(2(I4))',i,j
+  !enddo
+  !print*,'qktok2'
+  !do i=1,kpt_num
+  !  do j=1,kpt_num
+  !    k = qktok2(i,j)
+  !    print'(3(I4))',i,j,k
+  !  enddo
+  !enddo
+  !print*,'kconserv'
+  !do i=1,kpt_num
+  !  do j=1,kpt_num
+  !    do k=1,kpt_num
+  !      l = kconserv(i,j,k)
+  !      print'(4(I4))',i,j,k,l
+  !    enddo
+  !  enddo
+  !enddo
+
+  integer :: i,k,j,l
+  integer :: ki,kk,kj,kl
+  integer :: ii,ik,ij,il
+  integer :: kikk2,kjkl2,jl2,ik2
+  integer :: i_ao,j_ao,i_cd,kq
+
+  complex*16,allocatable :: ints_ik(:,:,:), ints_jl(:,:,:), ints_ikjl(:,:,:,:)
+
+  complex*16 :: integral
+  integer                        :: n_integrals_1, n_integrals_2
+  integer                        :: size_buffer
+  integer(key_kind),allocatable  :: buffer_i_1(:), buffer_i_2(:)
+  real(integral_kind),allocatable :: buffer_values_1(:), buffer_values_2(:)
+  double precision               :: tmp_re,tmp_im
+  integer                        :: ao_num_kpt_2
+
+  double precision               :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+  double precision               :: map_mb
+
+  logical :: use_map1
+  integer(keY_kind) :: idx_tmp
+  double precision :: sign
+
+  ao_num_kpt_2 = ao_num_per_kpt * ao_num_per_kpt
+
+  size_buffer = min(ao_num_per_kpt*ao_num_per_kpt*ao_num_per_kpt,16000000)
+  print*, 'Providing the ao_bielec integrals from 3-index cholesky integrals'
+  call write_time(6)
+!  call ezfio_set_integrals_bielec_disk_access_mo_integrals('Write')
+!  TOUCH read_mo_integrals read_ao_integrals write_mo_integrals write_ao_integrals 
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  !allocate( ints_jl(ao_num_per_kpt,ao_num_per_kpt,chol_num_max))
+
+  wall_0 = wall_1
+  ! ki + kj == kk + kl required for <ij|kl> to be nonzero
+  !TODO: change loops so that we only iterate over "correct" slices (i.e. ik block is stored directly, not as conj. transp.)
+  !      possible cases for (ik,jl) are (+,+), (+,-), (-,+), (-,-)
+  !      where + is the slice used as stored, and - is the conj. transp. of the stored data
+  !      (+,+) and (-,-) give the same information; we should always use (+,+)
+  !      (+,-) and (-,+) give the same information; we should always use (+,-)
+  do kQ = 1, kpt_num
+    do kl = 1, kpt_num
+      kj = qktok2(kQ,kl)
+      assert(kQ == qktok2(kj,kl))
+      if (kj>kl) cycle
+      call idx2_tri_int(kj,kl,kjkl2)
+      !TODO: verify the kj, kl as 4th index in expressions below
+      !if (kpt_sparse_map(kQ) > 0) then
+      !  ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))
+      !else
+      !  !do i_ao=1,ao_num_per_kpt
+      !  !  do j_ao=1,ao_num_per_kpt
+      !  !    do i_cd=1,chol_num_max
+      !  !      ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kj,-kpt_sparse_map(kQ)))
+      !  !    enddo
+      !  !  enddo
+      !  !enddo
+      !endif
+
+  !allocate( &
+  !  ints_ik(ao_num_per_kpt,ao_num_per_kpt,chol_num_max), &
+  !  ints_ikjl(ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt), &
+  !  buffer_i_1(size_buffer), &
+  !  buffer_i_2(size_buffer), &
+  !  buffer_values_1(size_buffer), &
+  !  buffer_values_2(size_buffer) &
+  !)
+
+      do kk=1,kl
+        ki = qktok2(minusk(kk),kQ)
+        assert(ki == kconserv(kl,kk,kj))
+        if (ki>kl) cycle
+      !  if ((kl == kj) .and. (ki > kk)) cycle
+        call idx2_tri_int(ki,kk,kikk2)
+        print*,kQ,kl,kj,kk,ki
+      !  if (kikk2 > kjkl2) cycle
+      !TODO: check this! (ki, kk slice index and transpose/notranspose)
+        !if (kpt_sparse_map(kQ) > 0) then
+        !  ints_ik = chol_ao_integrals_complex(:,:,:,ki,kpt_sparse_map(kQ))
+        !else
+        !  do i_ao=1,ao_num_per_kpt
+        !    do j_ao=1,ao_num_per_kpt
+        !      do i_cd=1,chol_num_max
+        !        ints_jl(i_ao,j_ao,i_cd) = dconjg(chol_ao_integrals_complex(j_ao,i_ao,i_cd,kk,-kpt_sparse_map(kQ)))
+        !      enddo
+        !    enddo
+        !  enddo
+        !endif
+
+        !call zgemm('N','T', ao_num_kpt_2, ao_num_kpt_2, chol_num(kQ), &
+        !       (1.d0,0.d0), ints_ik, ao_num_kpt_2, &
+        !       ints_jl, ao_num_kpt_2, &
+        !       (0.d0,0.d0), ints_ikjl, ao_num_kpt_2)
+
+        !n_integrals_1=0
+        !n_integrals_2=0
+        !do il=1,ao_num_per_kpt
+        !  l=il+(kl-1)*ao_num_per_kpt
+        !  do ij=1,ao_num_per_kpt
+        !    j=ij+(kj-1)*ao_num_per_kpt
+        !    if (j>l) exit
+        !    call idx2_tri_int(j,l,jl2)
+        !    do ik=1,ao_num_per_kpt
+        !      k=ik+(kk-1)*ao_num_per_kpt
+        !      if (k>l) exit
+        !      do ii=1,ao_num_per_kpt
+        !        i=ii+(ki-1)*ao_num_per_kpt
+        !        if ((j==l) .and. (i>k)) exit
+        !        call idx2_tri_int(i,k,ik2)
+        !        if (ik2 > jl2) exit
+        !        integral = ints_ikjl(ii,ik,ij,il)
+!       !         print*,i,k,j,l,real(integral),imag(integral)
+        !        if (cdabs(integral) < ao_integrals_threshold) then
+        !          cycle
+        !        endif
+        !        call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
+        !        tmp_re = dble(integral)
+        !        tmp_im = dimag(integral)
+        !        !if (use_map1) then
+        !        !  n_integrals_1 += 1
+        !        !  buffer_i_1(n_integrals_1)=idx_tmp
+        !        !  buffer_values_1(n_integrals_1)=tmp_re
+        !        !  if (sign.ne.0.d0) then 
+        !        !    n_integrals_1 += 1
+        !        !    buffer_i_1(n_integrals_1)=idx_tmp+1
+        !        !    buffer_values_1(n_integrals_1)=tmp_im*sign
+        !        !  endif 
+        !        !  if (n_integrals_1 >= size(buffer_i_1)-1) then
+        !        !    call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
+        !        !    n_integrals_1 = 0
+        !        !  endif
+        !        !else
+        !          !n_integrals_2 += 1
+        !          !buffer_i_2(n_integrals_2)=idx_tmp
+        !          !buffer_values_2(n_integrals_2)=tmp_re
+        !          !if (sign.ne.0.d0) then
+        !          !  n_integrals_2 += 1
+        !          !  buffer_i_2(n_integrals_2)=idx_tmp+1
+        !          !  buffer_values_2(n_integrals_2)=tmp_im*sign
+        !          !endif 
+        !          !if (n_integrals_2 >= size(buffer_i_2)-1) then
+        !          !  call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
+        !          !  n_integrals_2 = 0
+        !          !endif
+        !        endif
+
+        !      enddo !ii
+        !    enddo !ik
+        !  enddo !ij
+        !enddo !il
+
+        !if (n_integrals_1 > 0) then
+        !  call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
+        !endif
+        !if (n_integrals_2 > 0) then
+        !  call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
+        !endif
+      enddo !kk
+  !deallocate( &
+  !  ints_ik, &
+  !  ints_ikjl, &
+  !  buffer_i_1, &
+  !  buffer_i_2, &
+  !  buffer_values_1, &
+  !  buffer_values_2 &
+  !  )
+    enddo !kl
+    call wall_time(wall_2)
+    if (wall_2 - wall_0 > 1.d0) then
+      wall_0 = wall_2
+      !print*, 100.*float(kQ)/float(kpt_num), '% in ',               &
+      !    wall_2-wall_1,'s',map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB'
+    endif
+
+  enddo !kQ
+  !deallocate( ints_jl ) 
+
+  !call map_sort(ao_integrals_map)
+  !call map_unique(ao_integrals_map)
+  !call map_sort(ao_integrals_map_2)
+  !call map_unique(ao_integrals_map_2)
+  !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_1',ao_integrals_map)
+  !call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_complex_2',ao_integrals_map_2)
+  !call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+  
+  call wall_time(wall_2)
+  call cpu_time(cpu_2)
+
+  !integer*8                      :: get_ao_map_size, ao_map_size
+  !ao_map_size = get_ao_map_size()
+  
+  print*,'AO integrals provided:'
+  !print*,' Size of AO map           ', map_mb(ao_integrals_map),'+',map_mb(ao_integrals_map_2),'MB'
+  !print*,' Number of AO integrals: ',  ao_map_size
+  print*,' cpu  time :',cpu_2 - cpu_1, 's'
+  print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
+  
+end