Removed all vvv in CCSD

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -92,7 +92,7 @@ subroutine run_ccsd_space_orb
       call compute_H_vv_chol(nO,nV,tau_x,H_vv)
       call compute_H_vo_chol(nO,nV,t1,H_vo)
 
-      call compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+      call compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
       call compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
     else
       call compute_H_oo(nO,nV,t1,t2,tau,H_oo)
@@ -588,8 +588,6 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
 
   deallocate(W_vvov,T_vvoo)
 
-
-
   ! r1(u,beta) = r1(u,beta) - (2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)) * tau(i,j,a,beta)
   ! r1(u,beta) = r1(u,beta) - W(i,j,a,u) * tau(i,j,a,beta)
   !do beta = 1, nV

src/ccsd/ccsd_space_orb_sub_chol.irp.f

@@ -186,14 +186,13 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
   deallocate(X_ovov)
 
   integer :: iblock, block_size, nVmax
-  double precision, allocatable :: W_vvov(:,:,:,:), T_vvoo(:,:,:,:)
+  double precision, allocatable :: W_vvov(:,:,:,:), W_vvov_tmp(:,:,:,:), T_vvoo(:,:,:,:)
-  block_size = 8
+  block_size = 16
-  allocate(W_vvov(nV,nV,nO,block_size), T_vvoo(nV,nV,nO,nO))
+  allocate(W_vvov(nV,nV,nO,block_size), W_vvov_tmp(nV,nO,nV,block_size), T_vvoo(nV,nV,nO,nO))
 
   !$omp parallel &
-  !$omp shared(nO,nV,cc_space_v_vvov,W_vvov,T_vvoo,tau) &
+  !$omp private(u,i,b,a) &
-  !$omp private(b,beta,i,a) &
+  !$omp default(shared)
-  !$omp default(none)
   !$omp do
   do u = 1, nO
     do i = 1, nO
@@ -204,26 +203,32 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
       enddo
     enddo
   enddo
-  !$omp end do nowait
+  !$omp end do
   !$omp end parallel
 
   do iblock = 1, nV, block_size
     nVmax = min(block_size,nV-iblock+1)
+
+    call dgemm('T','N', nV*nO, nV*nVmax, cholesky_ao_num, 1.d0, &
+      cc_space_v_vo_chol            , cholesky_ao_num, &
+      cc_space_v_vv_chol(1,1,iblock), cholesky_ao_num, &
+      0.d0, W_vvov_tmp, nV*nO)
+
     !$omp parallel &
-    !$omp shared(nO,nV,cc_space_v_vvov,W_vvov,T_vvoo,tau,nVmax,iblock) &
     !$omp private(b,i,a,beta) &
-    !$omp default(none)
+    !$omp default(shared)
-    !$omp do collapse(2)
+    do beta = 1,  nVmax
-    do beta = iblock, iblock + nVmax - 1
       do i = 1, nO
+        !$omp do
         do b = 1, nV
           do a = 1, nV
-            W_vvov(a,b,i,beta-iblock+1) = 2d0 * cc_space_v_vvov(a,b,i,beta) - cc_space_v_vvov(b,a,i,beta)
+            W_vvov(a,b,i,beta) = 2d0 * W_vvov_tmp(a,i,b,beta) - W_vvov_tmp(b,i,a,beta)
-          enddo
-        enddo
           enddo
         enddo
         !$omp end do nowait
+      enddo
+    enddo
+    !$omp barrier
     !$omp end parallel
 
     call dgemm('T','N',nO,nVmax,nO*nV*nV, &
@@ -234,6 +239,7 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
 
   deallocate(W_vvov,T_vvoo)
 
+
   double precision, allocatable :: W_oovo(:,:,:,:)
   allocate(W_oovo(nO,nO,nV,nO))
 
@@ -462,7 +468,7 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   call compute_J1_chol(nO,nV,t1,t2,cc_space_v_ovvo,cc_space_v_ovoo, &
        cc_space_v_vvoo,J1)
   call compute_K1_chol(nO,nV,t1,t2,cc_space_v_ovoo,cc_space_v_vvoo, &
-       cc_space_v_ovov,cc_space_v_vvov,K1)
+       cc_space_v_ovov,K1)
 
   ! Residual
   !r2 = 0d0
@@ -1346,7 +1352,7 @@ end
 
 ! K1
 
