TC SPRING CLEANING: BEGINNING

plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -123,7 +123,7 @@ subroutine give_integrals_3_body_bi_ort_spin( n, sigma_n, l, sigma_l, k, sigma_k
   endif
 
   return
-end subroutine give_integrals_3_body_bi_ort_spin
+end
 
 ! ---
 

plugins/local/non_h_ints_mu/jast_1e_utils.irp.f

@@ -132,6 +132,7 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
   double precision, allocatable :: A(:,:,:,:), b(:), A_tmp(:,:,:,:)
   double precision, allocatable :: Pa(:,:), Pb(:,:), Pt(:,:)
   double precision, allocatable :: u1e_tmp(:), tmp(:,:,:)
+  double precision, allocatable :: tmp1(:,:,:), tmp2(:,:,:)
   double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:)
 
 
@@ -176,18 +177,19 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
   ! --- --- ---
   ! get A
 
-  allocate(tmp(n_points_final_grid,ao_num,ao_num))
+  allocate(tmp1(n_points_final_grid,ao_num,ao_num), tmp2(n_points_final_grid,ao_num,ao_num))
   allocate(A(ao_num,ao_num,ao_num,ao_num))
 
   !$OMP PARALLEL               &
   !$OMP DEFAULT (NONE)         &
   !$OMP PRIVATE (i, j, ipoint) &
-  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp)
+  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp1, tmp2)
   !$OMP DO COLLAPSE(2)
   do j = 1, ao_num
     do i = 1, ao_num
       do ipoint = 1, n_points_final_grid
-        tmp(ipoint,i,j) = dsqrt(final_weight_at_r_vector(ipoint)) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp1(ipoint,i,j) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp2(ipoint,i,j) =                                    aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
       enddo
     enddo
   enddo
@@ -195,7 +197,7 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
   !$OMP END PARALLEL
 
   call dgemm( "T", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0  &
-            , tmp(1,1,1), n_points_final_grid, tmp(1,1,1), n_points_final_grid  &
+            , tmp1(1,1,1), n_points_final_grid, tmp2(1,1,1), n_points_final_grid &
             , 0.d0, A(1,1,1,1), ao_num*ao_num)
 
   allocate(A_tmp(ao_num,ao_num,ao_num,ao_num))
@@ -207,13 +209,13 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
   allocate(b(ao_num*ao_num))
 
   do ipoint = 1, n_points_final_grid
-    u1e_tmp(ipoint) = dsqrt(final_weight_at_r_vector(ipoint)) * u1e_tmp(ipoint)
+    u1e_tmp(ipoint) = u1e_tmp(ipoint)
   enddo
 
-  call dgemv("T", n_points_final_grid, ao_num*ao_num, 1.d0, tmp(1,1,1), n_points_final_grid, u1e_tmp(1), 1, 0.d0, b(1), 1)
+  call dgemv("T", n_points_final_grid, ao_num*ao_num, 1.d0, tmp1(1,1,1), n_points_final_grid, u1e_tmp(1), 1, 0.d0, b(1), 1)
 
   deallocate(u1e_tmp)
-  deallocate(tmp)
+  deallocate(tmp1, tmp2)
 
   ! --- --- ---
   ! solve Ax = b

plugins/local/non_h_ints_mu/numerical_integ.irp.f

@@ -179,7 +179,7 @@ double precision function num_v_ij_erf_rk_cst_mu_env(i, j, ipoint)
     dx    = r1(1) - r2(1)
     dy    = r1(2) - r2(2)
     dz    = r1(3) - r2(3)
-    r12   = dsqrt( dx * dx + dy * dy + dz * dz )
+    r12   = dsqrt(dx*dx + dy*dy + dz*dz)
     if(r12 .lt. 1d-10) cycle
 
     tmp1  = (derf(mu_erf * r12) - 1.d0) / r12
@@ -228,7 +228,7 @@ subroutine num_x_v_ij_erf_rk_cst_mu_env(i, j, ipoint, integ)
     dx    = r1(1) - r2(1)
     dy    = r1(2) - r2(2)
     dz    = r1(3) - r2(3)
-    r12   = dsqrt( dx * dx + dy * dy + dz * dz )
+    r12   = dsqrt(dx*dx + dy*dy + dz*dz)
     if(r12 .lt. 1d-10) cycle
 
     tmp1  = (derf(mu_erf * r12) - 1.d0) / r12
@@ -530,7 +530,7 @@ subroutine num_int2_u_grad1u_total_env2(i, j, ipoint, integ)
     dx    = r1(1) - r2(1)
     dy    = r1(2) - r2(2)
     dz    = r1(3) - r2(3)
-    r12   = dsqrt( dx * dx + dy * dy + dz * dz )
+    r12   = dsqrt(dx*dx + dy*dy + dz*dz)
     if(r12 .lt. 1d-10) cycle
 
     tmp0 = env_nucl(r2)

plugins/local/non_h_ints_mu/tc_integ_num.irp.f

@@ -66,9 +66,7 @@
     !$OMP PARALLEL                   &
     !$OMP DEFAULT (NONE)             &
     !$OMP PRIVATE (i_blocks, ipoint) &
-    !$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, &
+    !$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, final_grid_points, tmp_grad1_u12, tmp_grad1_u12_squared)
-    !$OMP         final_grid_points, tmp_grad1_u12,        &
-    !$OMP         tmp_grad1_u12_squared)
     !$OMP DO 
     do i_blocks = 1, n_blocks
       ipoint = ii - 1 + i_blocks ! r1
@@ -102,9 +100,7 @@
     !$OMP PARALLEL                 &
     !$OMP DEFAULT (NONE)           &
     !$OMP PRIVATE (i_rest, ipoint) &
-    !$OMP SHARED (n_rest, n_points_extra_final_grid, ii, &
+    !$OMP SHARED (n_rest, n_points_extra_final_grid, ii, final_grid_points, tmp_grad1_u12, tmp_grad1_u12_squared)
-    !$OMP         final_grid_points, tmp_grad1_u12,      &
-    !$OMP         tmp_grad1_u12_squared)
     !$OMP DO 
     do i_rest = 1, n_rest
       ipoint = ii - 1 + i_rest ! r1

plugins/local/non_h_ints_mu/test_non_h_ints.irp.f

@@ -1125,6 +1125,7 @@ subroutine test_fit_coef_A1()
   double precision              :: accu, norm, diff
   double precision, allocatable :: A1(:,:)
   double precision, allocatable :: A2(:,:,:,:), tmp(:,:,:)
+  double precision, allocatable :: tmp1(:,:,:), tmp2(:,:,:)
 
   ! ---
 
@@ -1165,16 +1166,17 @@ subroutine test_fit_coef_A1()
 
   call wall_time(t1)
 
-  allocate(tmp(ao_num,ao_num,n_points_final_grid))
+  allocate(tmp1(ao_num,ao_num,n_points_final_grid), tmp2(ao_num,ao_num,n_points_final_grid))
   !$OMP PARALLEL               &
   !$OMP DEFAULT (NONE)         &
   !$OMP PRIVATE (i, j, ipoint) &
-  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp)
+  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp1, tmp2)
   !$OMP DO COLLAPSE(2)
   do j = 1, ao_num
     do i = 1, ao_num
       do ipoint = 1, n_points_final_grid
-        tmp(i,j,ipoint) = dsqrt(final_weight_at_r_vector(ipoint)) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp1(i,j,ipoint) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp2(i,j,ipoint) =                                    aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
       enddo
     enddo
   enddo
@@ -1184,9 +1186,9 @@ subroutine test_fit_coef_A1()
   allocate(A2(ao_num,ao_num,ao_num,ao_num))
 
   call dgemm( "N", "T", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
-            , tmp(1,1,1), ao_num*ao_num, tmp(1,1,1), ao_num*ao_num              &
+            , tmp1(1,1,1), ao_num*ao_num, tmp2(1,1,1), ao_num*ao_num            &
             , 0.d0, A2(1,1,1,1), ao_num*ao_num)
-  deallocate(tmp)
+  deallocate(tmp1, tmp2)
 
   call wall_time(t2)
   print*, ' WALL TIME FOR A2 (min) =', (t2-t1)/60.d0
@@ -1238,6 +1240,7 @@ subroutine test_fit_coef_inv()
   double precision, allocatable :: A1(:,:), A1_inv(:,:), A1_tmp(:,:)
   double precision, allocatable :: A2(:,:,:,:), tmp(:,:,:), A2_inv(:,:,:,:)
   double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A2_tmp(:,:,:,:)
+  double precision, allocatable :: tmp1(:,:,:), tmp2(:,:,:)
 
   cutoff_svd = 5d-8
 
@@ -1286,16 +1289,17 @@ subroutine test_fit_coef_inv()
 
   call wall_time(t1)
 
-  allocate(tmp(n_points_final_grid,ao_num,ao_num))
+  allocate(tmp1(n_points_final_grid,ao_num,ao_num), tmp2(n_points_final_grid,ao_num,ao_num))
   !$OMP PARALLEL               &
   !$OMP DEFAULT (NONE)         &
   !$OMP PRIVATE (i, j, ipoint) &
-  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp)
+  !$OMP SHARED (n_points_final_grid, ao_num, final_weight_at_r_vector, aos_in_r_array_transp, tmp1, tmp2)
   !$OMP DO COLLAPSE(2)
   do j = 1, ao_num
     do i = 1, ao_num
       do ipoint = 1, n_points_final_grid
-        tmp(ipoint,i,j) = dsqrt(final_weight_at_r_vector(ipoint)) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp1(ipoint,i,j) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
+        tmp2(ipoint,i,j) =                                    aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j)
       enddo
     enddo
   enddo
@@ -1305,10 +1309,10 @@ subroutine test_fit_coef_inv()
   allocate(A2(ao_num,ao_num,ao_num,ao_num))
 
   call dgemm( "T", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0  &
-            , tmp(1,1,1), n_points_final_grid, tmp(1,1,1), n_points_final_grid  &
+            , tmp1(1,1,1), n_points_final_grid, tmp2(1,1,1), n_points_final_grid &
             , 0.d0, A2(1,1,1,1), ao_num*ao_num)
 
-  deallocate(tmp)
+  deallocate(tmp1, tmp2)
 
   call wall_time(t2)
   print*, ' WALL TIME FOR A2 (min) =', (t2-t1)/60.d0

plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f

@@ -273,60 +273,6 @@ end
 
 ! ---
 
-subroutine lapack_diag_non_sym_right(n, A, WR, WI, VR)
-
-  implicit none
-
-  integer,          intent(in)  :: n
-  double precision, intent(in)  :: A(n,n)
-  double precision, intent(out) :: WR(n), WI(n), VR(n,n)
-
-  integer                       :: i, lda, ldvl, ldvr, LWORK, INFO
-  double precision, allocatable :: Atmp(:,:), WORK(:), VL(:,:)
-
-  lda  = n
-  ldvl = 1
-  ldvr = n
-
-  allocate( Atmp(n,n), VL(1,1) )
-  Atmp(1:n,1:n) = A(1:n,1:n)
-
-  allocate(WORK(1))
-  LWORK = -1
-  call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
-  if(INFO.gt.0)then
-    print*,'dgeev failed !!',INFO
-    stop
-  endif
-
-  LWORK = max(int(WORK(1)), 1) ! this is the optimal size of WORK 
-  deallocate(WORK)
-
-  allocate(WORK(LWORK))
-
-  ! Actual diagonalization 
-  call dgeev('N', 'V', n, Atmp, lda, WR, WI, VL, ldvl, VR, ldvr, WORK, LWORK, INFO)
-  if(INFO.ne.0) then
-    print*,'dgeev failed !!', INFO
-    stop
-  endif
-
-  deallocate(Atmp, WORK, VL)
-
-! print *, ' JOBL = F'
-! print *, ' eigenvalues'
-! do i = 1, n
-!   write(*, '(1000(F16.10,X))') WR(i), WI(i)
-! enddo
-! print *, ' right eigenvect' 
-! do i = 1, n
-!   write(*, '(1000(F16.10,X))') VR(:,i)
-! enddo
-
-end
-
-! ---
-
 subroutine non_hrmt_real_diag(n, A, leigvec, reigvec, n_real_eigv, eigval)
 
   BEGIN_DOC
@@ -1780,70 +1726,6 @@ end
 
 ! ---
 
-subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_nd, S, stop_ifnot)
-
-  implicit none
-  
-  integer,          intent(in)  :: n, m
-  double precision, intent(in)  :: Vl(n,m), Vr(n,m), W(n,n)
-  double precision, intent(in)  :: thr_d, thr_nd
-  logical,          intent(in)  :: stop_ifnot
-  double precision, intent(out) :: accu_d, accu_nd, S(m,m)
-
-  integer                       :: i, j
-  double precision, allocatable :: SS(:,:), tmp(:,:)
-
-  print *, ' check weighted bi-orthogonality'
-
-  ! ---
-
-  allocate(tmp(m,n))
-  call dgemm( 'T', 'N', m, n, n, 1.d0        &
-            , Vl, size(Vl, 1), W, size(W, 1) &
-            , 0.d0, tmp, size(tmp, 1) )
-  call dgemm( 'N', 'N', m, m, n, 1.d0            &
-            , tmp, size(tmp, 1), Vr, size(Vr, 1) &
-            , 0.d0, S, size(S, 1) )
-  deallocate(tmp)
-
-  !print *, ' overlap matrix:'
-  !do i = 1, m
-  !  write(*,'(1000(F16.10,X))') S(i,:)
-  !enddo
-
-  accu_d  = 0.d0
-  accu_nd = 0.d0
-  do i = 1, m
-    do j = 1, m
-      if(i==j) then
-        accu_d = accu_d + dabs(S(i,i))
-      else
-        accu_nd = accu_nd + S(j,i) * S(j,i)
-      endif
-    enddo
-  enddo
-  accu_nd = dsqrt(accu_nd)
-
-  print *, ' accu_nd = ', accu_nd
-  print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
-
-  ! ---
-
-  if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
-    print *, ' non bi-orthogonal vectors !'
-    print *, ' accu_nd = ', accu_nd
-    print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
-    !print *, ' overlap matrix:'
-    !do i = 1, m
-    !  write(*,'(1000(F16.10,X))') S(i,:)
-    !enddo
-    stop
-  endif
-
-end
-
-! ---
-
 subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot)
 
   implicit none

plugins/local/non_hermit_dav/new_routines.irp.f

