updated converters and fixed ao df ints

src/ao_one_e_ints/ao_overlap.irp.f

@@ -260,7 +260,7 @@ BEGIN_PROVIDER [ complex*16, S_half_inv_complex, (AO_num,AO_num) ]
   integer                         :: info, i, j, k
   double precision, parameter     :: threshold_overlap_AO_eigenvalues = 1.d-6
 
-  LDA = size(AO_overlap,1)
+  LDA = size(AO_overlap_complex,1)
   LDC = size(S_half_inv_complex,1)
 
   allocate(         &

src/ao_two_e_ints/df_ao_ints.irp.f

@@ -116,8 +116,18 @@ subroutine ao_map_fill_from_df
   do kl=1, kpt_num
     do kj=1, kl
       call idx2_tri_int(kj,kl,kjkl2)
-      print*,'kj,kl,kjkl2',kj,kl,kjkl2
+      !print*,'kj,kl,kjkl2',kj,kl,kjkl2
-      ints_jl = df_ao_integrals_complex(:,:,:,kjkl2)
+      if (kj < kl) then
+        do i_ao=1,ao_num_per_kpt
+          do j_ao=1,ao_num_per_kpt
+            do i_df=1,df_num
+              ints_jl(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kjkl2))
+            enddo
+          enddo
+        enddo
+      else
+        ints_jl = df_ao_integrals_complex(:,:,:,kjkl2)
+      endif
 
   !$OMP PARALLEL PRIVATE(i,k,j,l,ki,kk,ii,ik,ij,il,kikk2,jl2,ik2, &
       !$OMP  ints_ik, ints_ikjl, i_ao, j_ao, i_df, &
@@ -143,15 +153,12 @@ subroutine ao_map_fill_from_df
         ki=kconserv(kl,kk,kj)
         if ((kl == kj) .and. (ki > kk)) cycle
         call idx2_tri_int(ki,kk,kikk2)
-        print*,'ki,kk,kikk2',ki,kk,kikk2
         if (kikk2 > kjkl2) cycle
-        if (ki >= kk) then
+        if (ki < kk) then
-        !if (ki < kk) then !this didn't fix the problem
           do i_ao=1,ao_num_per_kpt
             do j_ao=1,ao_num_per_kpt
               do i_df=1,df_num
                 ints_ik(i_ao,j_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2))
-                !ints_ik(j_ao,i_ao,i_df) = dconjg(df_ao_integrals_complex(j_ao,i_ao,i_df,kikk2))
               enddo
             enddo
           enddo

src/hartree_fock/hf_energy.irp.f

@@ -11,50 +11,50 @@ BEGIN_PROVIDER [double precision, extra_e_contrib_density]
 
 END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, HF_energy]
+ BEGIN_PROVIDER [ double precision, hf_energy]
-&BEGIN_PROVIDER [ double precision, HF_two_electron_energy]
+&BEGIN_PROVIDER [ double precision, hf_two_electron_energy]
-&BEGIN_PROVIDER [ double precision, HF_one_electron_energy]
+&BEGIN_PROVIDER [ double precision, hf_one_electron_energy]
  implicit none
  BEGIN_DOC
  ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
  END_DOC
  integer :: i,j
- HF_energy = nuclear_repulsion
+ hf_energy = nuclear_repulsion
- HF_two_electron_energy = 0.d0
+ hf_two_electron_energy = 0.d0
- HF_one_electron_energy = 0.d0
+ hf_one_electron_energy = 0.d0
   if (is_complex) then
     complex*16 :: hf_1e_tmp, hf_2e_tmp
     hf_1e_tmp = (0.d0,0.d0)
     hf_2e_tmp = (0.d0,0.d0)
     do j=1,ao_num
       do i=1,ao_num
