Merge pull request #168 from QuantumPackage/cleaning_dft

fixed stupid bug with iorder=-1 in integration.irp.f

src/dft_utils_in_r/mo_in_r.irp.f

@@ -161,3 +161,23 @@
   enddo
  enddo
  END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, mos_grad_in_r_array_transp_3, (3,n_points_final_grid,mo_num)]
+ implicit none
+ BEGIN_DOC
+! Transposed gradients 
+! 
+ END_DOC
+ integer :: i,j,m
+ double precision :: mos_array(mo_num), r(3)
+ double precision :: mos_grad_array(3,mo_num)
+ do m = 1, 3
+  do j = 1, mo_num
+   do i = 1, n_points_final_grid
+    mos_grad_in_r_array_transp_3(m,i,j) = mos_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+
+

src/mo_two_e_ints/four_idx_novvvv.irp.f

@@ -51,6 +51,13 @@ end
 
 
 subroutine four_idx_novvvv
+  print*,'********'
+  print*,'********'
+  print*,'********'
+  print*,'WARNING :: Using four_idx_novvvv, and we are not sure that this routine is not bugged ...'
+  print*,'********'
+  print*,'********'
+  print*,'********'
   use map_module
   implicit none
   BEGIN_DOC

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -50,7 +50,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
   call cpu_time(cpu_1)
 
   if(no_vvvv_integrals)then
-    call four_idx_novvvv
+!    call four_idx_novvvv
+    call four_idx_novvvv_old
   else
     call add_integrals_to_map(full_ijkl_bitmask_4)
   endif

src/mo_two_e_ints/no_vvvv.irp.f

@@ -0,0 +1,88 @@
+
+subroutine four_idx_novvvv_old
+  use map_module
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Retransform MO integrals for next CAS-SCF step
+  END_DOC
+  integer(bit_kind)              :: mask_ijkl(N_int,4)
+  integer(bit_kind)              :: mask_ijk(N_int,3)
+
+  print*,'Using partial transformation'
+  print*,'It will not transform all integrals with at least 3 indices within the virtuals'
+    integer                        :: i,j,k,l
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I I I !!!!!!!!!!!!!!!!!!!!
+    ! (core+inact+act) ^ 4
+    ! <ii|ii>
+    print*, ''
+    print*, '<ii|ii>'
+    do i = 1,N_int
+      mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,2) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,3) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,4) =  core_inact_act_bitmask_4(i,1)
+    enddo
+    call add_integrals_to_map(mask_ijkl)
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I V V !!!!!!!!!!!!!!!!!!!!
+    ! (core+inact+act) ^ 2  (virt) ^2
+    ! <iv|iv>  = J_iv
+    print*, ''
+    print*, '<iv|iv>'
+    do i = 1,N_int
+      mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,2) =  virt_bitmask(i,1)
+      mask_ijkl(i,3) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,4) =  virt_bitmask(i,1)
+    enddo
+    call add_integrals_to_map(mask_ijkl)
+
+    ! (core+inact+act) ^ 2  (virt) ^2
+    ! <ii|vv> = (iv|iv)
+    print*, ''
+    print*, '<ii|vv>'
+    do i = 1,N_int
+      mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,3) = virt_bitmask(i,1)
+      mask_ijkl(i,4) = virt_bitmask(i,1)
+    enddo
+    call add_integrals_to_map(mask_ijkl)
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!!
+!    if(.not.no_vvv_integrals)then
+      print*, ''
+      print*, '<rv|sv> and <rv|vs>'
+      do i = 1,N_int
+        mask_ijk(i,1) =  virt_bitmask(i,1)
+        mask_ijk(i,2) =  virt_bitmask(i,1)
+        mask_ijk(i,3) =  virt_bitmask(i,1)
+      enddo
+      call add_integrals_to_map_three_indices(mask_ijk)
+!    endif
+
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I I I V !!!!!!!!!!!!!!!!!!!!
+    ! (core+inact+act) ^ 3  (virt) ^1
+    ! <iv|ii>
+    print*, ''
+    print*, '<iv|ii>'
+    do i = 1,N_int
+      mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,2) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,3) =  core_inact_act_bitmask_4(i,1)
+      mask_ijkl(i,4) =  virt_bitmask(i,1)
+    enddo
+    call add_integrals_to_map(mask_ijkl)
+    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I V V V !!!!!!!!!!!!!!!!!!!!
+    ! (core+inact+act) ^ 1  (virt) ^3
+    ! <iv|vv>
+!    if(.not.no_ivvv_integrals)then
+      print*, ''
+      print*, '<iv|vv>'
+      do i = 1,N_int
+        mask_ijkl(i,1) =  core_inact_act_bitmask_4(i,1)
+        mask_ijkl(i,2) =  virt_bitmask(i,1)
+        mask_ijkl(i,3) =  virt_bitmask(i,1)
+        mask_ijkl(i,4) =  virt_bitmask(i,1)
+      enddo
+      call add_integrals_to_map_no_exit_34(mask_ijkl)
+end

src/utils/integration.irp.f

@@ -55,6 +55,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   !        fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
   !               * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
   !               * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
+  !
+  ! WARNING ::: IF fact_k is too smal then: 
+  ! returns a "s" function centered in zero
+  ! with an inifinite exponent and a zero polynom coef
   END_DOC
   implicit none
   include 'constants.include.F'
@@ -82,10 +86,13 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   !DIR$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
   if (fact_k < thresh) then
+    ! IF fact_k is too smal then: 
+    ! returns a "s" function centered in zero
+    ! with an inifinite exponent and a zero polynom coef
     P_center = 0.d0
-    p = 1.d-10
+    p = 1.d+15
     P_new = 0.d0
-    iorder = -1
+    iorder = 0
     fact_k = 0.d0
     return
   endif