@@ -1,670 +0,0 @@
-subroutine non_hrmt_diag_split_degen_bi_orthog(n, A, leigvec, reigvec, n_real_eigv, eigval)
-
-  BEGIN_DOC
-  !
-  ! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors 
-  !
-  ! of a non hermitian matrix A(n,n)
-  !
-  ! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n" 
-  !
-  END_DOC
-
-  implicit none
-
-  integer,          intent(in)  :: n
-  double precision, intent(in)  :: A(n,n)
-  integer,          intent(out) :: n_real_eigv
-  double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
-  double precision, allocatable :: reigvec_tmp(:,:), leigvec_tmp(:,:)
-
-  integer                       :: i, j, n_degen,k , iteration
-  double precision              :: shift_current
-  double precision              :: r,thr,accu_d, accu_nd
-  integer,          allocatable :: iorder_origin(:),iorder(:)
-  double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:),S(:,:)
-  double precision, allocatable :: Aw(:,:),diag_elem(:),A_save(:,:)
-  double precision, allocatable :: im_part(:),re_part(:)
-  double precision :: accu,thr_cut, thr_norm=1d0
-
-
-  thr_cut = 1.d-15
-  print*,'Computing the left/right eigenvectors ...'
-  print*,'Using the degeneracy splitting algorithm'
- ! initialization
-  shift_current = 1.d-15
-  iteration = 0 
-  print*,'***** iteration = ',iteration
-
-
-  ! pre-processing the matrix :: sorting by diagonal elements
-  allocate(reigvec_tmp(n,n), leigvec_tmp(n,n))
-  allocate(diag_elem(n),iorder_origin(n),A_save(n,n))
-!  print*,'Aw'
-  do i = 1, n
-   iorder_origin(i) = i
-   diag_elem(i) = A(i,i)
-!   write(*,'(100(F16.10,X))')A(:,i)
-  enddo
-  call dsort(diag_elem, iorder_origin, n)
-  do i = 1, n
-   do j = 1, n
-    A_save(j,i) = A(iorder_origin(j),iorder_origin(i))
-   enddo
-  enddo
-
-  allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
-  allocate(im_part(n),iorder(n))
-  allocate( S(n,n) )
-
-
-  Aw = A_save
-  call cancel_small_elmts(aw,n,thr_cut)
-  call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-  do i = 1, n
-   im_part(i) = -dabs(WI(i))
-   iorder(i) = i
-  enddo
-  call dsort(im_part, iorder, n)
-  n_real_eigv = 0
-  do i = 1, n
-    if(dabs(WI(i)).lt.1.d-20)then
-      n_real_eigv += 1
-    else
-!      print*,'Found an imaginary component to eigenvalue'
-!      print*,'Re(i) + Im(i)',WR(i),WI(i)
-    endif
-  enddo
-  if(n_real_eigv.ne.n)then
-   shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-   print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-   print*,'Splitting the degeneracies by ',shift_current
-  else
-   print*,'All eigenvalues are real !'
-  endif
-
-
-  do while(n_real_eigv.ne.n)
-   iteration += 1
-   print*,'***** iteration = ',iteration
-   if(shift_current.gt.1.d-3)then
-    print*,'shift_current > 1.d-3 !!'
-    print*,'Your matrix intrinsically contains complex eigenvalues'
-    stop
-   endif
-   Aw = A_save
-   call cancel_small_elmts(Aw,n,thr_cut)
-   call split_matrix_degen(Aw,n,shift_current)
-   call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-   n_real_eigv = 0
-   do i = 1, n
-     if(dabs(WI(i)).lt.1.d-20)then
-       n_real_eigv+= 1
-     else
-!       print*,'Found an imaginary component to eigenvalue'
-!       print*,'Re(i) + Im(i)',WR(i),WI(i)
-     endif
-   enddo
-   if(n_real_eigv.ne.n)then
-    do i = 1, n
-     im_part(i) = -dabs(WI(i))
-     iorder(i) = i
-    enddo
-    call dsort(im_part, iorder, n)
-    shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-    print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-    print*,'Splitting the degeneracies by ',shift_current
-   else
-    print*,'All eigenvalues are real !'
-   endif
-  enddo
-  !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES 
-  do i = 1, n
-   eigval(i) = WR(i)
-   iorder(i) = i
-  enddo
-  call dsort(eigval,iorder,n)
-  do i = 1, n
-!   print*,'eigval(i) = ',eigval(i)
-   reigvec_tmp(:,i) = VR(:,iorder(i))
-   leigvec_tmp(:,i) = Vl(:,iorder(i))
-  enddo
-
-!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY
-  !                               check bi-orthogonality
-  call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else
-    print *, ' '
-    print *, ' bi-orthogonality: not imposed yet'
-    print *, ' '
-    print *, ' '
-    print *, ' orthog between degen eigenvect' 
-    print *, ' '
-    double precision, allocatable :: S_nh_inv_half(:,:)
-    allocate(S_nh_inv_half(n,n))
-    logical :: complex_root
-    deallocate(S_nh_inv_half)
-    call impose_orthog_degen_eigvec(n, eigval, reigvec_tmp)
-    call impose_orthog_degen_eigvec(n, eigval, leigvec_tmp)
-    call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-    if(accu_nd .lt. thresh_biorthog_nondiag) then
-      print *, ' bi-orthogonality: ok'
-    else 
-     print*,'New vectors not bi-orthonormals at ',accu_nd
-     call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S)
-     call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-     if(accu_nd .lt. thresh_biorthog_nondiag) then
-       print *, ' bi-orthogonality: ok'
-     else 
-      print*,'New vectors not bi-orthonormals at ',accu_nd
-      print*,'Must be a deep problem ...'
-      stop
-     endif
-    endif
-  endif
- 
-  !! EIGENVECTORS SORTED AND BI-ORTHONORMAL
-  do i = 1, n
-   do j = 1, n
-    VR(iorder_origin(j),i) = reigvec_tmp(j,i)
-    VL(iorder_origin(j),i) = leigvec_tmp(j,i)
-   enddo
-  enddo
-
-  !! RECOMPUTING THE EIGENVALUES 
-  eigval = 0.d0
-  do i = 1, n
-   iorder(i) = i
-   accu = 0.d0
-   do j = 1, n
-    accu += VL(j,i) * VR(j,i) 
-    do k = 1, n
-     eigval(i) +=  VL(j,i) * A(j,k) * VR(k,i) 
-    enddo
-   enddo
-   eigval(i) *= 1.d0/accu
-!   print*,'eigval(i) = ',eigval(i)
-  enddo
-  !! RESORT JUST TO BE SURE
-  call dsort(eigval, iorder, n)
-  do i = 1, n
-   do j = 1, n
-    reigvec(j,i) = VR(j,iorder(i))
-    leigvec(j,i) = VL(j,iorder(i))
-   enddo
-  enddo
-  print*,'Checking for final reigvec/leigvec'
-  shift_current = max(1.d-10,shift_current)
-  print*,'Thr for eigenvectors = ',shift_current
-  call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.)
-  call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else 
-   print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog'
-   print*,'Eigenvectors are not bi orthonormal ..'
-   print*,'accu_nd = ',accu_nd
-   stop
-  endif
-
-end 
-
-
-
-subroutine non_hrmt_diag_split_degen_s_inv_half(n, A, leigvec, reigvec, n_real_eigv, eigval)
-
-  BEGIN_DOC
-  !
-  ! routine which returns the sorted REAL EIGENVALUES ONLY and corresponding LEFT/RIGHT eigenvetors 
-  !
-  ! of a non hermitian matrix A(n,n)
-  !
-  ! n_real_eigv is the number of real eigenvalues, which might be smaller than the dimension "n" 
-  !
-  END_DOC
-
-  implicit none
-
-  integer,          intent(in)  :: n
-  double precision, intent(in)  :: A(n,n)
-  integer,          intent(out) :: n_real_eigv
-  double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
-  double precision, allocatable :: reigvec_tmp(:,:), leigvec_tmp(:,:)
-
-  integer                       :: i, j, n_degen,k , iteration
-  double precision              :: shift_current
-  double precision              :: r,thr,accu_d, accu_nd
-  integer,          allocatable :: iorder_origin(:),iorder(:)
-  double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:),S(:,:)
-  double precision, allocatable :: Aw(:,:),diag_elem(:),A_save(:,:)
-  double precision, allocatable :: im_part(:),re_part(:)
-  double precision :: accu,thr_cut, thr_norm=1.d0
-  double precision, allocatable :: S_nh_inv_half(:,:)
-  logical :: complex_root
-
-
-  thr_cut = 1.d-15
-  print*,'Computing the left/right eigenvectors ...'
-  print*,'Using the degeneracy splitting algorithm'
- ! initialization
-  shift_current = 1.d-15
-  iteration = 0 
-  print*,'***** iteration = ',iteration
-
-
-  ! pre-processing the matrix :: sorting by diagonal elements
-  allocate(reigvec_tmp(n,n), leigvec_tmp(n,n))
-  allocate(diag_elem(n),iorder_origin(n),A_save(n,n))
-!  print*,'Aw'
-  do i = 1, n
-   iorder_origin(i) = i
-   diag_elem(i) = A(i,i)
-!   write(*,'(100(F16.10,X))')A(:,i)
-  enddo
-  call dsort(diag_elem, iorder_origin, n)
-  do i = 1, n
-   do j = 1, n
-    A_save(j,i) = A(iorder_origin(j),iorder_origin(i))
-   enddo
-  enddo
-
-  allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
-  allocate(im_part(n),iorder(n))
-  allocate( S(n,n) )
-  allocate(S_nh_inv_half(n,n))
-
-
-  Aw = A_save
-  call cancel_small_elmts(aw,n,thr_cut)
-  call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-  do i = 1, n
-   im_part(i) = -dabs(WI(i))
-   iorder(i) = i
-  enddo
-  call dsort(im_part, iorder, n)
-  n_real_eigv = 0
-  do i = 1, n
-    if(dabs(WI(i)).lt.1.d-20)then
-      n_real_eigv += 1
-    else
-!      print*,'Found an imaginary component to eigenvalue'
-!      print*,'Re(i) + Im(i)',WR(i),WI(i)
-    endif
-  enddo
-  if(n_real_eigv.ne.n)then
-   shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-   print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-   print*,'Splitting the degeneracies by ',shift_current
-  else
-   print*,'All eigenvalues are real !'
-  endif
-
-
-  do while(n_real_eigv.ne.n)
-   iteration += 1
-   print*,'***** iteration = ',iteration
-   if(shift_current.gt.1.d-3)then
-    print*,'shift_current > 1.d-3 !!'
-    print*,'Your matrix intrinsically contains complex eigenvalues'
-    stop
-   endif
-   Aw = A_save
-!   thr_cut = shift_current
-   call cancel_small_elmts(Aw,n,thr_cut)
-   call split_matrix_degen(Aw,n,shift_current)
-   call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-   n_real_eigv = 0
-   do i = 1, n
-     if(dabs(WI(i)).lt.1.d-20)then
-       n_real_eigv+= 1
-     else
-!       print*,'Found an imaginary component to eigenvalue'
-!       print*,'Re(i) + Im(i)',WR(i),WI(i)
-     endif
-   enddo
-   if(n_real_eigv.ne.n)then
-    do i = 1, n
-     im_part(i) = -dabs(WI(i))
-     iorder(i) = i
-    enddo
-    call dsort(im_part, iorder, n)
-    shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-    print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-    print*,'Splitting the degeneracies by ',shift_current
-   else
-    print*,'All eigenvalues are real !'
-   endif
-  enddo
-  !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES 
-  do i = 1, n
-   eigval(i) = WR(i)
-   iorder(i) = i
-  enddo
-  call dsort(eigval,iorder,n)
-  do i = 1, n
-!   print*,'eigval(i) = ',eigval(i)
-   reigvec_tmp(:,i) = VR(:,iorder(i))
-   leigvec_tmp(:,i) = Vl(:,iorder(i))
-  enddo
-
-!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY
-  !                               check bi-orthogonality
-  call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else
-    print *, ' '
-    print *, ' bi-orthogonality: not imposed yet'
-    if(complex_root) then 
-     print *, ' '
-     print *, ' '
-     print *, ' orthog between degen eigenvect' 
-     print *, ' '
-     ! bi-orthonormalization using orthogonalization of left, right and then QR between left and right
-     call impose_orthog_degen_eigvec(n, eigval, reigvec_tmp) ! orthogonalization of reigvec
-     call impose_orthog_degen_eigvec(n, eigval, leigvec_tmp) ! orthogonalization of leigvec
-     call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S,  thresh_biorthog_diag, thresh_biorthog_nondiag, .false.) 
-
-     if(accu_nd .lt. thresh_biorthog_nondiag) then
-       print *, ' bi-orthogonality: ok'
-     else 
-      print*,'New vectors not bi-orthonormals at ', accu_nd
-      call get_inv_half_nonsymmat_diago(S, n, S_nh_inv_half, complex_root)
-      if(complex_root)then 
-       call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S) ! bi-orthonormalization using QR
-      else
-       print*,'S^{-1/2} exists !!'
-       call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization 
-      endif
-     endif
-    else ! the matrix S^{-1/2} exists
-     print*,'S^{-1/2} exists !!'
-     call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization 
-    endif
-    call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-    if(accu_nd .lt. thresh_biorthog_nondiag) then
-      print *, ' bi-orthogonality: ok'
-    else 
-     print*,'New vectors not bi-orthonormals at ',accu_nd
-     print*,'Must be a deep problem ...'
-     stop
-    endif
-  endif
- 
-  !! EIGENVECTORS SORTED AND BI-ORTHONORMAL
-  do i = 1, n
-   do j = 1, n
-    VR(iorder_origin(j),i) = reigvec_tmp(j,i)
-    VL(iorder_origin(j),i) = leigvec_tmp(j,i)
-   enddo
-  enddo
-
-  !! RECOMPUTING THE EIGENVALUES 
-  eigval = 0.d0
-  do i = 1, n
-   iorder(i) = i
-   accu = 0.d0
-   do j = 1, n
-    accu += VL(j,i) * VR(j,i) 
-    do k = 1, n
-     eigval(i) +=  VL(j,i) * A(j,k) * VR(k,i) 
-    enddo
-   enddo
-   eigval(i) *= 1.d0/accu
-!   print*,'eigval(i) = ',eigval(i)
-  enddo
-  !! RESORT JUST TO BE SURE
-  call dsort(eigval, iorder, n)
-  do i = 1, n
-   do j = 1, n
-    reigvec(j,i) = VR(j,iorder(i))
-    leigvec(j,i) = VL(j,iorder(i))
-   enddo
-  enddo
-  print*,'Checking for final reigvec/leigvec'
-  shift_current = max(1.d-10,shift_current)
-  print*,'Thr for eigenvectors = ',shift_current
-  call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.)
-  call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else 
-   print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog'
-   print*,'Eigenvectors are not bi orthonormal ..'
-   print*,'accu_nd = ',accu_nd
-   stop
-  endif
-
-end 
-
-
-subroutine non_hrmt_fock_mat(n, A, leigvec, reigvec, n_real_eigv, eigval)
-
-  BEGIN_DOC
-  !
-  ! routine returning the eigenvalues and left/right eigenvectors of the TC fock matrix 
-  !
-  END_DOC
-
-  implicit none
-
-  integer,          intent(in)  :: n
-  double precision, intent(in)  :: A(n,n)
-  integer,          intent(out) :: n_real_eigv
-  double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
-  double precision, allocatable :: reigvec_tmp(:,:), leigvec_tmp(:,:)
-
-  integer                       :: i, j, n_degen,k , iteration
-  double precision              :: shift_current
-  double precision              :: r,thr,accu_d, accu_nd
-  integer,          allocatable :: iorder_origin(:),iorder(:)
-  double precision, allocatable :: WR(:), WI(:), Vl(:,:), VR(:,:),S(:,:)
-  double precision, allocatable :: Aw(:,:),diag_elem(:),A_save(:,:)
-  double precision, allocatable :: im_part(:),re_part(:)
-  double precision :: accu,thr_cut
-  double precision, allocatable :: S_nh_inv_half(:,:)
-  logical :: complex_root
-  double precision :: thr_norm=1d0
-
-
-  thr_cut = 1.d-15
-  print*,'Computing the left/right eigenvectors ...'
-  print*,'Using the degeneracy splitting algorithm'
- ! initialization
-  shift_current = 1.d-15
-  iteration = 0 
-  print*,'***** iteration = ',iteration
-
-
-  ! pre-processing the matrix :: sorting by diagonal elements
-  allocate(reigvec_tmp(n,n), leigvec_tmp(n,n))
-  allocate(diag_elem(n),iorder_origin(n),A_save(n,n))
-!  print*,'Aw'
-  do i = 1, n
-   iorder_origin(i) = i
-   diag_elem(i) = A(i,i)
-!   write(*,'(100(F16.10,X))')A(:,i)
-  enddo
-  call dsort(diag_elem, iorder_origin, n)
-  do i = 1, n
-   do j = 1, n
-    A_save(j,i) = A(iorder_origin(j),iorder_origin(i))
-   enddo
-  enddo
-
-  allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
-  allocate(im_part(n),iorder(n))
-  allocate( S(n,n) )
-  allocate(S_nh_inv_half(n,n))
-
-
-  Aw = A_save
-  call cancel_small_elmts(aw,n,thr_cut)
-  call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-  do i = 1, n
-   im_part(i) = -dabs(WI(i))
-   iorder(i) = i
-  enddo
-  call dsort(im_part, iorder, n)
-  n_real_eigv = 0
-  do i = 1, n
-    if(dabs(WI(i)).lt.1.d-20)then
-      n_real_eigv += 1
-    else
-!      print*,'Found an imaginary component to eigenvalue'
-!      print*,'Re(i) + Im(i)',WR(i),WI(i)
-    endif
-  enddo
-  if(n_real_eigv.ne.n)then
-   shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-   print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-   print*,'Splitting the degeneracies by ',shift_current
-  else
-   print*,'All eigenvalues are real !'
-  endif
-
-
-  do while(n_real_eigv.ne.n)
-   iteration += 1
-   print*,'***** iteration = ',iteration
-   if(shift_current.gt.1.d-3)then
-    print*,'shift_current > 1.d-3 !!'
-    print*,'Your matrix intrinsically contains complex eigenvalues'
-    stop
-   endif
-   Aw = A_save
-!   thr_cut = shift_current
-   call cancel_small_elmts(Aw,n,thr_cut)
-   call split_matrix_degen(Aw,n,shift_current)
-   call lapack_diag_non_sym(n,Aw,WR,WI,VL,VR)
-   n_real_eigv = 0
-   do i = 1, n
-     if(dabs(WI(i)).lt.1.d-20)then
-       n_real_eigv+= 1
-     else
-!       print*,'Found an imaginary component to eigenvalue'
-!       print*,'Re(i) + Im(i)',WR(i),WI(i)
-     endif
-   enddo
-   if(n_real_eigv.ne.n)then
-    do i = 1, n
-     im_part(i) = -dabs(WI(i))
-     iorder(i) = i
-    enddo
-    call dsort(im_part, iorder, n)
-    shift_current = max(10.d0 * dabs(im_part(1)),shift_current*10.d0)
-    print*,'Largest imaginary part found in eigenvalues = ',im_part(1)
-    print*,'Splitting the degeneracies by ',shift_current
-   else
-    print*,'All eigenvalues are real !'
-   endif
-  enddo
-  !!!!!!!!!!!!!!!! SORTING THE EIGENVALUES 
-  do i = 1, n
-   eigval(i) = WR(i)
-   iorder(i) = i
-  enddo
-  call dsort(eigval,iorder,n)
-  do i = 1, n
-!   print*,'eigval(i) = ',eigval(i)
-   reigvec_tmp(:,i) = VR(:,iorder(i))
-   leigvec_tmp(:,i) = Vl(:,iorder(i))
-  enddo
-
-!!! ONCE ALL EIGENVALUES ARE REAL ::: CHECK BI-ORTHONORMALITY
-  !                               check bi-orthogonality
-  call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else
-    print *, ' '
-    print *, ' bi-orthogonality: not imposed yet'
-    print *, ' '
-    print *, ' '
-    print *, ' Using impose_unique_biorthog_degen_eigvec' 
-    print *, ' '
-    ! bi-orthonormalization using orthogonalization of left, right and then QR between left and right
-    call impose_unique_biorthog_degen_eigvec(n, eigval, mo_coef, leigvec_tmp, reigvec_tmp)
-    call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-    print*,'accu_nd = ',accu_nd
-    if(accu_nd .lt. thresh_biorthog_nondiag) then
-      print *, ' bi-orthogonality: ok'
-    else 
-     print*,'New vectors not bi-orthonormals at ',accu_nd
-     call get_inv_half_nonsymmat_diago(S, n, S_nh_inv_half,complex_root)
-     if(complex_root)then 
-      print*,'S^{-1/2} does not exits, using QR bi-orthogonalization'
-      call impose_biorthog_qr(n, n, leigvec_tmp, reigvec_tmp, S) ! bi-orthonormalization using QR
-     else
-      print*,'S^{-1/2} exists !!'
-      call bi_ortho_s_inv_half(n,leigvec_tmp,reigvec_tmp,S_nh_inv_half) ! use of S^{-1/2} bi-orthonormalization 
-     endif
-    endif
-    call check_biorthog(n, n, leigvec_tmp, reigvec_tmp, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-    if(accu_nd .lt. thresh_biorthog_nondiag) then
-      print *, ' bi-orthogonality: ok'
-    else 
-     print*,'New vectors not bi-orthonormals at ',accu_nd
-     print*,'Must be a deep problem ...'
-     stop
-    endif
-  endif
- 
-  !! EIGENVECTORS SORTED AND BI-ORTHONORMAL
-  do i = 1, n
-   do j = 1, n
-    VR(iorder_origin(j),i) = reigvec_tmp(j,i)
-    VL(iorder_origin(j),i) = leigvec_tmp(j,i)
-   enddo
-  enddo
-
-  !! RECOMPUTING THE EIGENVALUES 
-  eigval = 0.d0
-  do i = 1, n
-   iorder(i) = i
-   accu = 0.d0
-   do j = 1, n
-    accu += VL(j,i) * VR(j,i) 
-    do k = 1, n
-     eigval(i) +=  VL(j,i) * A(j,k) * VR(k,i) 
-    enddo
-   enddo
-   eigval(i) *= 1.d0/accu
-!   print*,'eigval(i) = ',eigval(i)
-  enddo
-  !! RESORT JUST TO BE SURE
-  call dsort(eigval, iorder, n)
-  do i = 1, n
-   do j = 1, n
-    reigvec(j,i) = VR(j,iorder(i))
-    leigvec(j,i) = VL(j,iorder(i))
-   enddo
-  enddo
-  print*,'Checking for final reigvec/leigvec'
-  shift_current = max(1.d-10,shift_current)
-  print*,'Thr for eigenvectors = ',shift_current
-  call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, shift_current, thr_norm, .false.)
-  call check_biorthog(n, n, leigvec, reigvec, accu_d, accu_nd, S, thresh_biorthog_diag, thresh_biorthog_nondiag, .false.)
-  print *, ' accu_nd bi-orthog = ', accu_nd
-  
-  if(accu_nd .lt. thresh_biorthog_nondiag) then
-    print *, ' bi-orthogonality: ok'
-  else 
-   print*,'Something went wrong in non_hrmt_diag_split_degen_bi_orthog'
-   print*,'Eigenvectors are not bi orthonormal ..'
-   print*,'accu_nd = ',accu_nd
-   stop
-  endif
-
-end 
-
-

