fixed transformation (#116)

smaller three to four index transformation

* minor fix

* fixed integral transformation; added complex fcidump; fixed kpts bitmasks

src/bitmask/core_inact_act_virt.irp.f

@@ -448,7 +448,7 @@ BEGIN_PROVIDER [ integer, n_core_orb_kpts, (kpt_num)]
 
   do k=1,kpt_num
     n_core_orb_kpts(k) = 0
-    kshift = (1-k)*mo_num_per_kpt
+    kshift = (k-1)*mo_num_per_kpt
     do i = 1, mo_num_per_kpt
       if(mo_class(i+kshift) == 'Core')then
         n_core_orb_kpts(k) += 1
@@ -469,7 +469,7 @@ BEGIN_PROVIDER [ integer, n_inact_orb_kpts, (kpt_num)]
 
   do k=1,kpt_num
     n_inact_orb_kpts(k) = 0
-    kshift = (1-k)*mo_num_per_kpt
+    kshift = (k-1)*mo_num_per_kpt
     do i = 1, mo_num_per_kpt
       if(mo_class(i+kshift) == 'Inactive')then
         n_inact_orb_kpts(k) += 1
@@ -490,7 +490,7 @@ BEGIN_PROVIDER [ integer, n_act_orb_kpts, (kpt_num)]
 
   do k=1,kpt_num
     n_act_orb_kpts(k) = 0
-    kshift = (1-k)*mo_num_per_kpt
+    kshift = (k-1)*mo_num_per_kpt
     do i = 1, mo_num_per_kpt
       if(mo_class(i+kshift) == 'Active')then
         n_act_orb_kpts(k) += 1
@@ -511,7 +511,7 @@ BEGIN_PROVIDER [ integer, n_virt_orb_kpts, (kpt_num)]
 
   do k=1,kpt_num
     n_virt_orb_kpts(k) = 0
-    kshift = (1-k)*mo_num_per_kpt
+    kshift = (k-1)*mo_num_per_kpt
     do i = 1, mo_num_per_kpt
       if(mo_class(i+kshift) == 'Virtual')then
         n_virt_orb_kpts(k) += 1
@@ -532,7 +532,7 @@ BEGIN_PROVIDER [ integer, n_del_orb_kpts, (kpt_num)]
 
   do k=1,kpt_num
     n_del_orb_kpts(k) = 0
-    kshift = (1-k)*mo_num_per_kpt
+    kshift = (k-1)*mo_num_per_kpt
     do i = 1, mo_num_per_kpt
       if(mo_class(i+kshift) == 'Deleted')then
         n_del_orb_kpts(k) += 1

src/mo_two_e_ints/df_mo_ints.irp.f

@@ -48,7 +48,8 @@ subroutine mo_map_fill_from_df_dot
   logical :: use_map1
   integer(key_kind) :: idx_tmp
   double precision :: sign
-  complex*16, external :: zdotc
+  !complex*16, external :: zdotc
+  complex*16, external :: zdotu
 
   mo_num_kpt_2 = mo_num_per_kpt * mo_num_per_kpt
 
@@ -145,7 +146,8 @@ subroutine mo_map_fill_from_df_dot
                 if ((j==l) .and. (i>k)) exit
                 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 = zdotc(df_num,ints_jl(1,ij,il),1,ints_ik(1,ii,ik),1)
+                integral = zdotu(df_num,ints_jl(1,ij,il),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/tools/fcidump.irp.f

@@ -18,6 +18,97 @@ program fcidump
 !   electrons
 !
   END_DOC
+  if (is_complex) then
+    call fcidump_complex
+  else
+    call fcidump_real
+  endif
+end
+
+subroutine fcidump_complex
+  implicit none
+  character*(128) :: output
+  integer :: i_unit_output,getUnitAndOpen
+  output=trim(ezfio_filename)//'.FCIDUMP'
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  integer :: i,j,k,l
+  integer :: i1,j1,k1,l1
+  integer :: i2,j2,k2,l2,ik2,jl2
+  integer :: ki,kj,kk,kl
+  integer :: ii,ij,ik,il
+  integer*8 :: m
+  character*(2), allocatable :: A(:)
+
+  write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, &
+   ', MS2=', (elec_alpha_num-elec_beta_num), ','
+  allocate (A(n_act_orb))
+  A = '1,'
+  write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb)
+  write(i_unit_output,*) 'ISYM=0,'
+  write(i_unit_output,*) '/'
+  deallocate(A)
+
+  integer(key_kind), allocatable :: keys(:)
+  double precision, allocatable  :: values(:)
+  integer(cache_map_size_kind)   :: n_elements, n_elements_max
+  PROVIDE mo_two_e_integrals_in_map
+
+  complex*16 :: get_two_e_integral_complex, integral
+
+  do kl=1,kpt_num
+    do kj=1,kl
+      do kk=1,kl
+        ki=kconserv(kl,kk,kj)
+        if (ki>kl) cycle
+        do l1=1,n_act_orb_kpts(kl)
+          il=list_act_kpts(l1,kl)
+          l = (kl-1)*mo_num_per_kpt + il
+          do j1=1,n_act_orb_kpts(kj)
+            ij=list_act_kpts(j1,kj)
+            j = (kj-1)*mo_num_per_kpt + ij
+            if (j>l) exit
+            call idx2_tri_int(j,l,jl2)
+            do k1=1,n_act_orb_kpts(kk)
+              ik=list_act_kpts(k1,kk)
+              k = (kk-1)*mo_num_per_kpt + ik
+              if (k>l) exit
+              do i1=1,n_act_orb_kpts(ki)
+                ii=list_act_kpts(i1,ki)
+                i = (ki-1)*mo_num_per_kpt + ii
+                if ((j==l) .and. (i>k)) exit
+                call idx2_tri_int(i,k,ik2)
+                if (ik2 > jl2) exit
+                integral = get_two_e_integral_complex(i,j,k,l,mo_integrals_map,mo_integrals_map_2)
+                if (cdabs(integral) > mo_integrals_threshold) then
+                  write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral), dimag(integral),i,k,j,l
+                endif
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  do kj=1,kpt_num
+    do j1=1,n_act_orb_kpts(kj)
+      ij = list_act_kpts(j1,kj)
+      j = (kj-1)*mo_num_per_kpt + ij
+      do i1=j1,n_act_orb_kpts(kj)
+        ii = list_act_kpts(i1,kj)
+        i = (kj-1)*mo_num_per_kpt + ii
+        integral = mo_one_e_integrals_kpts(ii,ij,kj) + core_fock_operator_complex(i,j)
+        if (cdabs(integral) > mo_integrals_threshold) then
+          write(i_unit_output,'(2(E25.15,X),4(I6,X))') dble(integral),dimag(integral), i,j,0,0
+        endif
+      enddo
+    enddo
+  enddo
+  write(i_unit_output,*) core_energy, 0, 0, 0, 0
+end
+subroutine fcidump_real
+  implicit none
   character*(128) :: output
   integer :: i_unit_output,getUnitAndOpen
   output=trim(ezfio_filename)//'.FCIDUMP'

src/utils/constants.include.F

@@ -7,6 +7,7 @@ double precision, parameter :: sqpi =  dsqrt(dacos(-1.d0))
 double precision, parameter :: pi_5_2 =  34.9868366552d0
 double precision, parameter :: dfour_pi =  4.d0*dacos(-1.d0)
 double precision, parameter :: dtwo_pi =  2.d0*dacos(-1.d0)
+double precision, parameter :: inv_pi =  1.d0/dacos(-1.d0)
 double precision, parameter :: inv_sq_pi =  1.d0/dsqrt(dacos(-1.d0))
 double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
 double precision, parameter :: thresh = 1.d-15