finished complex mapping, starting comples hartree fock

src/ao_two_e_ints/map_integrals.irp.f

@@ -258,6 +258,7 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
  complex(integral_kind)         :: integral
   integer(key_kind)              :: p,q,r,s,ik,jl
   logical :: ilek, jlel, iklejl
+  complex*16 :: get_ao_two_e_integral_periodic_simple
 
 
  !$OMP PARALLEL DO PRIVATE (ilek,jlel,p,q,r,s, ik,jl,iklejl, &
@@ -267,36 +268,8 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
      do j=ao_integrals_cache_min,ao_integrals_cache_max
        do i=ao_integrals_cache_min,ao_integrals_cache_max
          !DIR$ FORCEINLINE
-         call two_e_integrals_index(i,j,k,l,idx1)
+         integral = get_ao_two_e_integral_periodic_simple(i,j,k,l,&
-         ilek = (i.le.k)
+                    ao_integrals_map,ao_integrals_map_2)
-         jlel = (j.le.l)
-         idx1 = 2*idx1 - 1
-         if (ilek.eqv.jlel) then !map1
-           !TODO: merge these calls using map_get_2
-           call map_get(ao_integrals_map,idx1,tmp_re)
-           call map_get(ao_integrals_map,idx1+1,tmp_im)
-           if (ilek) then
-             integral = dcmplx(tmp_re,tmp_im)
-           else
-             integral = dcmplx(tmp_re,-tmp_im)
-           endif
-         else !map2
-           !TODO: merge these calls using map_get_2
-           call map_get(ao_integrals_map_2,idx1,tmp_re)
-           call map_get(ao_integrals_map_2,idx1+1,tmp_im)
-           p = min(i,k)
-           r = max(i,k)
-           ik = p+shiftr(r*r-r,1)
-           q = min(j,l)
-           s = max(j,l)
-           jl = q+shiftr(s*s-s,1)
-           iklejl = (ik.le.jl)
-           if (ilek.eqv.iklejl) then
-             integral = dcmplx(tmp_re,tmp_im)
-           else
-             integral = dcmplx(tmp_re,-tmp_im)
-           endif
-        endif
          
          ii = l-ao_integrals_cache_min
          ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
@@ -324,35 +297,61 @@ subroutine ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign)
   integer(key_kind), intent(out) :: idx
   logical, intent(out)           :: use_map1
   double precision, intent(out)  :: sign
-  integer(key_kind)              :: p,q,r,s,ik,jl
+  integer(key_kind)              :: p,q,r,s,ik,jl,ij,kl
-  logical :: ilek, jlel, iklejl, ikeqjl
+  logical :: iltk, jltl, ikltjl, ieqk, jeql, ikeqjl, ijeqkl
   ! i.le.k, j.le.l, tri(i,k).le.tri(j,l)
   !DIR$ FORCEINLINE
   call two_e_integrals_index_periodic(i,j,k,l,idx,ik,jl)
-  ilek = (i.le.k)
+  p = min(i,j)
-  jlel = (j.le.l)
+  r = max(i,j)
+  ij = p+shiftr(r*r-r,1)
+  q = min(k,l)
+  s = max(k,l)
+  kl = q+shiftr(s*s-s,1)
+
+
   idx = 2*idx-1
-  ikeqjl = (ik.eq.jl)
+
-  if (ilek.eqv.jlel) then !map1
+  if (ij==kl) then !real, map1
+    sign=0.d0
     use_map1=.True.
+  else
+    iltk = (i.lt.k)
+    jltl = (j.lt.l)
+    ieqk = (i.eq.k)
+    jeql = (j.eq.l)
+    ikltjl = (ik.lt.jl)
+    ikeqjl = (ik.eq.jl)
     if (ikeqjl) then
-      sign=0.d0
+      if (iltk) then
-    else if (ilek) then
-      sign=1.d0
-    else
-      sign=-1.d0
-    endif
-  else !map2
-    use_map1=.False.
-    if (ikeqjl) then
-      sign=0.d0
-    else
-      iklejl = (ik.le.jl)
-      if (ilek.eqv.iklejl) then
         sign=1.d0
+        use_map1=.False.
       else
         sign=-1.d0
+        use_map1=.False.
       endif
+    else if (ieqk) then
+      if (jltl) then
+        sign=1.d0
+        use_map1=.True.
+      else
+        sign=-1.d0
+        use_map1=.True.
+      endif
+    else if (jeql) then
+      if (iltk) then
+        sign=1.d0
+        use_map1=.True.
+      else
+        sign=-1.d0
+        use_map1=.True.
+      endif
+    else if (iltk.eqv.ikltjl) then
+      sign=1.d0
+      use_map1=.False.
+    else
+      sign=-1.d0
+      use_map1=.False.
     endif
   endif
 end
@@ -364,48 +363,33 @@ complex*16 function get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) resu
   ! Gets one AO bi-electronic integral from the AO map
   END_DOC
   integer, intent(in)            :: i,j,k,l