plugins/local/ortho_three_e_ints/mu_j_ints_usual_mos.irp.f

@@ -183,11 +183,3 @@ BEGIN_PROVIDER [ double precision, x_W_ij_erf_rk, ( n_points_final_grid,3,mo_num
 
 END_PROVIDER 
 
-BEGIN_PROVIDER [ double precision, sqrt_weight_at_r, (n_points_final_grid)]
- implicit none
- integer :: ipoint
- do ipoint = 1, n_points_final_grid
-  sqrt_weight_at_r(ipoint) = dsqrt(final_weight_at_r_vector(ipoint))
- enddo
-END_PROVIDER 
-

plugins/local/tc_bi_ortho/EZFIO.cfg

@@ -9,3 +9,14 @@ interface: ezfio
 doc: Coefficients for the right wave function
 type: double precision
 size: (determinants.n_det,determinants.n_states)
+
+[tc_gs_energy]
+type: Threshold
+doc: TC GS Energy
+interface: ezfio
+
+[tc_gs_var]
+type: Threshold
+doc: TC GS VAR
+interface: ezfio
+

plugins/local/tc_bi_ortho/print_tc_energy.irp.f

@@ -6,18 +6,9 @@ program print_tc_energy
 
   implicit none
 
-  print *, 'Hello world'
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
   read_wf = .True.
   touch read_wf
 
-
   PROVIDE j2e_type
   PROVIDE j1e_type
   PROVIDE env_type
@@ -26,6 +17,27 @@ program print_tc_energy
   print *, ' j1e_type = ', j1e_type
   print *, ' env_type = ', env_type
 
+
+  my_grid_becke  = .True.
+  PROVIDE tc_grid1_a tc_grid1_r
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  call write_int(6, my_n_pt_r_grid, 'radial  external grid over')
+  call write_int(6, my_n_pt_a_grid, 'angular external grid over')
+
+  if(tc_integ_type .eq. "numeric") then
+    my_extra_grid_becke  = .True.
+    PROVIDE tc_grid2_a tc_grid2_r
+    my_n_pt_r_extra_grid = tc_grid2_r
+    my_n_pt_a_extra_grid = tc_grid2_a
+    touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
+
+    call write_int(6, my_n_pt_r_extra_grid, 'radial  internal grid over')
+    call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
+  endif
+
   call write_tc_energy()
 
 end

plugins/local/tc_bi_ortho/print_tc_var.irp.f

@@ -6,7 +6,8 @@ program print_tc_var
 
   implicit none
 
-  print *, 'Hello world'
+  print *, ' TC VAR is available only for HF REF WF'
+  print *, ' DO NOT FORGET TO RUN A CISD CALCULATION BEF'
 
   my_grid_becke = .True.
   PROVIDE tc_grid1_a tc_grid1_r
@@ -17,7 +18,7 @@ program print_tc_var
   read_wf = .True.
   touch read_wf
 
-  call write_tc_var()
+  call write_tc_gs_var_HF()
 
 end
 

plugins/local/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f

@@ -38,9 +38,9 @@ subroutine main()
         call ezfio_has_cisd_energy(exists)
         if(.not.exists) then
 
-          call ezfio_has_tc_scf_bitc_energy(exists)
+          call ezfio_has_tc_scf_tcscf_energy(exists)
           if(exists) then
-            call ezfio_get_tc_scf_bitc_energy(e_ref)
+            call ezfio_get_tc_scf_tcscf_energy(e_ref)
           endif
 
         else
@@ -59,7 +59,7 @@ subroutine main()
 
   close(iunit)
 
-end subroutine main
+end
 
 ! --
 
@@ -89,7 +89,7 @@ subroutine write_lr_spindeterminants()
   call ezfio_set_spindeterminants_psi_left_coef_matrix_values(buffer)
   deallocate(buffer)
 
-end subroutine write_lr_spindeterminants
+end
 
 ! ---
 

plugins/local/tc_bi_ortho/tc_utils.irp.f

@@ -6,8 +6,63 @@ subroutine write_tc_energy()
   double precision              :: hmono, htwoe, hthree, htot
   double precision              :: E_TC, O_TC
   double precision              :: E_1e, E_2e, E_3e
+  double precision, allocatable :: E_TC_tmp(:), E_1e_tmp(:), E_2e_tmp(:), E_3e_tmp(:)
 
-  do k = 1, n_states
+  ! GS 
+  ! ---
+
+  allocate(E_TC_tmp(N_det), E_1e_tmp(N_det), E_2e_tmp(N_det), E_3e_tmp(N_det))
+
+  !$OMP PARALLEL                                                                &
+  !$OMP DEFAULT (NONE)                                                          &
+  !$OMP PRIVATE(i, j, hmono, htwoe, hthree, htot)                               &
+  !$OMP SHARED(N_det, psi_det, N_int, psi_l_coef_bi_ortho, psi_r_coef_bi_ortho, &
+  !$OMP        E_TC_tmp, E_1e_tmp, E_2e_tmp, E_3e_tmp)
+  !$OMP DO
+  do i = 1, N_det
+    E_TC_tmp(i) = 0.d0
+    E_1e_tmp(i) = 0.d0
+    E_2e_tmp(i) = 0.d0
+    E_3e_tmp(i) = 0.d0
+    do j = 1, N_det
+      call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+      E_TC_tmp(i) = E_TC_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htot
+      E_1e_tmp(i) = E_1e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * hmono
+      E_2e_tmp(i) = E_2e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htwoe
+      E_3e_tmp(i) = E_3e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * hthree
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  E_1e = 0.d0
+  E_2e = 0.d0
+  E_3e = 0.d0
+  E_TC = 0.d0
+  O_TC = 0.d0
+  do i = 1, N_det
+    E_1e = E_1e + E_1e_tmp(i)
+    E_2e = E_2e + E_2e_tmp(i)
+    E_3e = E_3e + E_3e_tmp(i)
+    E_TC = E_TC + E_TC_tmp(i)
+    O_TC = O_TC + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(i,1)
+  enddo
+
+  print *, ' state :', 1
+  print *, " E_TC = ", E_TC / O_TC
+  print *, " E_1e = ", E_1e / O_TC
+  print *, " E_2e = ", E_2e / O_TC
+  print *, " E_3e = ", E_3e / O_TC
+  print *, " O_TC = ", O_TC
+
+  call ezfio_set_tc_bi_ortho_tc_gs_energy(E_TC/O_TC)
+
+  ! ---
+
+  ! ES 
+  ! ---
+
+  do k = 2, n_states
 
     E_TC = 0.d0
     E_1e = 0.d0
@@ -37,6 +92,8 @@ subroutine write_tc_energy()
 
   enddo
 
+  deallocate(E_TC_tmp, E_1e_tmp, E_2e_tmp, E_3e_tmp)
+
 end
 
 ! ---
@@ -66,3 +123,25 @@ end
 
 ! ---
 
+subroutine write_tc_gs_var_HF()
+
+  implicit none
+  integer          :: i, j, k
+  double precision :: hmono, htwoe, hthree, htot
+  double precision :: SIGMA_TC
+
+  SIGMA_TC = 0.d0
+  do j = 2, N_det
+    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
+    SIGMA_TC = SIGMA_TC + htot * htot
+  enddo
+
+  print *, " SIGMA_TC = ", SIGMA_TC
+
+  call ezfio_set_tc_bi_ortho_tc_gs_var(SIGMA_TC)
+
+end
+
+! ---
+
+

plugins/local/tc_scf/EZFIO.cfg

@@ -1,6 +1,6 @@
-[bitc_energy]
+[tcscf_energy]
 type: Threshold
-doc: Energy bi-tc HF
+doc: TC-SCF ENERGY 
 interface: ezfio
 
 [converged_tcscf]

plugins/local/tc_scf/combine_lr_tcscf.irp.f

@@ -1,75 +0,0 @@
-
-! ---
-
-program combine_lr_tcscf
-
-  BEGIN_DOC
-  ! TODO : Put the documentation of the program here
-  END_DOC
-
-  implicit none
-
-  my_grid_becke  = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  bi_ortho = .True.
-  touch bi_ortho
-
-  call comb_orbitals()
-
-end
-
-! ---
-
-subroutine comb_orbitals()
-
-  implicit none
-  integer                       :: i, m, n, nn, mm
-  double precision              :: accu_d, accu_nd
-  double precision, allocatable :: R(:,:), L(:,:), Rnew(:,:), tmp(:,:), S(:,:)
-
-  n  = ao_num
-  m  = mo_num
-  nn = elec_alpha_num
-  mm = m - nn
-
-  allocate(L(n,m), R(n,m), Rnew(n,m), S(m,m))
-  L = mo_l_coef
-  R = mo_r_coef
-
-  call check_weighted_biorthog(n, m, ao_overlap, L, R, accu_d, accu_nd, S, .true.)
-
-  allocate(tmp(n,nn))
-  do i = 1, nn 
-    tmp(1:n,i) = R(1:n,i)
-  enddo
-  call impose_weighted_orthog_svd(n, nn, ao_overlap, tmp)
-  do i = 1, nn
-    Rnew(1:n,i) = tmp(1:n,i)
-  enddo
-  deallocate(tmp)
-
-  allocate(tmp(n,mm))
-  do i = 1, mm
-    tmp(1:n,i) = L(1:n,i+nn)
-  enddo
-  call impose_weighted_orthog_svd(n, mm, ao_overlap, tmp)
-  do i = 1, mm
-    Rnew(1:n,i+nn) = tmp(1:n,i)
-  enddo
-  deallocate(tmp)
-
-  call check_weighted_biorthog(n, m, ao_overlap, Rnew, Rnew, accu_d, accu_nd, S, .true.)
-
-  mo_r_coef = Rnew
-  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
-
-  deallocate(L, R, Rnew, S)
-
-end subroutine comb_orbitals
-
-! ---
-

plugins/local/tc_scf/diago_vartcfock.irp.f

@@ -1,96 +0,0 @@
-
-! ---
-
-BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_mo, (mo_num, mo_num)]
-
-  implicit none
-
-  integer                       :: i, j
-  integer                       :: liwork, lwork, n, info
-  integer,          allocatable :: iwork(:)
-  double precision, allocatable :: work(:), F(:,:), F_save(:,:)
-  double precision, allocatable :: diag(:)
-
-  PROVIDE mo_r_coef
-  PROVIDE Fock_matrix_vartc_mo_tot
-
-  allocate( F(mo_num,mo_num), F_save(mo_num,mo_num)  )
-  allocate (diag(mo_num) )
-
-  do j = 1, mo_num
-    do i = 1, mo_num
-      F(i,j) = Fock_matrix_vartc_mo_tot(i,j)
-    enddo
-  enddo
-
-  ! Insert level shift here
-  do i = elec_beta_num+1, elec_alpha_num
-    F(i,i) += 0.5d0 * level_shift_tcscf
-  enddo
-  do i = elec_alpha_num+1, mo_num
-    F(i,i) += level_shift_tcscf
-  enddo
-
-  n = mo_num
-  lwork = 1+6*n + 2*n*n
-  liwork = 3 + 5*n
-
-  allocate(work(lwork))
-  allocate(iwork(liwork) )
-
-  lwork = -1
-  liwork = -1
-
-  F_save = F
-  call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info)
-
-  if (info /= 0) then
-    print *,  irp_here//' DSYEVD failed : ', info
-    stop 1
-  endif
-  lwork = int(work(1))
-  liwork = iwork(1)
-  deallocate(iwork)
-  deallocate(work)
-
-  allocate(work(lwork))
-  allocate(iwork(liwork) )
-  call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info)
-  deallocate(iwork)
-
-  if (info /= 0) then
-    F = F_save
-    call dsyev('V', 'L', mo_num, F, size(F, 1), diag, work, lwork, info)
-
-    if (info /= 0) then
-      print *,  irp_here//' DSYEV failed : ', info
-      stop 1
-    endif
-  endif
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-      fock_vartc_eigvec_mo(j,i) = F(j,i)
-    enddo
-  enddo
-
-  deallocate(work, F, F_save, diag)
-
-END_PROVIDER
-
-! ---
-
-BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_ao, (ao_num, mo_num)]
-
-  implicit none
-
-  PROVIDE mo_r_coef
-
-  call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0                                             &
-            , mo_r_coef, size(mo_r_coef, 1), fock_vartc_eigvec_mo, size(fock_vartc_eigvec_mo, 1) &
-            , 0.d0, fock_vartc_eigvec_ao, size(fock_vartc_eigvec_ao, 1))
-
-END_PROVIDER
-
-! ---
-

plugins/local/tc_scf/diis_tcscf.irp.f

@@ -91,29 +91,15 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
   double precision, allocatable :: tmp(:,:)
   double precision, allocatable :: F(:,:)
 