-        hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_complex(i,j) * SCF_density_matrix_ao_alpha_complex(j,i) &
+        hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_complex(i,j) * scf_density_matrix_ao_alpha_complex(j,i) &
-                              +ao_two_e_integral_beta_complex(i,j)  * SCF_density_matrix_ao_beta_complex(j,i) )
+                              +ao_two_e_integral_beta_complex(i,j)  * scf_density_matrix_ao_beta_complex(j,i) )
-        hf_1e_tmp += ao_one_e_integrals_complex(i,j) * (SCF_density_matrix_ao_alpha_complex(j,i) &
+        hf_1e_tmp += ao_one_e_integrals_complex(i,j) * (scf_density_matrix_ao_alpha_complex(j,i) &
-                                                      + SCF_density_matrix_ao_beta_complex (j,i) )
+                                                      + scf_density_matrix_ao_beta_complex (j,i) )
       enddo
     enddo
     if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then
       print*,'HF_2e energy should be real:',irp_here
       stop -1
     else
-      HF_two_electron_energy = dble(hf_2e_tmp)
+      hf_two_electron_energy = dble(hf_2e_tmp)
     endif
     if (dabs(dimag(hf_1e_tmp)).gt.1.d-10) then
       print*,'HF_1e energy should be real:',irp_here
       stop -1
     else
-      HF_one_electron_energy = dble(hf_1e_tmp)
+      hf_one_electron_energy = dble(hf_1e_tmp)
     endif
   else
     do j=1,ao_num
       do i=1,ao_num
-        HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
+        hf_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * scf_density_matrix_ao_alpha(i,j) &
-                                           +ao_two_e_integral_beta(i,j)  * SCF_density_matrix_ao_beta(i,j) )
+                                           +ao_two_e_integral_beta(i,j)  * scf_density_matrix_ao_beta(i,j) )
-        HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+        hf_one_electron_energy += ao_one_e_integrals(i,j) * (scf_density_matrix_ao_alpha(i,j) + scf_density_matrix_ao_beta (i,j) )
       enddo
     enddo
   endif
- HF_energy += HF_two_electron_energy + HF_one_electron_energy
+ hf_energy += hf_two_electron_energy + hf_one_electron_energy
 END_PROVIDER
 

src/hartree_fock/print_e_scf.irp.f

@@ -0,0 +1,19 @@
+program print_e_scf
+  call run
+end
+
+subroutine run
+
+  use bitmasks
+  implicit none
+
+  call print_debug_scf_complex
+
+  print*,'hf 1e,2e,total energy'
+  print*,hf_one_electron_energy
+  print*,hf_two_electron_energy
+  print*,hf_energy
+  
+end
+
+

src/hartree_fock/scf.irp.f

@@ -98,6 +98,10 @@ subroutine run
     call roothaan_hall_scf
   endif
   call ezfio_set_hartree_fock_energy(SCF_energy)
+  print*,'hf 1e,2e,total energy'
+  print*,hf_one_electron_energy
+  print*,hf_two_electron_energy
+  print*,hf_energy
 
 end
 

src/scf_utils/print_debug_scf_complex.irp.f

@@ -15,26 +15,36 @@ subroutine print_debug_scf_complex
   do i=1,ao_num
     write(*,'(200(E24.15))') scf_density_matrix_ao_alpha_complex(i,:)
   enddo
-  write(*,'(A)') 'scf_density_matrix_ao_beta_complex'
+  write(*,'(A)') 'ao_one_e_integrals_complex'
   write(*,'(A)') '---------------'
   do i=1,ao_num
-    write(*,'(200(E24.15))') scf_density_matrix_ao_beta_complex(i,:)
+    write(*,'(200(E24.15))') ao_one_e_integrals_complex(i,:)
   enddo
   write(*,'(A)') 'ao_two_e_integral_alpha_complex'
   write(*,'(A)') '---------------'
   do i=1,ao_num
     write(*,'(200(E24.15))') ao_two_e_integral_alpha_complex(i,:)
   enddo