-  integer(key_kind)              :: idx1,idx2
+  integer(key_kind)              :: idx1,idx2,idx
   real(integral_kind)            :: tmp_re, tmp_im
   integer(key_kind)              :: idx_re,idx_im
   type(map_type), intent(inout)  :: map,map2
   integer                        :: ii
   complex(integral_kind)         :: tmp
   integer(key_kind)              :: p,q,r,s,ik,jl
-  logical :: ilek, jlel, iklejl
+  logical :: ilek, jlel, iklejl,use_map1
+  double precision :: sign
   ! a.le.c, b.le.d, tri(a,c).le.tri(b,d)
   PROVIDE ao_two_e_integrals_in_map
-         !DIR$ FORCEINLINE
+  call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign)
-         call two_e_integrals_index(i,j,k,l,idx1)
+  if (use_map1) then
-         ilek = (i.le.k)
+    call map_get(map,idx,tmp_re)
-         jlel = (j.le.l)
+    if (sign/=0.d0) then
-         idx1 = idx1*2-1
+      call map_get(map,idx+1,tmp_im)
-         if (ilek.eqv.jlel) then !map1
+      tmp_im *= sign
-           !TODO: merge these calls using map_get_2
+    else
-           call map_get(map,idx1,tmp_re)
+      tmp_im=0.d0
-           call map_get(map,idx1+1,tmp_im)
+    endif
-           if (ilek) then
+  else
-             tmp = dcmplx(tmp_re,tmp_im)
+    call map_get(map2,idx,tmp_re)
-           else
+    call map_get(map2,idx+1,tmp_im)
-             tmp = dcmplx(tmp_re,-tmp_im)
+    tmp_im *= sign
-           endif
+  endif
-         else !map2
+  tmp = dcmplx(tmp_re,tmp_im)
-           !TODO: merge these calls using map_get_2
+  result = tmp
-           call map_get(map2,idx1,tmp_re)
-           call map_get(map2,idx1+1,tmp_im)
-           p = min(i,k)
-           r = max(i,k)
-           ik = p+shiftr(r*r-r,1)
-           q = min(j,l)
-           s = max(j,l)
-           jl = q+shiftr(s*s-s,1)
-           iklejl = (ik.le.jl)
-           if (ilek.eqv.iklejl) then
-             tmp = dcmplx(tmp_re,tmp_im)
-           else
-             tmp = dcmplx(tmp_re,-tmp_im)
-           endif
-        endif
-    result = tmp
 end
 
 