-  !print *, ' Providing FQS_SQF_ao ...'
+  PROVIDE Fock_matrix_tc_ao_tot
-  !call wall_time(t0)
 
   allocate(F(ao_num,ao_num))
-  if(var_tc) then
-
-    do i = 1, ao_num
-      do j = 1, ao_num
-        F(j,i) = Fock_matrix_vartc_ao_tot(j,i)
-      enddo
-    enddo
-
-  else
-
-    PROVIDE Fock_matrix_tc_ao_tot
   do i = 1, ao_num
     do j = 1, ao_num
       F(j,i) = Fock_matrix_tc_ao_tot(j,i)
     enddo
   enddo
 
-  endif
-
   allocate(tmp(ao_num,ao_num))
 
   ! F x Q
@@ -140,9 +126,6 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
   deallocate(tmp)
   deallocate(F)
 
-  !call wall_time(t1)
-  !print *, ' Wall time for FQS_SQF_ao =', t1-t0
-
 END_PROVIDER
 
 ! ---
@@ -152,61 +135,13 @@ BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
   implicit none
   double precision :: t0, t1
 
-  !print*, ' Providing FQS_SQF_mo ...'
-  !call wall_time(t0)
-
   PROVIDE mo_r_coef mo_l_coef
   PROVIDE FQS_SQF_ao
 
   call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) &
                         , FQS_SQF_mo, size(FQS_SQF_mo, 1) )
 
-  !call wall_time(t1)
-  !print*, ' Wall time for FQS_SQF_mo =', t1-t0 
-
 END_PROVIDER
 
 ! ---
 
-! BEGIN_PROVIDER [ double precision, eigenval_Fock_tc_ao, (ao_num) ]
-!&BEGIN_PROVIDER [ double precision, eigenvec_Fock_tc_ao, (ao_num,ao_num) ]
-!
-!  BEGIN_DOC
-!  !
-!  ! Eigenvalues and eigenvectors of the Fock matrix over the ao basis
-!  !
-!  ! F' = X.T x F x X   where X = ao_overlap^(-1/2)
-!  ! 
-!  ! F'   x Cr' = Cr' x E ==> F   Cr = Cr x E with Cr = X x Cr'
-!  ! F'.T x Cl' = Cl' x E ==> F.T Cl = Cl x E with Cl = X x Cl'
-!  !
-!  END_DOC
-!
-!  implicit none
-!  double precision, allocatable :: tmp1(:,:), tmp2(:,:)
-!
-!  ! ---
-!  ! Fock matrix in orthogonal basis: F' = X.T x F x X
-!
-!  allocate(tmp1(ao_num,ao_num))
-!  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0                                                 &
-!            , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), S_half_inv, size(S_half_inv, 1) &
-!            , 0.d0, tmp1, size(tmp1, 1) )
-!
-!  allocate(tmp2(ao_num,ao_num))
-!  call dgemm( 'T', 'N', ao_num, ao_num, ao_num, 1.d0               &
-!            , S_half_inv, size(S_half_inv, 1), tmp1, size(tmp1, 1) &
-!            , 0.d0, tmp2, size(tmp2, 1) )
-!
-!  ! ---
-!
-!  ! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
-!  ! TODO
-!
-!  ! Back-transform eigenvectors: C =X.C'
-!
-!END_PROVIDER
-
-! ---
-
-~                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    

plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f

@@ -1,299 +0,0 @@
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
-
-  implicit none
-  integer                       :: a, b, i, j, ipoint
-  double precision              :: ti, tf
-  double precision              :: loc_1, loc_2, loc_3
-  double precision, allocatable :: Okappa(:), Jkappa(:,:)
-  double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:)
-  double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:), tmp_22(:,:,:)
-  double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:)
-
-  PROVIDE mo_l_coef mo_r_coef
-
-  !print *, ' PROVIDING fock_3e_uhf_mo_cs ...'
-  !call wall_time(ti)
-
-  ! ---
-
-  allocate(Jkappa(n_points_final_grid,3), Okappa(n_points_final_grid))
-  Jkappa = 0.d0
-  Okappa = 0.d0
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, i, tmp_omp_d1, tmp_omp_d2)               &
-  !$OMP SHARED (n_points_final_grid, elec_beta_num,               &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         int2_grad1_u12_bimo_t, Okappa, Jkappa)
-
-  allocate(tmp_omp_d2(n_points_final_grid,3), tmp_omp_d1(n_points_final_grid))
-  tmp_omp_d2 = 0.d0
-  tmp_omp_d1 = 0.d0
-
-  !$OMP DO
-  do i = 1, elec_beta_num
-    do ipoint = 1, n_points_final_grid
-      tmp_omp_d2(ipoint,1) += int2_grad1_u12_bimo_t(ipoint,1,i,i)
-      tmp_omp_d2(ipoint,2) += int2_grad1_u12_bimo_t(ipoint,2,i,i)
-      tmp_omp_d2(ipoint,3) += int2_grad1_u12_bimo_t(ipoint,3,i,i)
-      tmp_omp_d1(ipoint)   += mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    Jkappa(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    Jkappa(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    Jkappa(ipoint,3) += tmp_omp_d2(ipoint,3) 
-    Okappa(ipoint)   += tmp_omp_d1(ipoint)   
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp_omp_d2, tmp_omp_d1)
-
-  !$OMP END PARALLEL
-
-  ! ---
-
-  allocate(tmp_1(n_points_final_grid,4))
-
-  do ipoint = 1, n_points_final_grid
-    loc_1 = 2.d0 * Okappa(ipoint) 
-    tmp_1(ipoint,1) = loc_1 * Jkappa(ipoint,1)
-    tmp_1(ipoint,2) = loc_1 * Jkappa(ipoint,2)
-    tmp_1(ipoint,3) = loc_1 * Jkappa(ipoint,3)
-    tmp_1(ipoint,4) = Okappa(ipoint)
-  enddo
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, i, j, loc_1, tmp_omp_d2)                 &
-  !$OMP SHARED (n_points_final_grid, elec_beta_num,               &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         int2_grad1_u12_bimo_t, tmp_1)
-
-  allocate(tmp_omp_d2(n_points_final_grid,3))
-  tmp_omp_d2 = 0.d0
-
-  !$OMP DO COLLAPSE(2)
-  do i = 1, elec_beta_num
-    do j = 1, elec_beta_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i) 
-
-        tmp_omp_d2(ipoint,1) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,i,j) 
-        tmp_omp_d2(ipoint,2) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,i,j) 
-        tmp_omp_d2(ipoint,3) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,i,j) 
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    tmp_1(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    tmp_1(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    tmp_1(ipoint,3) += tmp_omp_d2(ipoint,3) 
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp_omp_d2)
-  !$OMP END PARALLEL
-
-  ! ---
-
-  if(tc_save_mem) then
-
-    allocate(tmp_22(n_points_final_grid,4,mo_num))
-    do a = 1, mo_num
-      !$OMP PARALLEL                                                  &
-      !$OMP DEFAULT (NONE)                                            &
-      !$OMP PRIVATE (ipoint, b, i)                                    &
-      !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, a,    &
-      !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-      !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector,  &
-      !$OMP         tmp_22)
-      !$OMP DO
-      do b = 1, mo_num
-        do ipoint = 1, n_points_final_grid
-          tmp_22(ipoint,1,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
-          tmp_22(ipoint,2,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
-          tmp_22(ipoint,3,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
-        enddo
-        tmp_22(:,4,b) = 0.d0
-        do i = 1, elec_beta_num
-          do ipoint = 1, n_points_final_grid
-            tmp_22(ipoint,4,b) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
-                                                                     + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
-                                                                     + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
-          enddo
-        enddo
-      enddo
-      !$OMP END DO
-      !$OMP END PARALLEL
-      call dgemv( 'T', 4*n_points_final_grid, mo_num, -2.d0        &
-                , tmp_22(1,1,1), size(tmp_22, 1) * size(tmp_22, 2) &
-                , tmp_1(1,1), 1                                    &
-                , 0.d0, fock_3e_uhf_mo_cs(1,a), 1)
-    enddo
-    deallocate(tmp_22)
-
-  else
-
-    allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
-    !$OMP PARALLEL                                                  &
-    !$OMP DEFAULT (NONE)                                            &
-    !$OMP PRIVATE (ipoint, a, b, i)                                 &
-    !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num,       &
-    !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-    !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector,  &
-    !$OMP         tmp_2)
-    !$OMP DO COLLAPSE(2)
-    do a = 1, mo_num
-      do b = 1, mo_num
-        do ipoint = 1, n_points_final_grid
-          tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
-          tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
-          tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
-        enddo
-        tmp_2(:,4,b,a) = 0.d0
-        do i = 1, elec_beta_num
-          do ipoint = 1, n_points_final_grid
-            tmp_2(ipoint,4,b,a) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
-                                                                      + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
-                                                                      + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
-          enddo
-        enddo
-      enddo
-    enddo
-    !$OMP END DO
-    !$OMP END PARALLEL
-    call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 &
-              , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2)  &
-              , tmp_1(1,1), 1                                    &
-              , 0.d0, fock_3e_uhf_mo_cs(1,1), 1)
-    deallocate(tmp_2)
-
-  endif
-
-  deallocate(tmp_1)
-
-  ! ---
-
-  allocate(tmp_3(n_points_final_grid,5,mo_num), tmp_4(n_points_final_grid,5,mo_num))
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, b, loc_1, loc_2)                         &
-  !$OMP SHARED (n_points_final_grid, mo_num,                      &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         final_weight_at_r_vector, Jkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    tmp_3(:,:,b) = 0.d0
-    tmp_4(:,:,b) = 0.d0
-    do ipoint = 1, n_points_final_grid
-      tmp_3(ipoint,1,b) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,b)
-
-      tmp_4(ipoint,1,b) = -2.d0 * mos_r_in_r_array_transp(ipoint,b) * ( Jkappa(ipoint,1) * Jkappa(ipoint,1) &
-                                                                      + Jkappa(ipoint,2) * Jkappa(ipoint,2) &
-                                                                      + Jkappa(ipoint,3) * Jkappa(ipoint,3) )
-      tmp_4(ipoint,5,b) = mos_r_in_r_array_transp(ipoint,b)
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, b, i, loc_1, loc_2)                      &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num,       &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,  &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         Jkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    do i = 1, elec_beta_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i)
-        loc_2 = mos_r_in_r_array_transp(ipoint,i)
-
-        tmp_3(ipoint,2,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,b,i)
-        tmp_3(ipoint,3,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,b,i)
-        tmp_3(ipoint,4,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,b,i)
-        tmp_3(ipoint,5,b) += 2.d0 * loc_1 * ( Jkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,b,i) &
-                                            + Jkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,b,i) &
-                                            + Jkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,b,i) )
-                                                                                                       
-        tmp_4(ipoint,2,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b)
-        tmp_4(ipoint,3,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b)
-        tmp_4(ipoint,4,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b)
-        tmp_4(ipoint,1,b) += 2.d0 * loc_2 * ( Jkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,i,b) &
-                                            + Jkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,i,b) &
-                                            + Jkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,i,b) )
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, b, i, j, loc_1, loc_2, loc_3)            &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num,       &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,  &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    do i = 1, elec_beta_num
-      do j = 1, elec_beta_num
-        do ipoint = 1, n_points_final_grid
-
-          loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,j)
-          loc_2 = mos_r_in_r_array_transp(ipoint,b)
-          loc_3 = mos_r_in_r_array_transp(ipoint,i)
-
-          tmp_3(ipoint,5,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,j) )
-
-          tmp_4(ipoint,1,b) += ( loc_2 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i)   &
-                                         + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i)   &
-                                         + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) ) &
-                               - loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,b)   &
-                                         + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,b)   &
-                                         + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,b) ) )
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  ! ---
-
-  call dgemm( 'T', 'N', mo_num, mo_num, 5*n_points_final_grid, 1.d0 &
-            , tmp_3(1,1,1), 5*n_points_final_grid                   &
-            , tmp_4(1,1,1), 5*n_points_final_grid                   &
-            , 1.d0, fock_3e_uhf_mo_cs(1,1), mo_num)
-
-  deallocate(tmp_3, tmp_4)
-  deallocate(Jkappa, Okappa)
-
-  ! ---
-
-  !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_mo_cs =', (tf - ti) / 60.d0
-
-END_PROVIDER 
-
-! ---
-

plugins/local/tc_scf/fock_3e_bi_ortho_os.irp.f