-  write(*,'(A)') 'ao_two_e_integral_beta_complex'
-  write(*,'(A)') '---------------'
-  do i=1,ao_num
-    write(*,'(200(E24.15))') ao_two_e_integral_beta_complex(i,:)
-  enddo
   write(*,'(A)') 'fock_matrix_ao_alpha_complex'
   write(*,'(A)') '---------------'
   do i=1,ao_num
     write(*,'(200(E24.15))') fock_matrix_ao_alpha_complex(i,:)
   enddo
+  write(*,'(A)') 'ao_overlap_complex'
+  write(*,'(A)') '---------------'
+  do i=1,ao_num
+    write(*,'(200(E24.15))') ao_overlap_complex(i,:)
+  enddo
+  write(*,'(A)') 'scf_density_matrix_ao_beta_complex'
+  write(*,'(A)') '---------------'
+  do i=1,ao_num
+    write(*,'(200(E24.15))') scf_density_matrix_ao_beta_complex(i,:)
+  enddo
+  write(*,'(A)') 'ao_two_e_integral_beta_complex'
+  write(*,'(A)') '---------------'
+  do i=1,ao_num
+    write(*,'(200(E24.15))') ao_two_e_integral_beta_complex(i,:)
+  enddo
   write(*,'(A)') 'fock_matrix_ao_beta_complex'
   write(*,'(A)') '---------------'
   do i=1,ao_num

src/scf_utils/scf_density_matrix_ao_complex.irp.f

@@ -1,4 +1,4 @@
-BEGIN_PROVIDER [complex*16, SCF_density_matrix_ao_alpha_complex, (ao_num,ao_num) ]
+BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_alpha_complex, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
    ! $C.C^t$ over $\alpha$ MOs
@@ -7,11 +7,11 @@ BEGIN_PROVIDER [complex*16, SCF_density_matrix_ao_alpha_complex, (ao_num,ao_num)
    call zgemm('N','C',ao_num,ao_num,elec_alpha_num,(1.d0,0.d0), &
         mo_coef_complex, size(mo_coef_complex,1), &
         mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), &
-        SCF_density_matrix_ao_alpha_complex, size(SCF_density_matrix_ao_alpha_complex,1))
+        scf_density_matrix_ao_alpha_complex, size(scf_density_matrix_ao_alpha_complex,1))
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_beta_complex,  (ao_num,ao_num) ]
+BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_beta_complex,  (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
    ! $C.C^t$ over $\beta$ MOs
@@ -20,21 +20,21 @@ BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_beta_complex,  (ao_num,ao_num
    call zgemm('N','C',ao_num,ao_num,elec_beta_num,(1.d0,0.d0), &
         mo_coef_complex, size(mo_coef_complex,1), &
         mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), &
-        SCF_density_matrix_ao_beta_complex, size(SCF_density_matrix_ao_beta_complex,1))
+        scf_density_matrix_ao_beta_complex, size(scf_density_matrix_ao_beta_complex,1))
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ complex*16, SCF_density_matrix_ao_complex, (ao_num,ao_num) ]
+BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_complex, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
    ! Sum of $\alpha$ and $\beta$ density matrices
    END_DOC
-   ASSERT (size(SCF_density_matrix_ao_complex,1) == size(SCF_density_matrix_ao_alpha_complex,1))
+   ASSERT (size(scf_density_matrix_ao_complex,1) == size(scf_density_matrix_ao_alpha_complex,1))
    if (elec_alpha_num== elec_beta_num) then
-     SCF_density_matrix_ao_complex = SCF_density_matrix_ao_alpha_complex + SCF_density_matrix_ao_alpha_complex
+     scf_density_matrix_ao_complex = scf_density_matrix_ao_alpha_complex + scf_density_matrix_ao_alpha_complex
    else
