Merge branch 'dev-stable' into dev-stable-tc-scf

src/bi_ort_ints/bi_ort_ints.irp.f

@@ -20,6 +20,8 @@ program bi_ort_ints
 ! call test_5idx2
  !call test_4idx
   call test_4idx2()
+  call test_5idx2
+ call test_5idx
 end
 
 subroutine test_5idx2
@@ -75,6 +77,8 @@ subroutine test_5idx
  k = 1
  n = 0
  accu = 0.d0
+ PROVIDE three_e_5_idx_direct_bi_ort_old
+
  do i = 1, mo_num
   do k = 1, mo_num
    do j = 1, mo_num
@@ -84,28 +88,28 @@ subroutine test_5idx
 !         stop
 !      endif
 
-!      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
+      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
-!      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
+      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
+      contrib = dabs(new - ref)
+      accu += contrib
+      if(contrib .gt. 1.d-10)then
+       print*,'direct'
+       print*,i,k,j,l,m
+       print*,ref,new,contrib
+       stop
+      endif
+!
+!      new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
 !      contrib = dabs(new - ref)
 !      accu += contrib
 !      if(contrib .gt. 1.d-10)then
-!       print*,'direct'
+!       print*,'exch12'
 !       print*,i,k,j,l,m
 !       print*,ref,new,contrib
 !       stop
 !      endif
 !
-      new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'exch12'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
 !
 !      new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
 !      ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
@@ -150,7 +154,7 @@ subroutine test_5idx
 !       print*,ref,new,contrib
 !       stop
 !      endif
-!
+
      enddo
     enddo
    enddo

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -15,90 +15,108 @@ end
   !
   ! matrix element of the -L  three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
   !
-  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
+  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> :: : notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
   !
   END_DOC
 
   implicit none
-  integer                       :: i, j, k, m, l
-  integer                       :: ipoint
-  double precision              :: wall1, wall0
-  double precision, allocatable :: grad_mli(:,:,:), orb_mat(:,:,:)
-  double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:)
-  double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:)
-  double precision, allocatable :: tmp_mat(:,:,:,:)
 
-  allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
+  integer                        :: i, j, k, m, l
+  double precision               :: wall1, wall0
+  integer                        :: ipoint
+  double precision, allocatable  :: grad_mli(:,:), orb_mat(:,:,:)
+  double precision, allocatable  :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:)
+  double precision, allocatable  :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:)
+  double precision, allocatable  :: tmp_mat(:,:,:)
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
   PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
 
+  call print_memory_usage
   print *, ' Providing the three_e_5_idx_bi_ort ...'
   call wall_time(wall0)
 
-  do m = 1, mo_num
+  three_e_5_idx_direct_bi_ort (:,:,:,:,:) = 0.d0
+  three_e_5_idx_cycle_1_bi_ort(:,:,:,:,:) = 0.d0
+  three_e_5_idx_cycle_2_bi_ort(:,:,:,:,:) = 0.d0
+  three_e_5_idx_exch23_bi_ort (:,:,:,:,:) = 0.d0
+  three_e_5_idx_exch13_bi_ort (:,:,:,:,:) = 0.d0
 
-    allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
+  call print_memory_usage
-    allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
 
-    !$OMP PARALLEL                                                  &
+  allocate(tmp_mat(mo_num,mo_num,mo_num))
-    !$OMP DEFAULT (NONE)                                            &
+  allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
-    !$OMP PRIVATE (i,l,ipoint)                                      &
-    !$OMP SHARED (m,mo_num,n_points_final_grid,                     &
-    !$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         grad_mli, orb_mat)
-    !$OMP DO COLLAPSE(2)
-    do i = 1, mo_num
-      do l = 1, mo_num
-         do ipoint = 1, n_points_final_grid
 
-           grad_mli(ipoint,l,i) = final_weight_at_r_vector(ipoint) * ( &
+  !$OMP PARALLEL DO PRIVATE (i,l,ipoint)
-                 int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + &
+  do i=1,mo_num
-                 int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + &
+    do l=1,mo_num
-                 int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) )
+      do ipoint=1, n_points_final_grid
 
-           orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
+        orb_mat(ipoint,l,i) = final_weight_at_r_vector(ipoint)       &
+            * mos_l_in_r_array_transp(ipoint,l)                      &
+            * mos_r_in_r_array_transp(ipoint,i)
 
-         enddo
       enddo
     enddo
