cd fix

bin/qp_convert_h5_to_ezfio

@@ -288,8 +288,8 @@ 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 = (False, False)
-    #dump_fci, dump_cd = (True, False)
+    dump_fci, dump_cd = (True, False)
     #dump_fci, dump_cd = (False, True)
     #dump_fci2 = True
     dump_fci2 = False
@@ -594,14 +594,33 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
                         kr = QKTok2[Qi,kp]-1
                         for ks in range(kpt_num):
                             kq = QKTok2[Qi,ks]-1
+                            # 3
+                            #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)
+                            # 4
+                            #if Qneg:
+                            #    A = LQ0[kr].transpose((1,0,2)).conj()
+                            #    B = LQ0[kq].transpose((1,0,2))
+                            #    W = np.einsum('prn,sqn->pqrs',A,B)
+                            #else:
+                            #    A = LQ0[kp]
+                            #    B = LQ0[ks].conj()
+                            #    W = np.einsum('prn,sqn->pqrs',A,B)
+                            # 5
                             if Qneg:
                                 A = LQ0[kr].transpose((1,0,2)).conj()
-                                B = LQ0[kq]
+                                B = LQ0[kq].transpose((1,0,2))
-                                W = np.einsum('prn,sqn->pqrs',A,B)
+                                W = np.einsum('prn,qsn->pqrs',A,B)
                             else:
                                 A = LQ0[kp]
-                                B = LQ0[ks].transpose((1,0,2)).conj()
+                                B = LQ0[ks].conj()
-                                W = np.einsum('rpn,qsn->pqrs',A,B)
+                                W = np.einsum('prn,qsn->pqrs',A,B)
                             p0 = kp*ao_num_per_kpt
                             r0 = kr*ao_num_per_kpt
                             q0 = kq*ao_num_per_kpt
@@ -622,6 +641,11 @@ def convert_kpts_cd(filename,qph5path,qmcpack=True,is_ao=True):
                 E1=0
                 E2j=0
                 E2k=0
+                print("Jij Kij")
+                for i in range(ao_num_tot):
+                    for j in range(ao_num_tot):
+                        print(f'{i:5d} {j:5d} {i:5d} {j:5d}  {Wfull[i,j,i,j].real:25.15E} {Wfull[i,j,i,j].imag:25.15E}')
+                        print(f'{i:5d} {j:5d} {j:5d} {i:5d}  {Wfull[i,j,j,i].real:25.15E} {Wfull[i,j,j,i].imag:25.15E}')
                 for imo, iocc in enumerate(mo_occ):
                     if iocc:
                         E1 += 2*H1[imo,imo]

src/ao_two_e_ints/cd_ao_ints.irp.f

@@ -91,14 +91,14 @@ subroutine ao_map_fill_from_chol
       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))
+        !ints_jl = chol_ao_integrals_complex(:,:,:,kl,kpt_sparse_map(kQ))
-        !ints_jl = dconjg(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) = 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))
+              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

src/mo_two_e_ints/cd_mo_ints.irp.f

@@ -160,7 +160,8 @@ subroutine mo_map_fill_from_chol_dot
                 call idx2_tri_int(i,k,ik2)
                 if (ik2 > jl2) exit
                 !integral = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
-                integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
+                !integral = zdotu(chol_num(kQ),ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
+                integral = zdotu(chol_num(kQ),ints_jl(1,il,ij),1,ints_ik(1,ii,ik),1)
 !                print*,i,k,j,l,real(integral),imag(integral)
                 if (cdabs(integral) < mo_integrals_threshold) then
                   cycle

src/utils_complex/dump_ao_2e_cplx.irp.f

@@ -15,7 +15,9 @@ subroutine run
       do k=1,ao_num 
         do l=1,ao_num 
           tmp_cmplx = get_ao_two_e_integral_complex(i,j,k,l,ao_integrals_map,ao_integrals_map_2)
+          if (cdabs(tmp_cmplx) .gt. 1E-10) then
             print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx
+          endif
         enddo
       enddo
     enddo