@@ -1,536 +0,0 @@
-
-! ---
-
- BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a_os, (mo_num, mo_num)]
-&BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b_os, (mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! Open Shell part of the Fock matrix from three-electron terms 
-  !
-  ! WARNING :: non hermitian if bi-ortho MOS used 
-  !
-  END_DOC
-
-  implicit none
-  integer                       :: a, b, i, j, ipoint
-  double precision              :: loc_1, loc_2, loc_3, loc_4
-  double precision              :: ti, tf
-  double precision, allocatable :: Okappa(:), Jkappa(:,:), Obarkappa(:), Jbarkappa(:,:)
-  double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:)
-  double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:)
-  double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:)
-
-  PROVIDE mo_l_coef mo_r_coef
-
-  !print *, ' Providing fock_3e_uhf_mo_a_os and fock_3e_uhf_mo_b_os ...'
-  !call wall_time(ti)
-
-  ! ---
-
-  allocate(Jkappa(n_points_final_grid,3), Okappa(n_points_final_grid))
-  allocate(Jbarkappa(n_points_final_grid,3), Obarkappa(n_points_final_grid))
-  Jkappa    = 0.d0
-  Okappa    = 0.d0
-  Jbarkappa = 0.d0
-  Obarkappa = 0.d0
-
-  !$OMP PARALLEL                                                    &
-  !$OMP DEFAULT (NONE)                                              &
-  !$OMP PRIVATE (ipoint, i, tmp_omp_d1, tmp_omp_d2)                 &
-  !$OMP SHARED (n_points_final_grid, elec_beta_num, elec_alpha_num, &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
-  !$OMP         int2_grad1_u12_bimo_t, Okappa, Jkappa, Obarkappa, Jbarkappa)
-
-  allocate(tmp_omp_d2(n_points_final_grid,3), tmp_omp_d1(n_points_final_grid))
-
-  tmp_omp_d2 = 0.d0
-  tmp_omp_d1 = 0.d0
-  !$OMP DO
-  do i = 1, elec_beta_num
-    do ipoint = 1, n_points_final_grid
-      tmp_omp_d2(ipoint,1) += int2_grad1_u12_bimo_t(ipoint,1,i,i)
-      tmp_omp_d2(ipoint,2) += int2_grad1_u12_bimo_t(ipoint,2,i,i)
-      tmp_omp_d2(ipoint,3) += int2_grad1_u12_bimo_t(ipoint,3,i,i)
-      tmp_omp_d1(ipoint)   += mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    Jkappa(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    Jkappa(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    Jkappa(ipoint,3) += tmp_omp_d2(ipoint,3) 
-    Okappa(ipoint)   += tmp_omp_d1(ipoint)   
-  enddo
-  !$OMP END CRITICAL
-
-  tmp_omp_d2 = 0.d0
-  tmp_omp_d1 = 0.d0
-  !$OMP DO
-  do i = elec_beta_num+1, elec_alpha_num
-    do ipoint = 1, n_points_final_grid
-      tmp_omp_d2(ipoint,1) += int2_grad1_u12_bimo_t(ipoint,1,i,i)
-      tmp_omp_d2(ipoint,2) += int2_grad1_u12_bimo_t(ipoint,2,i,i)
-      tmp_omp_d2(ipoint,3) += int2_grad1_u12_bimo_t(ipoint,3,i,i)
-      tmp_omp_d1(ipoint)   += mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    Jbarkappa(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    Jbarkappa(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    Jbarkappa(ipoint,3) += tmp_omp_d2(ipoint,3) 
-    Obarkappa(ipoint)   += tmp_omp_d1(ipoint)   
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp_omp_d2, tmp_omp_d1)
-  !$OMP END PARALLEL
-
-  ! ---
-
-  allocate(tmp_1(n_points_final_grid,4))
-
-  do ipoint = 1, n_points_final_grid
-
-    loc_1 = -2.d0 * Okappa   (ipoint) 
-    loc_2 = -2.d0 * Obarkappa(ipoint) 
-    loc_3 =         Obarkappa(ipoint) 
-
-    tmp_1(ipoint,1) = (loc_1 - loc_3) * Jbarkappa(ipoint,1) + loc_2 * Jkappa(ipoint,1)
-    tmp_1(ipoint,2) = (loc_1 - loc_3) * Jbarkappa(ipoint,2) + loc_2 * Jkappa(ipoint,2)
-    tmp_1(ipoint,3) = (loc_1 - loc_3) * Jbarkappa(ipoint,3) + loc_2 * Jkappa(ipoint,3)
-
-    tmp_1(ipoint,4) = Obarkappa(ipoint)
-  enddo
-
-
-  !$OMP PARALLEL                                                    &
-  !$OMP DEFAULT (NONE)                                              &
-  !$OMP PRIVATE (ipoint, i, j, loc_1, loc_2, tmp_omp_d2)            &
-  !$OMP SHARED (n_points_final_grid, elec_beta_num, elec_alpha_num, &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
-  !$OMP         int2_grad1_u12_bimo_t, tmp_1)
-
-  allocate(tmp_omp_d2(n_points_final_grid,3))
-
-  tmp_omp_d2 = 0.d0
-  !$OMP DO COLLAPSE(2)
-  do i = 1, elec_beta_num
-    do j = elec_beta_num+1, elec_alpha_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i) 
-        loc_2 = mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j) 
-
-        tmp_omp_d2(ipoint,1) += loc_1 * int2_grad1_u12_bimo_t(ipoint,1,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,1,j,i) 
-        tmp_omp_d2(ipoint,2) += loc_1 * int2_grad1_u12_bimo_t(ipoint,2,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,2,j,i) 
-        tmp_omp_d2(ipoint,3) += loc_1 * int2_grad1_u12_bimo_t(ipoint,3,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,3,j,i) 
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    tmp_1(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    tmp_1(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    tmp_1(ipoint,3) += tmp_omp_d2(ipoint,3) 
-  enddo
-  !$OMP END CRITICAL
-
-  tmp_omp_d2 = 0.d0
-  !$OMP DO COLLAPSE(2)
-  do i = elec_beta_num+1, elec_alpha_num
-    do j = elec_beta_num+1, elec_alpha_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i) 
-
-        tmp_omp_d2(ipoint,1) += loc_1 * int2_grad1_u12_bimo_t(ipoint,1,i,j)
-        tmp_omp_d2(ipoint,2) += loc_1 * int2_grad1_u12_bimo_t(ipoint,2,i,j)
-        tmp_omp_d2(ipoint,3) += loc_1 * int2_grad1_u12_bimo_t(ipoint,3,i,j)
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-  !$OMP CRITICAL
-  do ipoint = 1, n_points_final_grid
-    tmp_1(ipoint,1) += tmp_omp_d2(ipoint,1) 
-    tmp_1(ipoint,2) += tmp_omp_d2(ipoint,2) 
-    tmp_1(ipoint,3) += tmp_omp_d2(ipoint,3) 
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp_omp_d2)
-  !$OMP END PARALLEL
-
-  ! ---
-
-  allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, a, b)                                    &
-  !$OMP SHARED (n_points_final_grid, mo_num,                      &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector,  &
-  !$OMP         tmp_2)
-  !$OMP DO COLLAPSE(2)
-  do a = 1, mo_num
-    do b = 1, mo_num
-      do ipoint = 1, n_points_final_grid
-        tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
-        tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
-        tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, a, b, i)                                           &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         tmp_2)
-  !$OMP DO COLLAPSE(2)
-  do a = 1, mo_num
-    do b = 1, mo_num
-
-      tmp_2(:,4,b,a) = 0.d0
-      do i = 1, elec_beta_num
-        do ipoint = 1, n_points_final_grid
-          tmp_2(ipoint,4,b,a) += final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
-                                                                    + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
-                                                                    + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  ! ---
-
-  call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, 1.d0 &
-            , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) &
-            , tmp_1(1,1), 1                                   &
-            , 0.d0, fock_3e_uhf_mo_b_os(1,1), 1)
-
-  deallocate(tmp_1, tmp_2)
-
-  ! ---
-
-  allocate(tmp_3(n_points_final_grid,2,mo_num), tmp_4(n_points_final_grid,2,mo_num))
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, b, loc_1, loc_2)                         &
-  !$OMP SHARED (n_points_final_grid, mo_num,                      &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         final_weight_at_r_vector, Jkappa, Jbarkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    tmp_3(:,:,b) = 0.d0
-    tmp_4(:,:,b) = 0.d0
-    do ipoint = 1, n_points_final_grid
-
-      tmp_3(ipoint,1,b) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,b)
-
-      loc_1 = -2.0d0 * mos_r_in_r_array_transp(ipoint,b)
-
-      tmp_4(ipoint,1,b) = loc_1 * ( Jbarkappa(ipoint,1) * (Jkappa(ipoint,1) + 0.25d0 * Jbarkappa(ipoint,1)) &
-                                  + Jbarkappa(ipoint,2) * (Jkappa(ipoint,2) + 0.25d0 * Jbarkappa(ipoint,2)) &
-                                  + Jbarkappa(ipoint,3) * (Jkappa(ipoint,3) + 0.25d0 * Jbarkappa(ipoint,3)) )
-
-      tmp_4(ipoint,2,b) = mos_r_in_r_array_transp(ipoint,b)
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, b, i, loc_1, loc_2, loc_3, loc_4)                  &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,           &
-  !$OMP         Jkappa, Jbarkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-
-    do i = 1, elec_beta_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i)
-        loc_2 = mos_r_in_r_array_transp(ipoint,i)
-
-        tmp_3(ipoint,2,b) += loc_1 * ( Jbarkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,b,i) &
-                                     + Jbarkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,b,i) &
-                                     + Jbarkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,b,i) )
-                                                                                                       
-        tmp_4(ipoint,1,b) += loc_2 * ( Jbarkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,i,b) &
-                                     + Jbarkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,i,b) &
-                                     + Jbarkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,i,b) )
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, b, i, j, loc_1, loc_2, loc_3)                      &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,           &
-  !$OMP         tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    do i = 1, elec_beta_num
-      do j = elec_beta_num+1, elec_alpha_num
-        do ipoint = 1, n_points_final_grid
-
-          loc_2 = mos_r_in_r_array_transp(ipoint,b)
-
-          tmp_4(ipoint,1,b) += loc_2 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) )
-        enddo
-      enddo
-    enddo
-
-    do i = elec_beta_num+1, elec_alpha_num
-      do j = elec_beta_num+1, elec_alpha_num
-        do ipoint = 1, n_points_final_grid
-
-          loc_2 = 0.5d0 * mos_r_in_r_array_transp(ipoint,b)
-
-          tmp_4(ipoint,1,b) += loc_2 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) )
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  ! ---
-
-  call dgemm( 'T', 'N', mo_num, mo_num, 2*n_points_final_grid, 1.d0 &
-            , tmp_3(1,1,1), 2*n_points_final_grid                   &
-            , tmp_4(1,1,1), 2*n_points_final_grid                   &
-            , 1.d0, fock_3e_uhf_mo_b_os(1,1), mo_num)
-
-  deallocate(tmp_3, tmp_4)
-
-
-
-
-  ! ---
-
-  fock_3e_uhf_mo_a_os = fock_3e_uhf_mo_b_os
-
-  allocate(tmp_1(n_points_final_grid,1))
-
-  do ipoint = 1, n_points_final_grid
-    tmp_1(ipoint,1) = Obarkappa(ipoint) + 2.d0 * Okappa(ipoint) 
-  enddo
-
-  allocate(tmp_2(n_points_final_grid,1,mo_num,mo_num))
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, a, b, i)                                           &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         tmp_2)
-  !$OMP DO COLLAPSE(2)
-  do a = 1, mo_num
-    do b = 1, mo_num
-
-      tmp_2(:,1,b,a) = 0.d0
-      do i = elec_beta_num+1, elec_alpha_num
-        do ipoint = 1, n_points_final_grid
-          tmp_2(ipoint,1,b,a) += final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
-                                                                    + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
-                                                                    + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  call dgemv( 'T', n_points_final_grid, mo_num*mo_num, 1.d0   &
-            , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) &
-            , tmp_1(1,1), 1                                   &
-            , 1.d0, fock_3e_uhf_mo_a_os(1,1), 1)
-
-  deallocate(tmp_1, tmp_2)
-
-  ! ---
-
-  allocate(tmp_3(n_points_final_grid,8,mo_num), tmp_4(n_points_final_grid,8,mo_num))
-
-  !$OMP PARALLEL                                                  &
-  !$OMP DEFAULT (NONE)                                            &
-  !$OMP PRIVATE (ipoint, b)                                       &
-  !$OMP SHARED (n_points_final_grid, mo_num,                      &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
-  !$OMP         final_weight_at_r_vector, Jkappa, Jbarkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    tmp_3(:,:,b) = 0.d0
-    tmp_4(:,:,b) = 0.d0
-    do ipoint = 1, n_points_final_grid
-
-      tmp_3(ipoint,1,b) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,b)
-
-      tmp_4(ipoint,8,b) = mos_r_in_r_array_transp(ipoint,b)
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, b, i, loc_1, loc_2, loc_3, loc_4)                  &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,           &
-  !$OMP         Jkappa, Jbarkappa, tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-    do i = 1, elec_beta_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i)
-        loc_2 = mos_r_in_r_array_transp(ipoint,i)
-
-        tmp_3(ipoint,2,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,b,i)
-        tmp_3(ipoint,3,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,b,i)
-        tmp_3(ipoint,4,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,b,i)
-
-        tmp_4(ipoint,5,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b)
-        tmp_4(ipoint,6,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b)
-        tmp_4(ipoint,7,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b)
-      enddo
-    enddo
-
-    do i = elec_beta_num+1, elec_alpha_num
-      do ipoint = 1, n_points_final_grid
-
-        loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i)
-        loc_3 = 2.d0 * loc_1
-        loc_2 = mos_r_in_r_array_transp(ipoint,i)
-        loc_4 = 2.d0 * loc_2
-
-        tmp_3(ipoint,5,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,b,i)
-        tmp_3(ipoint,6,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,b,i)
-        tmp_3(ipoint,7,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,b,i)
-
-        tmp_3(ipoint,8,b) += loc_3 * ( (Jkappa(ipoint,1) + 0.5d0 * Jbarkappa(ipoint,1)) * int2_grad1_u12_bimo_t(ipoint,1,b,i) &
-                                     + (Jkappa(ipoint,2) + 0.5d0 * Jbarkappa(ipoint,2)) * int2_grad1_u12_bimo_t(ipoint,2,b,i) &
-                                     + (Jkappa(ipoint,3) + 0.5d0 * Jbarkappa(ipoint,3)) * int2_grad1_u12_bimo_t(ipoint,3,b,i) )
-                                                                                                       
-        tmp_4(ipoint,1,b) += loc_4 * ( (Jkappa(ipoint,1) + 0.5d0 * Jbarkappa(ipoint,1)) * int2_grad1_u12_bimo_t(ipoint,1,i,b) &
-                                     + (Jkappa(ipoint,2) + 0.5d0 * Jbarkappa(ipoint,2)) * int2_grad1_u12_bimo_t(ipoint,2,i,b) &
-                                     + (Jkappa(ipoint,3) + 0.5d0 * Jbarkappa(ipoint,3)) * int2_grad1_u12_bimo_t(ipoint,3,i,b) )
-
-        tmp_4(ipoint,2,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b)
-        tmp_4(ipoint,3,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b)
-        tmp_4(ipoint,4,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b)
-
-        tmp_4(ipoint,5,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b)
-        tmp_4(ipoint,6,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b)
-        tmp_4(ipoint,7,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b)
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !$OMP PARALLEL                                                            &
-  !$OMP DEFAULT (NONE)                                                      &
-  !$OMP PRIVATE (ipoint, b, i, j, loc_1, loc_2, loc_3)                      &
-  !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, &
-  !$OMP         final_weight_at_r_vector, int2_grad1_u12_bimo_t,            &
-  !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp,           &
-  !$OMP         tmp_3, tmp_4)
-  !$OMP DO
-  do b = 1, mo_num
-
-    do i = 1, elec_beta_num
-      do j = elec_beta_num+1, elec_alpha_num
-        do ipoint = 1, n_points_final_grid
-
-          loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,j)
-          loc_2 = mos_r_in_r_array_transp(ipoint,b)
-          loc_3 = mos_r_in_r_array_transp(ipoint,i)
-
-          tmp_3(ipoint,8,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,j) )
-
-          tmp_4(ipoint,1,b) -= loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,b) )
-
-          loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i)
-          loc_3 = mos_r_in_r_array_transp(ipoint,j)
-
-          tmp_3(ipoint,8,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,b,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,b,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) )
-
-          tmp_4(ipoint,1,b) -= loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,j,i) * int2_grad1_u12_bimo_t(ipoint,1,i,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,j,i) * int2_grad1_u12_bimo_t(ipoint,2,i,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,j,i) * int2_grad1_u12_bimo_t(ipoint,3,i,b) )
-        enddo
-      enddo
-    enddo
-
-    do i = elec_beta_num+1, elec_alpha_num
-      do j = elec_beta_num+1, elec_alpha_num
-        do ipoint = 1, n_points_final_grid
-
-          loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,j)
-          loc_2 = 0.5d0 * mos_r_in_r_array_transp(ipoint,b)
-          loc_3 = mos_r_in_r_array_transp(ipoint,i)
-
-          tmp_3(ipoint,8,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,j) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,j) )
-
-          tmp_4(ipoint,1,b) -= loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,b) &
-                                       + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,b) )
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  ! ---
-
-  call dgemm( 'T', 'N', mo_num, mo_num, 8*n_points_final_grid, 1.d0 &
-            , tmp_3(1,1,1), 8*n_points_final_grid                   &
-            , tmp_4(1,1,1), 8*n_points_final_grid                   &
-            , 1.d0, fock_3e_uhf_mo_a_os(1,1), mo_num)
-
-  deallocate(tmp_3, tmp_4)
-  deallocate(Jkappa, Okappa)
-
-  !call wall_time(tf)
-  !print *, ' Wall time for fock_3e_uhf_mo_a_os and fock_3e_uhf_mo_b_os =', tf - ti
-
-END_PROVIDER 
-
-! ---
-

plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f

@@ -1,77 +0,0 @@
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! Fock matrix alpha from three-electron terms 
-  !
-  ! WARNING :: non hermitian if bi-ortho MOS used 
-  !
-  END_DOC
-
-  implicit none
-  double precision :: ti, tf
-
-  PROVIDE mo_l_coef mo_r_coef
-
-  !print *, ' Providing fock_3e_uhf_mo_a ...'
-  !call wall_time(ti)
-
-  ! CLOSED-SHELL PART
-  PROVIDE fock_3e_uhf_mo_cs
-  fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs
-
-  if(elec_alpha_num .ne. elec_beta_num) then
-
-    ! OPEN-SHELL PART
-    PROVIDE fock_3e_uhf_mo_a_os
-
-    fock_3e_uhf_mo_a += fock_3e_uhf_mo_a_os
-  endif
-
-  !call wall_time(tf)
-  !print *, ' Wall time for fock_3e_uhf_mo_a (min) =', (tf - ti)/60.d0
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! Fock matrix beta from three-electron terms 
-  !
-  ! WARNING :: non hermitian if bi-ortho MOS used 
-  !
-  END_DOC
-
-  implicit none
-  double precision :: ti, tf
-
-  PROVIDE mo_l_coef mo_r_coef
-
-  !print *, ' Providing and fock_3e_uhf_mo_b ...'
-  !call wall_time(ti)
-
-  ! CLOSED-SHELL PART
-  PROVIDE fock_3e_uhf_mo_cs
-  fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs
-
-  if(elec_alpha_num .ne. elec_beta_num) then
-
-    ! OPEN-SHELL PART
-    PROVIDE fock_3e_uhf_mo_b_os
-
-    fock_3e_uhf_mo_b += fock_3e_uhf_mo_b_os
-  endif
-
-  !call wall_time(tf)
-  !print *, ' Wall time for fock_3e_uhf_mo_b =', tf - ti
-
-END_PROVIDER 
-
-! ---
-

plugins/local/tc_scf/fock_3e_bi_ortho_uhf_old.irp.f