-     ASSERT (size(SCF_density_matrix_ao_complex,1) == size(SCF_density_matrix_ao_beta_complex ,1))
+     ASSERT (size(scf_density_matrix_ao_complex,1) == size(scf_density_matrix_ao_beta_complex ,1))
-     SCF_density_matrix_ao_complex = SCF_density_matrix_ao_alpha_complex + SCF_density_matrix_ao_beta_complex
+     scf_density_matrix_ao_complex = scf_density_matrix_ao_alpha_complex + scf_density_matrix_ao_beta_complex
    endif
 
 END_PROVIDER

src/utils_complex/Gen_Ezfio_from_integral_complex_3idx.sh

@@ -17,6 +17,9 @@ echo 'Create EZFIO'
 qp_edit -c $ezfio &> /dev/null
 #cp $ezfio/{ao,mo}_basis/ao_md5 
 
+qp_run import_kconserv $ezfio 
+#qp_run import_ao_2e_complex $ezfio 
+#qp_run dump_ao_2e_from_df $ezfio
 #Read the integral
 #echo 'Read Integral'
 

src/utils_complex/MolPyscfToQPkpts.py

@@ -508,6 +508,7 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
 
     natom = cell.natm
     nelec = cell.nelectron
+    neleca,nelecb = cell.nelec
     atom_xyz = mf.cell.atom_coords()
     if not(mf.cell.unit.startswith(('B','b','au','AU'))):
         atom_xyz /= nist.BOHR # always convert to au
@@ -537,8 +538,8 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
     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_alpha_num']=neleca*Nk 
-    qph5['electrons'].attrs['elec_beta_num']=nelec*Nk
+    qph5['electrons'].attrs['elec_beta_num']=nelecb*Nk
     qph5['mo_basis'].attrs['mo_num']=Nk*nmo
     qph5['ao_basis'].attrs['ao_num']=Nk*nao
     qph5['nuclei'].attrs['nucl_num']=Nk*natom
@@ -572,6 +573,18 @@ def pyscf2QP2(cell,mf, kpts, kmesh=None, cas_idx=None, int_threshold = 1E-8,
     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)
+    
+    with open('C.qp','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))
 
     # ___                                              
     #  |  ._ _|_  _   _  ._ _. |  _   |\/|  _  ._   _  

src/utils_complex/import_ao_2e_complex.irp.f

@@ -0,0 +1,89 @@
+program import_ao_2e_complex
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j,k,l
+  double precision :: integral
+  complex*16, allocatable       :: C(:,:)
+  double precision :: tmp_re, tmp_im
+
+  integer             :: n_integrals_1, n_integrals_2 
+  integer(key_kind), allocatable   :: buffer_i_1(:), buffer_i_2(:) 
+  real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:)
+  logical :: use_map1
+  integer(key_kind) :: idx_tmp
+  double precision :: sign
+
+
+!  call ezfio_set_ao_basis_ao_num(ao_num)
+
+  allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3))
+  allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3))
+  iunit = getunitandopen('W.qp','r')
+  n_integrals_1=0
+  n_integrals_2=0
+  buffer_values_1 = 0.d0
+  buffer_values_2 = 0.d0
+  do 
+    read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im
+    call ao_two_e_integral_complex_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
+!    print'(4(I4),(L3),(I6),(F7.1))',i,j,k,l,use_map1,idx_tmp,sign
+    if (use_map1) then
+      n_integrals_1 += 1
+      buffer_i_1(n_integrals_1)=idx_tmp
+      buffer_values_1(n_integrals_1)=tmp_re
+!      print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp,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
+!        print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp+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
+!      print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp,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
+!        print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp+1,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
+  13 continue
+  close(iunit)
+  
+  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
+
+  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')
+  
+end

src/utils_complex/import_kconserv.irp.f

@@ -22,7 +22,7 @@ subroutine run
   allocate(A(kpt_num,kpt_num,kpt_num))
   
   A = 0
-  iunit = getunitandopen('kconserv','r')
+  iunit = getunitandopen('K.qp','r')
   do 
     read (iunit,*,end=10) i,j,k,l
     A(i,j,k) = l