-    !$OMP END DO
+  enddo
-    !$OMP END PARALLEL
+  !$OMP END PARALLEL DO
 
-    call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
+  tmp_mat = 0.d0
-        orb_mat, n_points_final_grid,  &
+  call print_memory_usage
-        grad_mli, n_points_final_grid,  0.d0, &
+!
-        tmp_mat, mo_num*mo_num)
+  do m = 1, mo_num
 
-    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
+    allocate(grad_mli(n_points_final_grid,mo_num))
-    do i = 1, mo_num
+
+    do i=1,mo_num
+      !$OMP PARALLEL DO PRIVATE (l,ipoint)
+      do l=1,mo_num
+        do ipoint=1, n_points_final_grid
+
+          grad_mli(ipoint,l) =                                       &
+              int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) +&
+              int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) +&
+              int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
+
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, n_points_final_grid, 1.d0,&
+          orb_mat, n_points_final_grid,                              &
+          grad_mli, n_points_final_grid,  0.d0,                      &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL PRIVATE(j,k,l)
+      !$OMP DO
       do k = 1, mo_num
         do j = 1, mo_num
           do l = 1, mo_num
-            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
+            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k)
           enddo
         enddo
       enddo
+      !$OMP END DO
+      !$OMP DO
+      do j = 1, mo_num
+        do l = 1, mo_num
+          do k = 1, mo_num
+            three_e_5_idx_direct_bi_ort(m,k,i,l,j) = three_e_5_idx_direct_bi_ort(m,k,i,l,j) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
     enddo
-    !$OMP END PARALLEL DO
 
-    deallocate(orb_mat,grad_mli)
+    deallocate(grad_mli)
 
     allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
-    allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
+    allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
-    allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
 
-    !$OMP PARALLEL                                                  &
+    !$OMP PARALLEL DO PRIVATE (i,l,ipoint)
-    !$OMP DEFAULT (NONE)                                            &
-    !$OMP PRIVATE (i,l,ipoint)                                      &
-    !$OMP SHARED (m,mo_num,n_points_final_grid,                     &
-    !$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         rm_grad_ik, lm_grad_ik, rk_grad_im, lk_grad_mi)
-    !$OMP DO COLLAPSE(2)
     do i=1,mo_num
       do l=1,mo_num
         do ipoint=1, n_points_final_grid
@@ -107,128 +125,118 @@ end
           lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
           lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
 
-          rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
-          rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
-          rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
-
-          rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
-          rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
-          rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
-        enddo
-      enddo
-    enddo
-    !$OMP END DO
-    !$OMP END PARALLEL
-    call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
-        lm_grad_ik, 3*n_points_final_grid,        &
-        rm_grad_ik, 3*n_points_final_grid,  0.d0, &
-        tmp_mat, mo_num*mo_num)
-
-    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-    do i = 1, mo_num
-      do k = 1, mo_num
-        do j = 1, mo_num
-          do l = 1, mo_num
-            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
-          enddo
-        enddo
-      enddo
-    enddo
-    !$OMP END PARALLEL DO
-  
-    call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
-        lm_grad_ik, 3*n_points_final_grid,        &
-        rk_grad_im, 3*n_points_final_grid,  0.d0, &
-        tmp_mat, mo_num*mo_num)
-  
-    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-    do i = 1, mo_num
-      do k = 1, mo_num
-        do j = 1, mo_num
-          do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,j,k)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = - tmp_mat(k,j,i,l)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = - tmp_mat(k,i,j,l)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = - tmp_mat(l,j,i,k)
-          enddo
-        enddo
-      enddo
-    enddo
-    !$OMP END PARALLEL DO
-  
-    deallocate(lm_grad_ik)
-  
-    allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
-  
-   !$OMP PARALLEL                                                  &
-   !$OMP DEFAULT (NONE)                                            &
-   !$OMP PRIVATE (i,l,ipoint)                                      &
-   !$OMP SHARED (m,mo_num,n_points_final_grid,                     &
-   !$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         lk_grad_mi)
-   !$OMP DO COLLAPSE(2)
-    do i=1,mo_num
-      do l=1,mo_num
-        do ipoint=1, n_points_final_grid
-  
           lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint)
           lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint)
           lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint)