@@ -1,490 +0,0 @@
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs_old, (mo_num, mo_num)]
-
-  implicit none
-  integer                       :: a, b, i, j
-  double precision              :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
-  double precision              :: ti, tf
-  double precision, allocatable :: tmp(:,:)
-
-  PROVIDE mo_l_coef mo_r_coef
-  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
-
-  !print *, ' PROVIDING fock_3e_uhf_mo_cs_old ...'
-  !call wall_time(ti)
-
-  fock_3e_uhf_mo_cs_old = 0.d0
-
-  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
-  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
-  !$OMP SHARED  (mo_num, elec_beta_num, fock_3e_uhf_mo_cs_old)
-
-  allocate(tmp(mo_num,mo_num))
-  tmp = 0.d0
-
-  !$OMP DO
-  do a = 1, mo_num
-    do b = 1, mo_num
-   
-      do j = 1, elec_beta_num
-        do i = 1, elec_beta_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij &
-                              +        I_bij_ija &
-                              +        I_bij_jai &
-                              - 2.d0 * I_bij_aji &
-                              - 2.d0 * I_bij_iaj &
-                              - 2.d0 * I_bij_jia )
-
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do a = 1, mo_num
-    do b = 1, mo_num
-      fock_3e_uhf_mo_cs_old(b,a) += tmp(b,a)
-    enddo
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp)
-  !$OMP END PARALLEL
-
-  !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_mo_cs_old =', tf - ti
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a_old, (mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! ALPHA part of the Fock matrix from three-electron terms 
-  !
-  ! WARNING :: non hermitian if bi-ortho MOS used 
-  !
-  END_DOC
-
-  implicit none
-  integer                       :: a, b, i, j, o
-  double precision              :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
-  double precision              :: ti, tf
-  double precision, allocatable :: tmp(:,:)
-
-  PROVIDE mo_l_coef mo_r_coef
-  PROVIDE fock_3e_uhf_mo_cs
-
-  !print *, ' Providing fock_3e_uhf_mo_a_old ...'
-  !call wall_time(ti)
-
-  o = elec_beta_num + 1
-  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
-
-  PROVIDE fock_3e_uhf_mo_cs_old
-  fock_3e_uhf_mo_a_old = fock_3e_uhf_mo_cs_old
-
-  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
-  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
-  !$OMP SHARED  (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_a_old)
-
-  allocate(tmp(mo_num,mo_num))
-  tmp = 0.d0
-
-  !$OMP DO
-  do a = 1, mo_num
-    do b = 1, mo_num
-
-      ! ---
-
-      do j = o, elec_alpha_num
-        do i = 1, elec_beta_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
-                              +        I_bij_ija &
-                              +        I_bij_jai &
-                              -        I_bij_aji &
-                              -        I_bij_iaj &
-                              - 2.d0 * I_bij_jia )
-
-        enddo
-      enddo
-
-      ! ---
-
-      do j = 1, elec_beta_num
-        do i = o, elec_alpha_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
-                              +        I_bij_ija &
-                              +        I_bij_jai &
-                              -        I_bij_aji &
-                              - 2.d0 * I_bij_iaj &
-                              -        I_bij_jia )
-
-        enddo
-      enddo
-
-      ! ---
-
-      do j = o, elec_alpha_num
-        do i = o, elec_alpha_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( I_bij_aij &
-                              + I_bij_ija &
-                              + I_bij_jai &
-                              - I_bij_aji &
-                              - I_bij_iaj &
-                              - I_bij_jia )
-
-        enddo
-      enddo
-
-      ! ---
-
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do a = 1, mo_num
-    do b = 1, mo_num
-      fock_3e_uhf_mo_a_old(b,a) += tmp(b,a)
-    enddo
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp)
-  !$OMP END PARALLEL
-
-  !call wall_time(tf)
-  !print *, ' Wall time for fock_3e_uhf_mo_a_old =', tf - ti
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b_old, (mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! BETA part of the Fock matrix from three-electron terms 
-  !
-  ! WARNING :: non hermitian if bi-ortho MOS used 
-  !
-  END_DOC
-
-  implicit none
-  integer                       :: a, b, i, j, o
-  double precision              :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
-  double precision              :: ti, tf
-  double precision, allocatable :: tmp(:,:)
-
-  PROVIDE mo_l_coef mo_r_coef
-
-  !print *, ' PROVIDING fock_3e_uhf_mo_b_old ...'
-  !call wall_time(ti)
-
-  o = elec_beta_num + 1
-  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
-
-  PROVIDE fock_3e_uhf_mo_cs_old
-  fock_3e_uhf_mo_b_old = fock_3e_uhf_mo_cs_old
-
-  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
-  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
-  !$OMP SHARED  (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_b_old)
-
-  allocate(tmp(mo_num,mo_num))
-  tmp = 0.d0
-
-  !$OMP DO
-  do a = 1, mo_num
-    do b = 1, mo_num
-
-      ! ---
-
-      do j = o, elec_alpha_num
-        do i = 1, elec_beta_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
-                              -        I_bij_aji &
-                              -        I_bij_iaj )
-
-        enddo
-      enddo
-
-      ! ---
-
-      do j = 1, elec_beta_num
-        do i = o, elec_alpha_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
-                              -        I_bij_aji &
-                              -        I_bij_jia )
-
-        enddo
-      enddo
-
-      ! ---
-
-      do j = o, elec_alpha_num
-        do i = o, elec_alpha_num
-
-          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
-          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
-          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
-          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
-          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
-          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
-
-          tmp(b,a) -= 0.5d0 * ( I_bij_aij &
-                              - I_bij_aji )
-
-        enddo
-      enddo
-
-      ! ---
-
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do a = 1, mo_num
-    do b = 1, mo_num
-      fock_3e_uhf_mo_b_old(b,a) += tmp(b,a)
-    enddo
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(tmp)
-  !$OMP END PARALLEL
-
-  !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_mo_b_old =', tf - ti
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)]
-
-  BEGIN_DOC
-  !
-  ! Equations (B6) and (B7)
-  !
-  ! g <--> gamma
-  ! d <--> delta
-  ! e <--> eta
-  ! k <--> kappa
-  !
-  END_DOC
-
-  implicit none
-  integer                       :: g, d, e, k, mu, nu
-  double precision              :: dm_ge_a, dm_ge_b, dm_ge
-  double precision              :: dm_dk_a, dm_dk_b, dm_dk
-  double precision              :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
-  double precision              :: ti, tf
-  double precision, allocatable :: f_tmp(:,:)
-
-  !print *, ' PROVIDING fock_3e_uhf_ao_a ...'
-  !call wall_time(ti)
-
-  fock_3e_uhf_ao_a = 0.d0
-
-  !$OMP PARALLEL DEFAULT (NONE)                                                                &
-  !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp,  &
-  !$OMP          i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
-  !$OMP SHARED  (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a)
-
-  allocate(f_tmp(ao_num,ao_num))
-  f_tmp = 0.d0
-
-  !$OMP DO
-  do g = 1, ao_num
-    do e = 1, ao_num
-      dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
-      dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
-      dm_ge   = dm_ge_a + dm_ge_b
-      do d = 1, ao_num
-        do k = 1, ao_num
-          dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
-          dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
-          dm_dk   = dm_dk_a + dm_dk_b
-          do mu = 1, ao_num
-            do nu = 1, ao_num
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
-              f_tmp(mu,nu) -= 0.5d0 * ( dm_ge   * dm_dk   * i_mugd_nuek &
-                                      + dm_ge_a * dm_dk_a * i_mugd_eknu &
-                                      + dm_ge_a * dm_dk_a * i_mugd_knue &
-                                      - dm_ge_a * dm_dk   * i_mugd_enuk &
-                                      - dm_ge   * dm_dk_a * i_mugd_kenu &
-                                      - dm_ge_a * dm_dk_a * i_mugd_nuke &
-                                      - dm_ge_b * dm_dk_b * i_mugd_nuke )
-            enddo
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do mu = 1, ao_num
-    do nu = 1, ao_num
-      fock_3e_uhf_ao_a(mu,nu) += f_tmp(mu,nu)
-    enddo
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(f_tmp)
-  !$OMP END PARALLEL
-
-  !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)]
-
-  BEGIN_DOC
-  !
-  ! Equations (B6) and (B7)
-  !
-  ! g <--> gamma
-  ! d <--> delta
-  ! e <--> eta
-  ! k <--> kappa
-  !
-  END_DOC
-
-  implicit none
-  integer                       :: g, d, e, k, mu, nu
-  double precision              :: dm_ge_a, dm_ge_b, dm_ge
-  double precision              :: dm_dk_a, dm_dk_b, dm_dk
-  double precision              :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
-  double precision              :: ti, tf
-  double precision, allocatable :: f_tmp(:,:)
-
-  !print *, ' PROVIDING fock_3e_uhf_ao_b ...'
-  !call wall_time(ti)
-
-  fock_3e_uhf_ao_b = 0.d0
-
-  !$OMP PARALLEL DEFAULT (NONE)                                                                &
-  !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp,  &
-  !$OMP          i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
-  !$OMP SHARED  (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b)
-
-  allocate(f_tmp(ao_num,ao_num))
-  f_tmp = 0.d0
-
-  !$OMP DO
-  do g = 1, ao_num
-    do e = 1, ao_num
-      dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
-      dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
-      dm_ge   = dm_ge_a + dm_ge_b
-      do d = 1, ao_num
-        do k = 1, ao_num
-          dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
-          dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
-          dm_dk   = dm_dk_a + dm_dk_b
-          do mu = 1, ao_num
-            do nu = 1, ao_num
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
-              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
-              f_tmp(mu,nu) -= 0.5d0 * ( dm_ge   * dm_dk   * i_mugd_nuek &
-                                      + dm_ge_b * dm_dk_b * i_mugd_eknu &
-                                      + dm_ge_b * dm_dk_b * i_mugd_knue &
-                                      - dm_ge_b * dm_dk   * i_mugd_enuk &
-                                      - dm_ge   * dm_dk_b * i_mugd_kenu &
-                                      - dm_ge_b * dm_dk_b * i_mugd_nuke &
-                                      - dm_ge_a * dm_dk_a * i_mugd_nuke )
-            enddo
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  !$OMP CRITICAL
-  do mu = 1, ao_num
-    do nu = 1, ao_num
-      fock_3e_uhf_ao_b(mu,nu) += f_tmp(mu,nu)
-    enddo
-  enddo
-  !$OMP END CRITICAL
-
-  deallocate(f_tmp)
-  !$OMP END PARALLEL
-
-  !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti
-
-END_PROVIDER 
-
-! ---
-

plugins/local/tc_scf/fock_tc_mo_tot.irp.f

@@ -1,4 +1,6 @@
 
+! ---
+
  BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ]
 &BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)]
 
@@ -23,9 +25,6 @@
   integer          :: i, j, n
   double precision :: t0, t1
 
-  !print*, ' Providing Fock_matrix_tc_mo_tot ...'
-  !call wall_time(t0)
-
   if(elec_alpha_num == elec_beta_num) then
 
     PROVIDE Fock_matrix_tc_mo_alpha
@@ -158,8 +157,8 @@
     Fock_matrix_tc_mo_tot += fock_3_mat
   endif
 
-  !call wall_time(t1)
-  !print*, ' Wall time for Fock_matrix_tc_mo_tot =', t1-t0
-
 END_PROVIDER
 
+! ---
+
+

plugins/local/tc_scf/fock_vartc.irp.f

@@ -1,287 +0,0 @@
-
-! ---
-
- BEGIN_PROVIDER [ double precision, two_e_vartc_integral_alpha, (ao_num, ao_num)]
-&BEGIN_PROVIDER [ double precision, two_e_vartc_integral_beta , (ao_num, ao_num)]
-
-  implicit none
-  integer                       :: i, j, k, l
-  double precision              :: density, density_a, density_b, I_coul, I_kjli
-  double precision              :: t0, t1
-  double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
-
-  two_e_vartc_integral_alpha = 0.d0
-  two_e_vartc_integral_beta  = 0.d0
-
- !$OMP PARALLEL DEFAULT (NONE)                                                                        &
- !$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli)              &
- !$OMP SHARED  (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_tc_tot, &
- !$OMP         two_e_vartc_integral_alpha, two_e_vartc_integral_beta)
-
-  allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num))
-  tmp_a = 0.d0
-  tmp_b = 0.d0
-
- !$OMP DO
-  do j = 1, ao_num
-    do l = 1, ao_num
-      density_a = TCSCF_density_matrix_ao_alpha(l,j)
-      density_b = TCSCF_density_matrix_ao_beta (l,j)
-      density   = density_a + density_b                      
-      do i = 1, ao_num
-        do k = 1, ao_num
-
-          I_coul = density * ao_two_e_tc_tot(k,i,l,j)
-          I_kjli = ao_two_e_tc_tot(k,j,l,i)
-
-          tmp_a(k,i) += I_coul - density_a * I_kjli
-          tmp_b(k,i) += I_coul - density_b * I_kjli
-        enddo
-      enddo
-    enddo
-  enddo
- !$OMP END DO NOWAIT
-
- !$OMP CRITICAL
-  do i = 1, ao_num
-    do j = 1, ao_num
-      two_e_vartc_integral_alpha(j,i) += tmp_a(j,i)
-      two_e_vartc_integral_beta (j,i) += tmp_b(j,i)
-    enddo
-  enddo
- !$OMP END CRITICAL
-
-  deallocate(tmp_a, tmp_b)
- !$OMP END PARALLEL
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_alpha, (ao_num, ao_num)]
-
-  implicit none
-
-  Fock_matrix_vartc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_vartc_integral_alpha 
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_beta, (ao_num, ao_num)]
-
-  implicit none
-
-  Fock_matrix_vartc_ao_beta = ao_one_e_integrals_tc_tot + two_e_vartc_integral_beta 
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_alpha, (mo_num, mo_num) ]
-
-  implicit none
-
-  call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_alpha, size(Fock_matrix_vartc_ao_alpha, 1) &
-                        , Fock_matrix_vartc_mo_alpha, size(Fock_matrix_vartc_mo_alpha, 1) )
-  if(three_body_h_tc) then
-    Fock_matrix_vartc_mo_alpha += fock_3e_uhf_mo_a
-  endif
-
-END_PROVIDER
-
-! ---
-
-BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_beta, (mo_num,mo_num) ]
-
-  implicit none
-
-  call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_beta, size(Fock_matrix_vartc_ao_beta, 1) &
-                        , Fock_matrix_vartc_mo_beta, size(Fock_matrix_vartc_mo_beta, 1) )
-  if(three_body_h_tc) then
-    Fock_matrix_vartc_mo_beta += fock_3e_uhf_mo_b
-  endif
-
-END_PROVIDER
-
-! ---
-
-BEGIN_PROVIDER [ double precision, grad_vartc]
-
-  implicit none
-  integer          :: i, k
-  double precision :: grad_left, grad_right
-
-  grad_left  = 0.d0
-  grad_right = 0.d0
-
-  do i = 1, elec_beta_num ! doc --> SOMO
-    do k = elec_beta_num+1, elec_alpha_num
-      grad_left  = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
-      grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
-    enddo
-  enddo
-
-  do i = 1, elec_beta_num ! doc --> virt 
-    do k = elec_alpha_num+1, mo_num
-      grad_left  = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
-      grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
-    enddo
-  enddo
-
-  do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt 
-    do k = elec_alpha_num+1, mo_num
-      grad_left  = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
-      grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
-    enddo
-  enddo
-
-  grad_vartc = grad_left + grad_right
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_tot, (ao_num, ao_num) ]
-
-  implicit none
-
-  call mo_to_ao_bi_ortho( Fock_matrix_vartc_mo_tot, size(Fock_matrix_vartc_mo_tot, 1) &
-                        , Fock_matrix_vartc_ao_tot, size(Fock_matrix_vartc_ao_tot, 1) )
-
-END_PROVIDER
-
-! ---
-
- BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_tot, (mo_num,mo_num) ]
-&BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_diag_mo_tot, (mo_num)]
-
-  implicit none
-  integer :: i, j, n
-
-  if(elec_alpha_num == elec_beta_num) then
-    Fock_matrix_vartc_mo_tot = Fock_matrix_vartc_mo_alpha
-  else
-
-    do j = 1, elec_beta_num
-      ! F-K
-      do i = 1, elec_beta_num !CC
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
-            - (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-      ! F+K/2
-      do i = elec_beta_num+1, elec_alpha_num  !CA
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-      ! F
-      do i = elec_alpha_num+1, mo_num !CV
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
-      enddo
-    enddo
-
-    do j = elec_beta_num+1, elec_alpha_num
-      ! F+K/2
-      do i = 1, elec_beta_num !AC
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-      ! F
-      do i = elec_beta_num+1, elec_alpha_num !AA
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
-      enddo
-      ! F-K/2
-      do i = elec_alpha_num+1, mo_num !AV
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-    enddo
-    
-    do j = elec_alpha_num+1, mo_num
-      ! F
-      do i = 1, elec_beta_num !VC
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
-      enddo
-      ! F-K/2
-      do i = elec_beta_num+1, elec_alpha_num !VA
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-      ! F+K
-      do i = elec_alpha_num+1, mo_num !VV
-         Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) &
-             + (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
-      enddo
-    enddo
-    if(three_body_h_tc)then
-      ! C-O
-      do j = 1, elec_beta_num
-        do i = elec_beta_num+1, elec_alpha_num
-          Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
-          Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
-        enddo
-      enddo
-      ! C-V
-      do j = 1, elec_beta_num
-        do i = elec_alpha_num+1, mo_num
-          Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
-          Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
-        enddo
-      enddo
-      ! O-V
-      do j = elec_beta_num+1, elec_alpha_num
-        do i = elec_alpha_num+1, mo_num
-          Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
-          Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
-        enddo
-      enddo
-    endif
-
-  endif
-
-  do i = 1, mo_num
-    Fock_matrix_vartc_diag_mo_tot(i) = Fock_matrix_vartc_mo_tot(i,i)
-  enddo
-
-  if(frozen_orb_scf)then
-    integer :: iorb, jorb
-    do i = 1, n_core_orb
-     iorb = list_core(i)
-     do j = 1, n_act_orb
-      jorb = list_act(j)
-      Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
-      Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
-     enddo
-    enddo
-  endif
-
-  if(no_oa_or_av_opt)then
-    do i = 1, n_act_orb
-      iorb = list_act(i)
-      do j = 1, n_inact_orb
-        jorb = list_inact(j)
-        Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
-        Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
-      enddo
-      do j = 1, n_virt_orb
-        jorb = list_virt(j)
-        Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
-        Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
-      enddo
-      do j = 1, n_core_orb
-        jorb = list_core(j)
-        Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
-        Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0                                                                                                                 
-      enddo
-    enddo
-  endif
-
-  !call check_sym(Fock_matrix_vartc_mo_tot, mo_num)
-  !do i = 1, mo_num
-  !  write(*,'(100(F15.8, I4))') Fock_matrix_vartc_mo_tot(i,:)
-  !enddo
-
-END_PROVIDER
-
-! ---
-

plugins/local/tc_scf/rh_tcscf_diis.irp.f

@@ -234,7 +234,7 @@ subroutine rh_tcscf_diis()
     call unlock_io
 
     if(er_delta .lt. 0.d0) then
-      call ezfio_set_tc_scf_bitc_energy(etc_tot)
+      call ezfio_set_tc_scf_tcscf_energy(etc_tot)
       call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
       call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
       write(json_unit, json_true_fmt) 'saved'
@@ -263,7 +263,7 @@ subroutine rh_tcscf_diis()
 
   deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, E_DIIS)
 