@@ -428,11 +412,12 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map,map2) result(resu
   ! a.le.c, b.le.d, tri(a,c).le.tri(b,d)
   PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min
   !DIR$ FORCEINLINE
-  if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
+!  if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
-    tmp = (0.d0,0.d0)
+!    tmp = (0.d0,0.d0)
-  else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then
+!  else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then
-    tmp = (0.d0,0.d0)
+!    tmp = (0.d0,0.d0)
-  else
+!  else
+  if (.True.) then
     ii = l-ao_integrals_cache_min
     ii = ior(ii, k-ao_integrals_cache_min)
     ii = ior(ii, j-ao_integrals_cache_min)

src/mo_basis/mos.irp.f

@@ -280,6 +280,35 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
   deallocate(T)
 end
 
+subroutine ao_to_mo_complex(A_ao,LDA_ao,A_mo,LDA_mo)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the AO basis to the MO basis
+  ! where A is complex in the AO basis
+  !
+  ! Ct.A_ao.C
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  complex*16, intent(in)   :: A_ao(LDA_ao,ao_num)
+  complex*16, intent(out)  :: A_mo(LDA_mo,mo_num)
+  complex*16, allocatable  :: T(:,:)
+  
+  allocate ( T(ao_num,mo_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call zgemm('N','N', ao_num, mo_num, ao_num,                    &
+      (1.d0,0.d0), A_ao,LDA_ao,                                      &
+      mo_coef_complex, size(mo_coef_complex,1),                                      &
+      (0.d0,0.d0), T, size(T,1))
+  
+  call zgemm('C','N', mo_num, mo_num, ao_num,                &
+      (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1),                          &
+      T, ao_num,                                                     &   
+      (0.d0,0.d0), A_mo, size(A_mo,1))
+  
+  deallocate(T)
+end
+
 
 subroutine mix_mo_jk(j,k)
   implicit none

src/mo_one_e_ints/ao_to_mo.irp.f

@@ -50,6 +50,61 @@ subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
   deallocate(T)
 end
 
+subroutine mo_to_ao_complex(A_mo,LDA_mo,A_ao,LDA_ao)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the MO basis to the AO basis
+  !
+  ! (S.C).A_mo.(S.C)t
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  complex*16, intent(in)         :: A_mo(LDA_mo,mo_num)
+  complex*16, intent(out)        :: A_ao(LDA_ao,ao_num)
+  complex*16, allocatable        :: T(:,:)
+  
+  allocate ( T(mo_num,ao_num) )
+  
+  call zgemm('N','C', mo_num, ao_num, mo_num,                &
+      (1.d0,0.d0), A_mo,size(A_mo,1),                                &
+      S_mo_coef_complex, size(S_mo_coef_complex,1),                                  &
+      (0.d0,0.d0), T, size(T,1))
+  
+  call zgemm('N','N', ao_num, ao_num, mo_num,                    &
+      (1.d0,0.d0), S_mo_coef_complex, size(S_mo_coef_complex,1),                     &
+      T, size(T,1),                                                  &
+      (0.d0,0.d0), A_ao, size(A_ao,1))
+  
+  deallocate(T)
+end
+
+subroutine mo_to_ao_no_overlap_complex(A_mo,LDA_mo,A_ao,LDA_ao)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the MO basis to the S^-1 AO basis
+  ! Useful for density matrix
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  complex*16, intent(in)         :: A_mo(LDA_mo,mo_num)
+  complex*16, intent(out)        :: A_ao(LDA_ao,ao_num)
+  complex*16, allocatable        :: T(:,:)
+  
+  allocate ( T(mo_num,ao_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call zgemm('N','C', mo_num, ao_num, mo_num,                &
+      (1.d0,0.d0), A_mo,size(A_mo,1),                                &
+      mo_coef_complex, size(mo_coef_complex,1),                                      &
+      (0.d0,0.d0), T, size(T,1))
+  
+  call zgemm('N','N', ao_num, ao_num, mo_num,                    &
+      (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1),                          &
+      T, size(T,1),                                                  &
+      (0.d0,0.d0), A_ao, size(A_ao,1))
+  
+  deallocate(T)
+end
+
+
 BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num, mo_num) ]
  implicit none
  BEGIN_DOC
@@ -63,4 +118,17 @@ BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num, mo_num) ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ complex*16, S_mo_coef_complex, (ao_num, mo_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Product S.C where S is the overlap matrix in the AO basis and C the mo_coef matrix.
+  END_DOC
+
+  call zgemm('N','N',ao_num, mo_num, ao_num, (1.d0,0.d0), &
+        ao_overlap_complex, size(ao_overlap_complex,1), &
+        mo_coef_complex, size(mo_coef_complex,1), &
+        (0.d0,0.d0), &
+        S_mo_coef_complex, size(S_mo_coef_complex,1))
+
+END_PROVIDER
 

src/scf_utils/fock_matrix.irp.f

@@ -101,18 +101,45 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_num,mo_num) ]
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
+   if (is_periodic) then
+     print*,'error',irp_here
+     stop -1
+   else
+     call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
                  Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1))
+   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Fock matrix on the MO basis
+   END_DOC
+   call ao_to_mo_complex(Fock_matrix_ao_alpha_complex,size(Fock_matrix_ao_alpha_complex,1), &
+                 Fock_matrix_mo_alpha_complex,size(Fock_matrix_mo_alpha_complex,1))
+END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_num,mo_num) ]
    implicit none
    BEGIN_DOC
    ! Fock matrix on the MO basis
    END_DOC
-   call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
+   if (is_periodic) then
+     print*,'error',irp_here
+     stop -1
+   else
+     call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
                  Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1))
+   endif
+END_PROVIDER
+
+BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Fock matrix on the MO basis
+   END_DOC
+   call ao_to_mo_complex(Fock_matrix_ao_beta_complex,size(Fock_matrix_ao_beta_complex,1), &
+                 Fock_matrix_mo_beta_complex,size(Fock_matrix_mo_beta_complex,1))
 END_PROVIDER
 
 
@@ -143,6 +170,33 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Fock matrix in AO basis set
+ END_DOC
+
+ if(frozen_orb_scf)then
+   call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1),              &
+       Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1))
+ else
+   if ( (elec_alpha_num == elec_beta_num).and.                       &
+         (level_shift == 0.) )                                       &
+         then
+     integer                        :: i,j
+     do j=1,ao_num
+       do i=1,ao_num
+         Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(i,j)
+       enddo
+     enddo
+   else
+     call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1),            &
+         Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1))
+   endif
+ endif
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, SCF_energy ]
  implicit none
  BEGIN_DOC

src/utils_periodic/export_integrals_ao_periodic.irp.f

@@ -153,6 +153,14 @@ provide ao_two_e_integrals_in_map
     call map_get(ao_integrals_map_2,idx_tmp+1,tmp6)
     print*,tmp3,tmp4
     print*,tmp5,tmp6
+    integer*8 :: ii
+          ii = l-ao_integrals_cache_min
+          ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
+          ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
+          ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
+    print*,'cache(pbc)=', ao_integrals_cache_periodic(ii)
+    print*,'cache(old)=', ao_integrals_cache(ii)
+    print*
 !    if (use_map1) then
 !      n_integrals_1 += 1
 !      buffer_i_1(n_integrals_1-1)=idx_tmp