-  
+
-        enddo
-      enddo
-    enddo
-    !$OMP END DO
-    !$OMP END PARALLEL
-  
-    call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
-        lk_grad_mi, 3*n_points_final_grid,  &
-        rm_grad_ik, 3*n_points_final_grid,  0.d0, &
-        tmp_mat, mo_num*mo_num)
-  
-    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-    do i = 1, mo_num
-      do k = 1, mo_num
-        do j = 1, mo_num
-          do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l,i)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k,i)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(k,i,l,j)
-          enddo
         enddo
       enddo
     enddo
     !$OMP END PARALLEL DO
-  
+
-    call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
+    allocate(rm_grad_ik(n_points_final_grid,3,mo_num))
-        lk_grad_mi, 3*n_points_final_grid,  &
+    allocate(rk_grad_im(n_points_final_grid,3,mo_num))
-        rk_grad_im, 3*n_points_final_grid,  0.d0, &
+
-        tmp_mat, mo_num*mo_num)
+    do i=1,mo_num
-  
+      !$OMP PARALLEL DO PRIVATE (l,ipoint)
-    !$OMP PARALLEL DO PRIVATE(i,j,k,l)
+      do l=1,mo_num
-    do i = 1, mo_num
+        do ipoint=1, n_points_final_grid
+
+          rm_grad_ik(ipoint,1,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
+          rm_grad_ik(ipoint,2,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
+          rm_grad_ik(ipoint,3,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
+
+          rk_grad_im(ipoint,1,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
+          rk_grad_im(ipoint,2,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
+          rk_grad_im(ipoint,3,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
+
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lm_grad_ik, 3*n_points_final_grid,                         &
+          rm_grad_ik, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
       do k = 1, mo_num
         do j = 1, mo_num
           do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,j,i,k)
+            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,i,j,l)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,j,i,l)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,i,j,k)
           enddo
         enddo
       enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lm_grad_ik, 3*n_points_final_grid,                         &
+          rk_grad_im, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) = three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) - tmp_mat(l,k,j)
+            three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) = three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) - tmp_mat(k,j,l)
+            three_e_5_idx_exch23_bi_ort (m,i,j,k,l) = three_e_5_idx_exch23_bi_ort (m,i,j,k,l) - tmp_mat(k,l,j)
+            three_e_5_idx_exch13_bi_ort (m,l,j,i,k) = three_e_5_idx_exch13_bi_ort (m,l,j,i,k) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lk_grad_mi, 3*n_points_final_grid,                         &
+          rm_grad_ik, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l)
+            three_e_5_idx_cycle_2_bi_ort(m,l,i,k,j) = three_e_5_idx_cycle_2_bi_ort(m,l,i,k,j) - tmp_mat(l,j,k)
+            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k)
+            three_e_5_idx_exch13_bi_ort (m,l,i,k,j) = three_e_5_idx_exch13_bi_ort (m,l,i,k,j) - tmp_mat(k,j,l)
+          enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lk_grad_mi, 3*n_points_final_grid,                         &
+          rk_grad_im, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) = three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) - tmp_mat(l,j,k)
+            three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) = three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) - tmp_mat(k,l,j)
+            three_e_5_idx_exch23_bi_ort (m,i,j,k,l) = three_e_5_idx_exch23_bi_ort (m,i,j,k,l) - tmp_mat(k,j,l)
+            three_e_5_idx_exch13_bi_ort (m,l,j,i,k) = three_e_5_idx_exch13_bi_ort (m,l,j,i,k) - tmp_mat(l,k,j)
+          enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
     enddo
-    !$OMP END PARALLEL DO
+
-  
-    deallocate(lk_grad_mi)
     deallocate(rm_grad_ik)
     deallocate(rk_grad_im)
+    deallocate(lk_grad_mi)
+    deallocate(lm_grad_ik)
+
   enddo
 
   deallocate(tmp_mat)
 
+  deallocate(orb_mat)
+
   call wall_time(wall1)
   print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0
   call print_memory_usage()

src/scf_utils/fock_matrix.irp.f

@@ -5,6 +5,90 @@
    ! Fock matrix on the MO basis.
    ! For open shells, the ROHF Fock Matrix is ::
    !