-  call ezfio_set_tc_scf_bitc_energy(TC_HF_energy)
+  call ezfio_set_tc_scf_tcscf_energy(TC_HF_energy)
   call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
   call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
 

plugins/local/tc_scf/rh_tcscf_simple.irp.f

@@ -91,7 +91,7 @@ subroutine rh_tcscf_simple()
       e_delta = dabs(etc_tot - e_save)
       e_save  = etc_tot
 
-      call ezfio_set_tc_scf_bitc_energy(etc_tot)
+      call ezfio_set_tc_scf_tcscf_energy(etc_tot)
 
       call wall_time(t1)
       write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &

plugins/local/tc_scf/rh_vartcscf_simple.irp.f

@@ -1,89 +0,0 @@
-! ---
-
-subroutine rh_vartcscf_simple()
-
-  implicit none
-  integer          :: i, j, it, dim_DIIS
-  double precision :: t0, t1
-  double precision :: e_save, e_delta, rho_delta
-  double precision :: etc_tot, etc_1e, etc_2e, etc_3e
-  double precision :: er_DIIS
-
-
-  it       = 0
-  e_save   = 0.d0
-  dim_DIIS = 0
-
-  ! ---
-
-  PROVIDE level_shift_tcscf
-  PROVIDE mo_r_coef
-
-  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')                      &
-    '====', '================', '================', '================', '================', '================' &
-          , '================', '================', '====', '========'
-  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')                      &
-    ' it ', '  SCF TC Energy ', '      E(1e)     ', '      E(2e)     ', '      E(3e)     ', '   energy diff  ' &
-          , '    DIIS error  ', '  level shift   ', 'DIIS', '  WT (m)'
-  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')                      &
-    '====', '================', '================', '================', '================', '================' &
-          , '================', '================', '====', '========'
-
-
-  ! first iteration (HF orbitals)
-  call wall_time(t0)
-
-  etc_tot = VARTC_HF_energy
-  etc_1e  = VARTC_HF_one_e_energy
-  etc_2e  = VARTC_HF_two_e_energy
-  etc_3e  = 0.d0
-  if(three_body_h_tc) then
-    etc_3e = diag_three_elem_hf
-  endif
-  er_DIIS = maxval(abs(FQS_SQF_mo))
-  e_delta = dabs(etc_tot - e_save)
-  e_save  = etc_tot
-
-  call wall_time(t1)
-  write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
-    it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
-
-  do while(er_DIIS .gt. dsqrt(thresh_tcscf))
-    call wall_time(t0)
-
-    it += 1
-    if(it > n_it_tcscf_max) then
-      print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
-      stop
-    endif
-
-    mo_r_coef = fock_vartc_eigvec_ao
-    mo_l_coef = mo_r_coef
-    call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
-    call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
-    TOUCH mo_l_coef mo_r_coef
-
-    etc_tot = VARTC_HF_energy
-    etc_1e  = VARTC_HF_one_e_energy
-    etc_2e  = VARTC_HF_two_e_energy
-    etc_3e  = 0.d0
-    if(three_body_h_tc) then
-      etc_3e = diag_three_elem_hf
-    endif
-    er_DIIS = maxval(abs(FQS_SQF_mo))
-    e_delta = dabs(etc_tot - e_save)
-    e_save  = etc_tot
-
-    call ezfio_set_tc_scf_bitc_energy(etc_tot)
-
-    call wall_time(t1)
-    write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
-      it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
-  enddo
-
-  print *, ' VAR-TCSCF Simple converged !'
-
-end
-
-! ---
-

plugins/local/tc_scf/tc_scf.irp.f

@@ -13,7 +13,6 @@ program tc_scf
   PROVIDE j1e_type
   PROVIDE j2e_type
   PROVIDE tcscf_algorithm
-  PROVIDE var_tc
 
   print *, ' TC-SCF with:'
   print *, ' j1e_type = ', j1e_type
@@ -45,22 +44,6 @@ program tc_scf
   !call create_guess()
   !call orthonormalize_mos()
 
-
-  if(var_tc) then
-
-    print *, ' VAR-TC'
-
-    if(tcscf_algorithm == 'DIIS') then
-      print*, ' NOT implemented yet'
-    elseif(tcscf_algorithm == 'Simple') then
-      call rh_vartcscf_simple()
-    else
-      print *, ' not implemented yet', tcscf_algorithm
-      stop
-    endif
-
-  else
-
   if(tcscf_algorithm == 'DIIS') then
     call rh_tcscf_diis()
   elseif(tcscf_algorithm == 'Simple') then
@@ -83,7 +66,6 @@ program tc_scf
   endif
   call print_energy_and_mos(good_angles)
 
-  endif
 
   write(json_unit,json_array_close_fmtx)
   call json_close

plugins/local/tc_scf/tc_scf_energy.irp.f

@@ -11,11 +11,8 @@
   integer          :: i, j
   double precision :: t0, t1
 
-  !print*, ' Providing TC energy ...'
-  !call wall_time(t0)
-
   PROVIDE mo_l_coef mo_r_coef
-  PROVIDE two_e_tc_non_hermit_integral_alpha two_e_tc_non_hermit_integral_beta
+  PROVIDE two_e_tc_integral_alpha two_e_tc_integral_beta
 
   TC_HF_energy = nuclear_repulsion
   TC_HF_one_e_energy = 0.d0
@@ -23,8 +20,8 @@
 
   do j = 1, ao_num
     do i = 1, ao_num
-      TC_HF_two_e_energy += 0.5d0 * ( two_e_tc_non_hermit_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
+      TC_HF_two_e_energy += 0.5d0 * ( two_e_tc_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
-                                    + two_e_tc_non_hermit_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) )
+                                    + two_e_tc_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) )
       TC_HF_one_e_energy += ao_one_e_integrals_tc_tot(i,j) &
                           * (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) )
     enddo
@@ -33,38 +30,6 @@
   TC_HF_energy += TC_HF_one_e_energy + TC_HF_two_e_energy
   TC_HF_energy += diag_three_elem_hf
 
-  !call wall_time(t1)
-  !print*, ' Wall time for TC energy=', t1-t0
-
-END_PROVIDER
-
-! ---
-
- BEGIN_PROVIDER [ double precision, VARTC_HF_energy]
-&BEGIN_PROVIDER [ double precision, VARTC_HF_one_e_energy]
-&BEGIN_PROVIDER [ double precision, VARTC_HF_two_e_energy]
-
-  implicit none
-  integer :: i, j
-
-  PROVIDE mo_r_coef
-
-  VARTC_HF_energy = nuclear_repulsion
-  VARTC_HF_one_e_energy = 0.d0
-  VARTC_HF_two_e_energy = 0.d0
-
-  do j = 1, ao_num
-    do i = 1, ao_num
-      VARTC_HF_two_e_energy += 0.5d0 * ( two_e_vartc_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
-                                       + two_e_vartc_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) )
-      VARTC_HF_one_e_energy += ao_one_e_integrals_tc_tot(i,j) &
-                             * (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) )
-    enddo
-  enddo
-
-  VARTC_HF_energy += VARTC_HF_one_e_energy + VARTC_HF_two_e_energy
-  VARTC_HF_energy += diag_three_elem_hf
-
 END_PROVIDER
 
 ! ---

plugins/local/tc_scf/test_int.irp.f