-subroutine compute_K1_chol(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
+subroutine compute_K1_chol(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,K1)
 
   implicit none
 
@@ -1354,7 +1360,7 @@ subroutine compute_K1_chol(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
   double precision, intent(in)  :: t1(nO, nV)
   double precision, intent(in)  :: t2(nO, nO, nV, nV)
   double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO), v_ovov(nO,nV,nO,nV)
-  double precision, intent(in)  :: v_vvov(nV,nV,nO,nV), v_ovoo(nO,nV,nO,nO)
+  double precision, intent(in)  :: v_ovoo(nO,nV,nO,nO)
   double precision, intent(out) :: K1(nO, nV, nO, nV)
 
   double precision, allocatable :: X(:,:,:,:), Y(:,:,:,:), Z(:,:,:,:)
@@ -1412,11 +1418,6 @@ subroutine compute_K1_chol(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
   double precision, allocatable :: K1tmp(:,:,:,:), t1v(:,:,:)
   allocate(K1tmp(nO,nO,nV,nV), t1v(cholesky_ao_num,nO,nO))
 
-!  call dgemm('N','N',nO,nV*nO*nV,nV, &
-!             1d0, t1    , size(t1,1), &
-!                  v_vvov, size(v_vvov,1), &
-!             1d0, K1    , size(K1,1))
-
   call dgemm('N','T', cholesky_ao_num*nO, nO, nV, 1.d0, &
     cc_space_v_ov_chol, cholesky_ao_num*nO, t1, nO, 0.d0, &
     t1v, cholesky_ao_num*nO)

src/utils_cc/mo_integrals_cc.irp.f

@@ -190,7 +190,40 @@ BEGIN_PROVIDER [double precision, cc_space_v_oooo, (cc_nOa, cc_nOa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_oooo,1)
+    n2 = size(cc_space_v_oooo,2)
+    n3 = size(cc_space_v_oooo,3)
+    n4 = size(cc_space_v_oooo,4)
+
+    double precision, allocatable :: buffer(:,:,:,:)
+    allocate(buffer(n1,n3,n2,n4))
+
+    call dgemm('T','N', n1*n3, n2*n4, cholesky_ao_num, 1.d0, &
+         cc_space_v_oo_chol, cholesky_ao_num, &
+         cc_space_v_oo_chol, cholesky_ao_num, 0.d0, buffer, n1*n3)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_oooo(i1,i2,i3,i4) = buffer(i1,i3,i2,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(buffer)
+
+  else
     call gen_v_space(cc_nOa,cc_nOa,cc_nOa,cc_nOa, cc_list_occ,cc_list_occ,cc_list_occ,cc_list_occ, cc_space_v_oooo)
+  endif
 
 END_PROVIDER
 
@@ -200,7 +233,40 @@ BEGIN_PROVIDER [double precision, cc_space_v_vooo, (cc_nVa, cc_nOa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_vooo,1)
+    n2 = size(cc_space_v_vooo,2)
+    n3 = size(cc_space_v_vooo,3)
+    n4 = size(cc_space_v_vooo,4)
+
+    double precision, allocatable :: buffer(:,:,:,:)
+    allocate(buffer(n1,n3,n2,n4))
+
+    call dgemm('T','N', n1*n3, n2*n4, cholesky_ao_num, 1.d0, &
+         cc_space_v_vo_chol, cholesky_ao_num, &
+         cc_space_v_oo_chol, cholesky_ao_num, 0.d0, buffer, n1*n3)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_vooo(i1,i2,i3,i4) = buffer(i1,i3,i2,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(buffer)
+
+  else
     call gen_v_space(cc_nVa,cc_nOa,cc_nOa,cc_nOa, cc_list_vir,cc_list_occ,cc_list_occ,cc_list_occ, cc_space_v_vooo)
+  endif
 
 END_PROVIDER
 
@@ -210,7 +276,32 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovoo, (cc_nOa, cc_nVa, cc_nOa, cc_n
 
   implicit none
 
+
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_ovoo,1)
+    n2 = size(cc_space_v_ovoo,2)
+    n3 = size(cc_space_v_ovoo,3)
+    n4 = size(cc_space_v_ovoo,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_ovoo(i1,i2,i3,i4) = cc_space_v_vooo(i2,i1,i4,i3)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nVa,cc_nOa,cc_nOa, cc_list_occ,cc_list_vir,cc_list_occ,cc_list_occ, cc_space_v_ovoo)
+  endif
 
 END_PROVIDER
 
@@ -220,7 +311,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovo, (cc_nOa, cc_nOa, cc_nVa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_oovo,1)
+    n2 = size(cc_space_v_oovo,2)
+    n3 = size(cc_space_v_oovo,3)
+    n4 = size(cc_space_v_oovo,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_oovo(i1,i2,i3,i4) = cc_space_v_vooo(i3,i2,i1,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nOa,cc_nVa,cc_nOa, cc_list_occ,cc_list_occ,cc_list_vir,cc_list_occ, cc_space_v_oovo)
+  endif
 
 END_PROVIDER
 
@@ -230,7 +345,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_ooov, (cc_nOa, cc_nOa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_oovo,1)
+    n2 = size(cc_space_v_oovo,2)
+    n3 = size(cc_space_v_oovo,3)
+    n4 = size(cc_space_v_oovo,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_ooov(i1,i2,i3,i4) = cc_space_v_ovoo(i1,i4,i3,i2)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nOa,cc_nOa,cc_nVa, cc_list_occ,cc_list_occ,cc_list_occ,cc_list_vir, cc_space_v_ooov)
+  endif
 
 END_PROVIDER
 
@@ -240,7 +379,40 @@ BEGIN_PROVIDER [double precision, cc_space_v_vvoo, (cc_nVa, cc_nVa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_vvoo,1)
+    n2 = size(cc_space_v_vvoo,2)
+    n3 = size(cc_space_v_vvoo,3)
+    n4 = size(cc_space_v_vvoo,4)
+
+    double precision, allocatable :: buffer(:,:,:,:)
+    allocate(buffer(n1,n3,n2,n4))
+
+    call dgemm('T','N', n1*n3, n2*n4, cholesky_ao_num, 1.d0, &
+         cc_space_v_vo_chol, cholesky_ao_num, &
+         cc_space_v_vo_chol, cholesky_ao_num, 0.d0, buffer, n1*n3)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_vvoo(i1,i2,i3,i4) = buffer(i1,i3,i2,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(buffer)
+
+  else
     call gen_v_space(cc_nVa,cc_nVa,cc_nOa,cc_nOa, cc_list_vir,cc_list_vir,cc_list_occ,cc_list_occ, cc_space_v_vvoo)
+  endif
 
 END_PROVIDER
 
@@ -250,7 +422,40 @@ BEGIN_PROVIDER [double precision, cc_space_v_vovo, (cc_nVa, cc_nOa, cc_nVa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_vovo,1)
+    n2 = size(cc_space_v_vovo,2)
+    n3 = size(cc_space_v_vovo,3)
+    n4 = size(cc_space_v_vovo,4)
+
+    double precision, allocatable :: buffer(:,:,:,:)
+    allocate(buffer(n1,n3,n2,n4))
+
+    call dgemm('T','N', n1*n3, n2*n4, cholesky_ao_num, 1.d0, &
+         cc_space_v_vv_chol, cholesky_ao_num, &
+         cc_space_v_oo_chol, cholesky_ao_num, 0.d0, buffer, n1*n3)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_vovo(i1,i2,i3,i4) = buffer(i1,i3,i2,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(buffer)
+
+  else
     call gen_v_space(cc_nVa,cc_nOa,cc_nVa,cc_nOa, cc_list_vir,cc_list_occ,cc_list_vir,cc_list_occ, cc_space_v_vovo)
+  endif
 
 END_PROVIDER
 
@@ -260,7 +465,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_voov, (cc_nVa, cc_nOa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_voov,1)
+    n2 = size(cc_space_v_voov,2)
+    n3 = size(cc_space_v_voov,3)
+    n4 = size(cc_space_v_voov,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_voov(i1,i2,i3,i4) = cc_space_v_vvoo(i1,i4,i3,i2)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nVa,cc_nOa,cc_nOa,cc_nVa, cc_list_vir,cc_list_occ,cc_list_occ,cc_list_vir, cc_space_v_voov)
+  endif
 
 END_PROVIDER
 
@@ -270,7 +499,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovvo, (cc_nOa, cc_nVa, cc_nVa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_ovvo,1)
+    n2 = size(cc_space_v_ovvo,2)
+    n3 = size(cc_space_v_ovvo,3)
+    n4 = size(cc_space_v_ovvo,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_ovvo(i1,i2,i3,i4) = cc_space_v_vvoo(i3,i2,i1,i4)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nVa,cc_nVa,cc_nOa, cc_list_occ,cc_list_vir,cc_list_vir,cc_list_occ, cc_space_v_ovvo)
+  endif
 
 END_PROVIDER
 
@@ -280,7 +533,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovov, (cc_nOa, cc_nVa, cc_nOa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_ovov,1)
+    n2 = size(cc_space_v_ovov,2)
+    n3 = size(cc_space_v_ovov,3)
+    n4 = size(cc_space_v_ovov,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_ovov(i1,i2,i3,i4) = cc_space_v_vovo(i2,i1,i4,i3)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nVa,cc_nOa,cc_nVa, cc_list_occ,cc_list_vir,cc_list_occ,cc_list_vir, cc_space_v_ovov)
+  endif
 
 END_PROVIDER
 
@@ -290,7 +567,31 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovv, (cc_nOa, cc_nOa, cc_nVa, cc_n
 
   implicit none
 
+  if (do_ao_cholesky) then
+
+    integer :: i1, i2, i3, i4
+    integer :: n1, n2, n3, n4
+    
+    n1 = size(cc_space_v_oovv,1)
+    n2 = size(cc_space_v_oovv,2)
+    n3 = size(cc_space_v_oovv,3)
+    n4 = size(cc_space_v_oovv,4)
+
+    !$OMP PARALLEL DO PRIVATE(i1,i2,i3,i4) COLLAPSE(2)
+    do i4 = 1, n4
+      do i3 = 1, n3
+        do i2 = 1, n2
+          do i1 = 1, n1
+            cc_space_v_oovv(i1,i2,i3,i4) = cc_space_v_vvoo(i3,i4,i1,i2)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+  else
     call gen_v_space(cc_nOa,cc_nOa,cc_nVa,cc_nVa, cc_list_occ,cc_list_occ,cc_list_vir,cc_list_vir, cc_space_v_oovv)
+  endif
 
 END_PROVIDER