+   !       |  Rcc  |  F^b  |  Fcv  |
+   !       |-----------------------|
+   !       |  F^b  |  Roo  |  F^a  |
+   !       |-----------------------|
+   !       |  Fcv  |  F^a  |  Rvv  |
+   !
+   ! C: Core, O: Open, V: Virtual 
+   ! 
+   ! Rcc = Acc Fcc^a + Bcc Fcc^b
+   ! Roo = Aoo Foo^a + Boo Foo^b
+   ! Rvv = Avv Fvv^a + Bvv Fvv^b
+   ! Fcv = (F^a + F^b)/2
+   ! 
+   ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO)
+   ! A,B: Coupling parameters
+   !
+   ! J. Chem. Phys. 133, 141102 (2010), https://doi.org/10.1063/1.3503173
+   ! Coupling parameters from J. Chem. Phys. 125, 204110 (2006); https://doi.org/10.1063/1.2393223.
+   !       cc   oo   vv
+   !  A  -0.5  0.5  1.5
+   !  B   1.5  0.5 -0.5
+   ! 
+   END_DOC
+   integer                        :: i,j,n
+   if (elec_alpha_num == elec_beta_num) then
+     Fock_matrix_mo = Fock_matrix_mo_alpha
+   else
+     ! Core
+     do j = 1, elec_beta_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) = - 0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 1.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+     enddo
+     ! Open
+     do j = elec_beta_num+1, elec_alpha_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j)
+       enddo
+     enddo
+     ! Virtual
+     do j = elec_alpha_num+1, mo_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) =   1.5d0 * fock_matrix_mo_alpha(i,j) &
+                               - 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+     enddo
+   endif
+
+   ! Old
+   ! BEGIN_DOC
+   ! Fock matrix on the MO basis.
+   ! For open shells, the ROHF Fock Matrix is ::
+   !
    !       |   F-K    |  F + K/2  |    F     |
    !       |---------------------------------|
    !       | F + K/2  |     F     |  F - K/2 |
@@ -16,64 +100,64 @@
    !
    ! K = Fb - Fa
    !
-   END_DOC
+   ! END_DOC
-   integer                        :: i,j,n
+   !integer                        :: i,j,n
-   if (elec_alpha_num == elec_beta_num) then
+   !if (elec_alpha_num == elec_beta_num) then
-     Fock_matrix_mo = Fock_matrix_mo_alpha
+   !  Fock_matrix_mo = Fock_matrix_mo_alpha
-   else
+   !else
 
-     do j=1,elec_beta_num
+   !  do j=1,elec_beta_num
-       ! F-K
+   !    ! F-K
-       do i=1,elec_beta_num !CC
+   !    do i=1,elec_beta_num !CC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-       ! F+K/2
+   !    ! F+K/2
-       do i=elec_beta_num+1,elec_alpha_num  !CA
+   !    do i=elec_beta_num+1,elec_alpha_num  !CA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-       ! F
+   !    ! F
-       do i=elec_alpha_num+1, mo_num !CV
+   !    do i=elec_alpha_num+1, mo_num !CV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
+   !    enddo
-     enddo
+   !  enddo
 
-     do j=elec_beta_num+1,elec_alpha_num
+   !  do j=elec_beta_num+1,elec_alpha_num
-       ! F+K/2
+   !    ! F+K/2
-       do i=1,elec_beta_num !AC
+   !    do i=1,elec_beta_num !AC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-       ! F
+   !    ! F
-       do i=elec_beta_num+1,elec_alpha_num !AA
+   !    do i=elec_beta_num+1,elec_alpha_num !AA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
+   !    enddo
-       ! F-K/2
+   !    ! F-K/2
-       do i=elec_alpha_num+1, mo_num !AV
+   !    do i=elec_alpha_num+1, mo_num !AV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-     enddo
+   !  enddo
 
-     do j=elec_alpha_num+1, mo_num
+   !  do j=elec_alpha_num+1, mo_num
-       ! F
+   !    ! F
-       do i=1,elec_beta_num !VC
+   !    do i=1,elec_beta_num !VC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
+   !    enddo
-       ! F-K/2
+   !    ! F-K/2
-       do i=elec_beta_num+1,elec_alpha_num !VA
+   !    do i=elec_beta_num+1,elec_alpha_num !VA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-       ! F+K
+   !    ! F+K
-       do i=elec_alpha_num+1,mo_num !VV
+   !    do i=elec_alpha_num+1,mo_num !VV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) &
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) &
-             + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !          + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
+   !    enddo
-     enddo
+   !  enddo
 
-   endif
+   !endif
 
    do i = 1, mo_num
      Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
@@ -115,8 +199,6 @@
 
 END_PROVIDER
 
-
-
 BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_num,mo_num) ]
    implicit none
    BEGIN_DOC