@@ -1,970 +0,0 @@
-program test_ints
-
-  BEGIN_DOC
-  ! TODO : Put the documentation of the program here
-  END_DOC
-
-  implicit none
-
-  print *, ' starting test_ints ...'
-
-  my_grid_becke  = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  my_extra_grid_becke = .True.
-  my_n_pt_r_extra_grid = 30
-  my_n_pt_a_extra_grid = 50 ! small extra_grid for quick debug
-  touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
-
-!! OK 
-! call routine_int2_u_grad1u_env2 
-! OK
-! call routine_v_ij_erf_rk_cst_mu_env
-! OK 
-! call routine_x_v_ij_erf_rk_cst_mu_env
-! OK
-! call routine_int2_u2_env2
-! OK
-! call routine_int2_u_grad1u_x_env2
-! OK 
-! call routine_int2_grad1u2_grad2u2_env2
-! call routine_int2_u_grad1u_env2
-! call test_int2_grad1_u12_ao_test
-! call routine_v_ij_u_cst_mu_env_test
-! call test_grid_points_ao
- !call test_int_gauss
-
-  !call test_fock_3e_uhf_ao()
-  !call test_fock_3e_uhf_mo()
-
-  !call test_two_e_tc_non_hermit_integral()
-
-!!PROVIDE TC_HF_energy VARTC_HF_energy
-!!print *, '    TC_HF_energy = ',    TC_HF_energy
-!!print *, ' VARTC_HF_energy = ', VARTC_HF_energy
-
-  call test_fock_3e_uhf_mo_cs()
-  call test_fock_3e_uhf_mo_a()
-  call test_fock_3e_uhf_mo_b()
-
-end
-
-! ---
-
-subroutine routine_test_env
- implicit none
- integer :: i,icount,j
- icount = 0
- do i = 1, List_env1s_square_size
-  if(dabs(List_env1s_square_coef(i)).gt.1.d-10)then
-   print*,''
-   print*,List_env1s_square_expo(i),List_env1s_square_coef(i)
-   print*,List_env1s_square_cent(1:3,i)
-   print*,''
-   icount += 1
-  endif
-  
- enddo
- print*,'List_env1s_square_coef,icount = ',List_env1s_square_size,icount
- do i = 1, ao_num
-  do j = 1, ao_num
-   do icount = 1, List_comb_thr_b3_size(j,i)
-    print*,'',j,i
-    print*,List_comb_thr_b3_expo(icount,j,i),List_comb_thr_b3_coef(icount,j,i)
-    print*,List_comb_thr_b3_cent(1:3,icount,j,i)
-    print*,''
-   enddo
-!   enddo
-  enddo
- enddo
- print*,'max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size,List_env1s_square_size
-
-end
-
-subroutine routine_int2_u_grad1u_env2
- implicit none
- integer :: i,j,ipoint,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += int2_u_grad1u_env2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += int2_u_grad1u_env2(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_int2_u_grad1u_env2'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-subroutine routine_v_ij_erf_rk_cst_mu_env
- implicit none
- integer :: i,j,ipoint,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += v_ij_erf_rk_cst_mu_env_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += v_ij_erf_rk_cst_mu_env(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_v_ij_erf_rk_cst_mu_env'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-
-subroutine routine_x_v_ij_erf_rk_cst_mu_env
- implicit none
- integer :: i,j,ipoint,k,l,m
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      do m = 1, 3
-       array(j,i,l,k)     += x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-       array_ref(j,i,l,k) += x_v_ij_erf_rk_cst_mu_env     (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      enddo
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
-
- print*,'******'
- print*,'******'
- print*,'routine_x_v_ij_erf_rk_cst_mu_env'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-
-
-subroutine routine_v_ij_u_cst_mu_env_test
- implicit none
- integer :: i,j,ipoint,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += v_ij_u_cst_mu_env_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += v_ij_u_cst_mu_env_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_v_ij_u_cst_mu_env_test'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-end
-
-subroutine routine_int2_grad1u2_grad2u2_env2
- implicit none
- integer :: i,j,ipoint,k,l
- integer :: ii , jj
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
- double precision, allocatable :: ints(:,:,:)
- allocate(ints(ao_num, ao_num, n_points_final_grid))
-! do ipoint = 1, n_points_final_grid
-!  do i = 1, ao_num
-!   do j = 1, ao_num
-!    read(33,*)ints(j,i,ipoint)
-!   enddo
-!  enddo
-! enddo
-
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += int2_grad1u2_grad2u2_env2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-!     !array(j,i,l,k)     += int2_grad1u2_grad2u2_env2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-!      array_ref(j,i,l,k)     += int2_grad1u2_grad2u2_env2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-!     !array(j,i,l,k) += ints(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-!       array_ref(j,i,l,k) += int2_grad1u2_grad2u2_env2(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += ints(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-!      if(dabs(int2_grad1u2_grad2u2_env2_test(j,i,ipoint)).gt.1.d-6)then
-!       if(dabs(int2_grad1u2_grad2u2_env2_test(j,i,ipoint) - int2_grad1u2_grad2u2_env2_test(j,i,ipoint)).gt.1.d-6)then
-!        print*,j,i,ipoint
-!        print*,int2_grad1u2_grad2u2_env2_test(j,i,ipoint) , int2_grad1u2_grad2u2_env2_test(j,i,ipoint), dabs(int2_grad1u2_grad2u2_env2_test(j,i,ipoint) - int2_grad1u2_grad2u2_env2_test(j,i,ipoint))
-!        print*,int2_grad1u2_grad2u2_env2_test(i,j,ipoint) , int2_grad1u2_grad2u2_env2_test(i,j,ipoint), dabs(int2_grad1u2_grad2u2_env2_test(i,j,ipoint) - int2_grad1u2_grad2u2_env2_test(i,j,ipoint))
-!        stop
-!       endif
-!      endif
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- double precision :: e_ref, e_new
- accu_relat = 0.d0
- accu_abs   = 0.d0
- e_ref = 0.d0
- e_new = 0.d0
- do ii = 1, elec_alpha_num
-  do jj = ii, elec_alpha_num
-   do k = 1, ao_num
-    do l = 1, ao_num
-     do i = 1, ao_num
-      do j = 1, ao_num
-       e_ref += mo_coef(j,ii) * mo_coef(i,ii) * array_ref(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj)
-       e_new += mo_coef(j,ii) * mo_coef(i,ii) * array(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj)
-       contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-       accu_abs += contrib
-!       if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-!        accu_relat += contrib/dabs(array_ref(j,i,l,k))
-!       endif
-      enddo
-     enddo
-    enddo
-   enddo
-
-  enddo
- enddo
- print*,'e_ref = ',e_ref
- print*,'e_new = ',e_new
-! print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
-! print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-subroutine routine_int2_u2_env2
- implicit none
- integer :: i,j,ipoint,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += int2_u2_env2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += int2_u2_env2(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_int2_u2_env2'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-
-subroutine routine_int2_u_grad1u_x_env2
- implicit none
- integer :: i,j,ipoint,k,l,m
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      do m = 1, 3
-       array(j,i,l,k)     += int2_u_grad1u_x_env2_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-       array_ref(j,i,l,k) += int2_u_grad1u_x_env2     (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      enddo
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_int2_u_grad1u_x_env2'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-  
-
-end
-
-subroutine routine_v_ij_u_cst_mu_env
- implicit none
- integer :: i,j,ipoint,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do ipoint = 1, n_points_final_grid
-  weight = final_weight_at_r_vector(ipoint)
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      array(j,i,l,k)     += v_ij_u_cst_mu_env_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += v_ij_u_cst_mu_env_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'routine_v_ij_u_cst_mu_env'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-
-end
-
-! ---
-
-subroutine test_fock_3e_uhf_ao()
-
-  implicit none
-  integer                       :: i, j
-  double precision              :: diff_tot, diff_ij, thr_ih, norm
-  double precision, allocatable :: fock_3e_uhf_ao_a_mo(:,:), fock_3e_uhf_ao_b_mo(:,:)
-
-  thr_ih = 1d-7
-
-  PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth
-  PROVIDE fock_3e_uhf_ao_a fock_3e_uhf_ao_b
-
-  ! ---
-
-  allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num))
-  call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a   , size(fock_3e_uhf_ao_a   , 1) &
-                        , fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) )
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      diff_ij = dabs(fock_3e_uhf_ao_a_mo(j,i) - fock_a_tot_3e_bi_orth(j,i))
-      if(diff_ij .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i)
-        print *, ' UHF  : ', fock_3e_uhf_ao_a_mo  (j,i)
-        !stop
-      endif
-
-      norm     += dabs(fock_a_tot_3e_bi_orth(j,i))
-      diff_tot += diff_ij
-    enddo
-  enddo
-  print *, ' diff on F_a = ', diff_tot / norm
-  print *, ' '
-
-  deallocate(fock_3e_uhf_ao_a_mo)
-
-  ! ---
-
-  allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num))
-  call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b   , size(fock_3e_uhf_ao_b   , 1) &
-                        , fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) )
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      diff_ij = dabs(fock_3e_uhf_ao_b_mo(j,i) - fock_b_tot_3e_bi_orth(j,i))
-      if(diff_ij .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i)
-        print *, ' UHF  : ', fock_3e_uhf_ao_b_mo  (j,i)
-        !stop
-      endif
-
-      norm     += dabs(fock_b_tot_3e_bi_orth(j,i))
-      diff_tot += diff_ij
-    enddo
-  enddo
-  print *, ' diff on F_b = ', diff_tot/norm
-  print *, ' '
-
-  deallocate(fock_3e_uhf_ao_b_mo)
-
-  ! ---
-
-end subroutine test_fock_3e_uhf_ao()
-
-! ---
-
-subroutine test_fock_3e_uhf_mo()
-
-  implicit none
-  integer          :: i, j
-  double precision :: diff_tot, diff_ij, thr_ih, norm
-
-  thr_ih = 1d-12
-
-  PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth
-  PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b
-
-  ! ---
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i))
-      if(diff_ij .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i)
-        print *, ' UHF  : ', fock_3e_uhf_mo_a     (j,i)
-        !stop
-      endif
-
-      norm     += dabs(fock_a_tot_3e_bi_orth(j,i))
-      diff_tot += diff_ij
-    enddo
-  enddo
-  print *, ' diff on F_a = ', diff_tot / norm
-  print *, '      norm_a = ', norm
-  print *, ' '
-
-  ! ---
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i))
-      if(diff_ij .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i)
-        print *, ' UHF  : ', fock_3e_uhf_mo_b     (j,i)
-        !stop
-      endif
-
-      norm     += dabs(fock_b_tot_3e_bi_orth(j,i))
-      diff_tot += diff_ij
-    enddo
-  enddo
-  print *, ' diff on F_b = ', diff_tot/norm
-  print *, '      norm_b = ', norm
-  print *, ' '
-
-  ! ---
-
-end
-
-! ---
-
-subroutine test_grid_points_ao
- implicit none
- integer :: i,j,ipoint,icount,icount_good, icount_bad,icount_full
- double precision :: thr
- thr = 1.d-10
-! print*,'max_n_pts_grid_ao_prod = ',max_n_pts_grid_ao_prod
-! print*,'n_pts_grid_ao_prod'
- do i = 1, ao_num
-  do j = i, ao_num
-  icount = 0
-  icount_good = 0
-  icount_bad = 0
-  icount_full = 0
-  do ipoint = 1, n_points_final_grid
-!   if(dabs(int2_u_grad1u_x_env2_test(j,i,ipoint,1)) & 
-!    + dabs(int2_u_grad1u_x_env2_test(j,i,ipoint,2)) &
-!    + dabs(int2_u_grad1u_x_env2_test(j,i,ipoint,3)) )
-!   if(dabs(int2_u2_env2_test(j,i,ipoint)).gt.thr)then
-!    icount += 1
-!   endif
-   if(dabs(v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then
-    icount_full += 1
-   endif
-   if(dabs(v_ij_u_cst_mu_env_test(j,i,ipoint)).gt.thr)then
-    icount += 1
-    if(dabs(v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then
-    icount_good += 1
-    else
-    print*,j,i,ipoint
-    print*,dabs(v_ij_u_cst_mu_env_test(j,i,ipoint)), dabs(v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint))/dabs(v_ij_u_cst_mu_env_test(j,i,ipoint))
-    icount_bad  += 1
-    endif
-   endif
-!   if(dabs(v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint)).gt.thr)then
-!   endif
-  enddo
-   print*,''
-   print*,j,i
-   print*,icount,icount_full, icount_bad!,n_pts_grid_ao_prod(j,i)
-   print*,dble(icount)/dble(n_points_final_grid),dble(icount_full)/dble(n_points_final_grid)
-!          dble(n_pts_grid_ao_prod(j,i))/dble(n_points_final_grid)
-!   if(icount.gt.n_pts_grid_ao_prod(j,i))then
-!    print*,'pb !!'
-!   endif
-  enddo
- enddo
-end
-
-subroutine test_int_gauss
- implicit none
- integer :: i,j
- print*,'center'
- do i = 1, ao_num
-  do j = i, ao_num
-   print*,j,i
-   print*,ao_prod_sigma(j,i),ao_overlap_abs_grid(j,i)
-   print*,ao_prod_center(1:3,j,i)
-  enddo
- enddo
- print*,''
- double precision :: weight, r(3),integral_1,pi,center(3),f_r,alpha,distance,integral_2
- center = 0.d0
- pi = dacos(-1.d0)
- integral_1 = 0.d0
- integral_2 = 0.d0
- alpha = 0.75d0
- do i = 1,  n_points_final_grid
-  ! you get x, y and z of the ith grid point
-  r(1) = final_grid_points(1,i)
-  r(2) = final_grid_points(2,i)
-  r(3) = final_grid_points(3,i)
-  weight = final_weight_at_r_vector(i)
-  distance = dsqrt( (r(1) - center(1))**2 +  (r(2) - center(2))**2 + (r(3) - center(3))**2 )
-  f_r = dexp(-alpha * distance*distance)
-  ! you add the contribution of the grid point to the integral
-  integral_1 += f_r * weight
-  integral_2 += f_r * distance * weight
- enddo
- print*,'integral_1      =',integral_1
- print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5
- print*,'integral_2      =',integral_2
- print*,'(pi/alpha)**1.5 =',2.d0*pi / (alpha)**2
-
-
-end
-
-! ---
-
-subroutine test_two_e_tc_non_hermit_integral()
-
-  implicit none
-  integer          :: i, j
-  double precision :: diff_tot, diff, thr_ih, norm
-
-  thr_ih = 1d-10
-
-  PROVIDE two_e_tc_non_hermit_integral_beta two_e_tc_non_hermit_integral_alpha
-  PROVIDE two_e_tc_non_hermit_integral_seq_beta two_e_tc_non_hermit_integral_seq_alpha
-
-  ! ---
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, ao_num
-    do j = 1, ao_num
-
-      diff = dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i) - two_e_tc_non_hermit_integral_alpha(j,i))
-      if(diff .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' seq         : ', two_e_tc_non_hermit_integral_seq_alpha(j,i)
-        print *, ' //          : ', two_e_tc_non_hermit_integral_alpha    (j,i)
-        !stop
-      endif
-
-      norm     += dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i))
-      diff_tot += diff
-    enddo
-  enddo
-
-  print *, ' diff tot a = ', diff_tot / norm
-  print *, '     norm a = ', norm
-  print *, ' '
-
-  ! ---
-
-  norm     = 0.d0
-  diff_tot = 0.d0
-  do i = 1, ao_num
-    do j = 1, ao_num
-
-      diff = dabs(two_e_tc_non_hermit_integral_seq_beta(j,i) - two_e_tc_non_hermit_integral_beta(j,i))
-      if(diff .gt. thr_ih) then
-        print *, ' difference on ', j, i
-        print *, ' seq         : ', two_e_tc_non_hermit_integral_seq_beta(j,i)
-        print *, ' //          : ', two_e_tc_non_hermit_integral_beta    (j,i)
-        !stop
-      endif
-
-      norm     += dabs(two_e_tc_non_hermit_integral_seq_beta(j,i))
-      diff_tot += diff
-    enddo
-  enddo
-
-  print *, ' diff tot b = ', diff_tot / norm
-  print *, '     norm b = ', norm
-  print *, ' '
-
-  ! ---
-
-  return
-
-end
-
-! ---
-
-subroutine test_int2_grad1_u12_ao_test
- implicit none
- integer :: i,j,ipoint,m,k,l
- double precision :: weight,accu_relat, accu_abs, contrib
- double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
- allocate(array(ao_num, ao_num, ao_num, ao_num))
- array = 0.d0
- allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
- array_ref = 0.d0
- do m = 1, 3
-  do ipoint = 1, n_points_final_grid
-   weight = final_weight_at_r_vector(ipoint)
-   do k = 1, ao_num
-    do l = 1, ao_num
-     do i = 1, ao_num
-      do j = 1, ao_num
-       array(j,i,l,k)     += int2_grad1_u12_ao_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-       array_ref(j,i,l,k) += int2_grad1_u12_ao(j,i,ipoint,m)      * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      enddo
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
-
- accu_relat = 0.d0
- accu_abs   = 0.d0
-  do k = 1, ao_num
-   do l = 1, ao_num
-    do i = 1, ao_num
-     do j = 1, ao_num
-      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
-      accu_abs += contrib
-      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
-       accu_relat += contrib/dabs(array_ref(j,i,l,k))
-      endif
-     enddo
-    enddo
-   enddo
-  enddo
- print*,'******'
- print*,'******'
- print*,'test_int2_grad1_u12_ao_test'
- print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
- print*,'accu_relat = ',accu_relat/dble(ao_num)**4
-end
-
-! ---
-
-subroutine test_fock_3e_uhf_mo_cs()
-
-  implicit none
-  integer          :: i, j
-  double precision :: I_old, I_new
-  double precision :: diff_tot, diff, thr_ih, norm
-
-!  double precision :: t0, t1
-!  print*, ' Providing fock_a_tot_3e_bi_orth ...'
-!  call wall_time(t0)
-!  PROVIDE fock_a_tot_3e_bi_orth
-!  call wall_time(t1)
-!  print*, ' Wall time for fock_a_tot_3e_bi_orth =', t1 - t0 
-
-  PROVIDE fock_3e_uhf_mo_cs fock_3e_uhf_mo_cs_old
-
-  thr_ih   = 1d-8
-  norm     = 0.d0
-  diff_tot = 0.d0
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      I_old = fock_3e_uhf_mo_cs_old(j,i)
-      I_new = fock_3e_uhf_mo_cs    (j,i)
-
-      diff = dabs(I_old - I_new)
-      if(diff .gt. thr_ih) then
-        print *, ' problem in fock_3e_uhf_mo_cs on ', j, i
-        print *, ' old value = ', I_old
-        print *, ' new value = ', I_new
-        !stop
-      endif
-
-      norm     += dabs(I_old)
-      diff_tot += diff
-    enddo
-  enddo
-
-  print *, ' diff tot (%) = ', 100.d0 * diff_tot / norm
-
-  return
-end
-
-! ---
-
-subroutine test_fock_3e_uhf_mo_a()
-
-  implicit none
-  integer          :: i, j
-  double precision :: I_old, I_new
-  double precision :: diff_tot, diff, thr_ih, norm
-
-  PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_a_old
-
-  thr_ih   = 1d-8
-  norm     = 0.d0
-  diff_tot = 0.d0
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      I_old = fock_3e_uhf_mo_a_old(j,i)
-      I_new = fock_3e_uhf_mo_a    (j,i)
-
-      diff = dabs(I_old - I_new)
-      if(diff .gt. thr_ih) then
-        print *, ' problem in fock_3e_uhf_mo_a on ', j, i
-        print *, ' old value = ', I_old
-        print *, ' new value = ', I_new
-        !stop
-      endif
-
-      norm     += dabs(I_old)
-      diff_tot += diff
-    enddo
-  enddo
-
-  print *, ' diff tot (%) = ', 100.d0 * diff_tot / norm
-
-  return
-end
-
-! ---
-
-subroutine test_fock_3e_uhf_mo_b()
-
-  implicit none
-  integer          :: i, j
-  double precision :: I_old, I_new
-  double precision :: diff_tot, diff, thr_ih, norm
-
-  PROVIDE fock_3e_uhf_mo_b fock_3e_uhf_mo_b_old
-
-  thr_ih   = 1d-8
-  norm     = 0.d0
-  diff_tot = 0.d0
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-
-      I_old = fock_3e_uhf_mo_b_old(j,i)
-      I_new = fock_3e_uhf_mo_b    (j,i)
-
-      diff = dabs(I_old - I_new)
-      if(diff .gt. thr_ih) then
-        print *, ' problem in fock_3e_uhf_mo_b on ', j, i
-        print *, ' old value = ', I_old
-        print *, ' new value = ', I_new
-        !stop
-      endif
-
-      norm     += dabs(I_old)
-      diff_tot += diff
-    enddo
-  enddo
-
-  print *, ' diff tot (%) = ', 100.d0 * diff_tot / norm
-
-  return
-end
-
-! ---
-

src/becke_numerical_grid/extra_grid_vector.irp.f

@@ -70,17 +70,6 @@ END_PROVIDER
         index_final_points_extra(2,i_count) = i
         index_final_points_extra(3,i_count) = j
         index_final_points_extra_reverse(k,i,j) = i_count
-
-        if(final_weight_at_r_vector_extra(i_count) .lt. 0.d0) then
-          print *, ' !!! WARNING !!!'
-          print *, ' negative weight !!!!'
-          print *, i_count, final_weight_at_r_vector_extra(i_count)
-          if(dabs(final_weight_at_r_vector_extra(i_count)) .lt. 1d-10) then
-            final_weight_at_r_vector_extra(i_count) = 0.d0
-          else
-            stop
-          endif
-        endif
       enddo
     enddo
   enddo

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -67,17 +67,6 @@ END_PROVIDER
         index_final_points(2,i_count) = i
         index_final_points(3,i_count) = j
         index_final_points_reverse(k,i,j) = i_count
-
-        if(final_weight_at_r_vector(i_count) .lt. 0.d0) then
-          print *, ' !!! WARNING !!!'
-          print *, ' negative weight !!!!'
-          print *, i_count, final_weight_at_r_vector(i_count)
-          if(dabs(final_weight_at_r_vector(i_count)) .lt. 1d-10) then
-            final_weight_at_r_vector(i_count) = 0.d0
-          else
-            stop
-          endif
-        endif 
       enddo
     enddo
   enddo

src/utils/util.irp.f

@@ -576,7 +576,7 @@ logical function is_same_spin(sigma_1, sigma_2)
     is_same_spin = .false.
   endif
 
-end function is_same_spin
+end
 
 ! ---
    
@@ -596,7 +596,7 @@ function Kronecker_delta(i, j) result(delta)
     delta = 0.d0
   endif
 
-end function Kronecker_delta
+end
 
 ! ---
 
@@ -634,7 +634,81 @@ subroutine diagonalize_sym_matrix(N, A, e)
     print*,'Problem in diagonalize_sym_matrix (dsyev)!!'
   endif
 
-end subroutine diagonalize_sym_matrix
+end
+
+! ---
+
+
+subroutine give_degen(A, n, shift, list_degen, n_degen_list)
+
+  BEGIN_DOC
+  ! returns n_degen_list :: the number of degenerated SET of elements (i.e. with |A(i)-A(i+1)| below shift)
+  !
+  ! for each of these sets, list_degen(1,i) = first degenerate element of the set i, 
+  !
+  !                         list_degen(2,i) = last degenerate element of the set i.
+  END_DOC
+
+  implicit none
+
+  double precision, intent(in)  :: A(n)
+  double precision, intent(in)  :: shift
+  integer,          intent(in)  :: n
+  integer,          intent(out) :: list_degen(2,n), n_degen_list
+
+  integer                       :: i, j, n_degen, k
+  logical                       :: keep_on
+  double precision, allocatable :: Aw(:)
+
+  list_degen = -1
+  allocate(Aw(n))
+  Aw = A
+  i=1
+  k = 0
+  do while(i.lt.n)
+   if(dabs(Aw(i)-Aw(i+1)).lt.shift)then
+    k+=1
+    j=1
+    list_degen(1,k) = i
+    keep_on = .True.
+    do while(keep_on)
+     if(i+j.gt.n)then
+      keep_on = .False.
+      exit
+     endif
+     if(dabs(Aw(i)-Aw(i+j)).lt.shift)then
+      j+=1
+     else
+      keep_on=.False.
+      exit
+     endif
+    enddo
+    n_degen = j
+    list_degen(2,k) = list_degen(1,k)-1 + n_degen
+    j=0
+    keep_on = .True.
+    do while(keep_on)
+     if(i+j+1.gt.n)then
+      keep_on = .False.
+      exit
+     endif
+     if(dabs(Aw(i+j)-Aw(i+j+1)).lt.shift)then
+      Aw(i+j) += (j-n_degen/2) * shift
+      j+=1
+     else
+      keep_on = .False.
+      exit
+     endif
+    enddo
+    Aw(i+n_degen-1) += (n_degen-1-n_degen/2) * shift
+    i+=n_degen
+   else
+    i+=1
+   endif
+  enddo
+  n_degen_list = k
+
+end
 
 ! ---