diff --git a/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f b/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
index f22611b..0de043a 100644
--- a/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
+++ b/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
@@ -489,11 +489,15 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   !$omp end do
   !$omp end parallel
 
+  double precision, allocatable :: J1(:,:,:,:)
+  allocate(J1(nO,nV,nV,nO))
+  J1 = 0.d0
+
   call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,tau, &
       cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
       cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
-      cc_space_v_oovo, cc_space_v_ovvo, &
-      cc_space_f_vo, H_vv, g_occ, r2)
+      cc_space_v_oovo, cc_space_v_ovvo, cc_space_v_ovoo, &
+      cc_space_f_vo, H_vv, g_occ, J1, r2)
 
 
 !---
@@ -504,8 +508,6 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
 
   double precision, allocatable :: X_vovo(:,:,:,:), Y_oovo(:,:,:,:)
 
-  double precision, allocatable :: J1(:,:,:,:)
-  allocate(J1(nO,nV,nV,nO))
   call compute_J1_chol(nO,nV,t1,t2,cc_space_v_ovvo,cc_space_v_ovoo, &
        cc_space_v_vvoo,J1)
 
@@ -735,90 +737,15 @@ subroutine compute_J1_chol(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvoo,J1)
   double precision, intent(in)  :: t2(nO, nO, nV, nV)
   double precision, intent(in)  :: v_ovvo(nO,nV,nV,nO), v_ovoo(nO,nV,nO,nO)
   double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
-  double precision, intent(out) :: J1(nO, nV, nV, nO)
+  double precision, intent(inout) :: J1(nO, nV, nV, nO)
 
   integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
 
   double precision, allocatable :: X_ovoo(:,:,:,:), Y_ovov(:,:,:,:)
   allocate(X_ovoo(nO,nV,nO,nO),Y_ovov(nO,nV,nO,nV))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,J1,v_ovvo,v_ovoo,X_ovoo) &
-  !$omp private(i,j,a,u,beta) &
-  !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do beta = 1, nV
-      do a = 1, nV
-        do u = 1, nO
-          J1(u,a,beta,i) = v_ovvo(u,a,beta,i)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
-
-  !$omp do collapse(2)
-  do j = 1, nO
-    do i = 1, nO
-      do a = 1, nV
-        do u = 1, nO
-          X_ovoo(u,a,i,j) = v_ovoo(u,a,j,i)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-  call dgemm('N','N',nO*nV*nO,nV,nO, &
-            -1d0, X_ovoo, size(X_ovoo,1) * size(X_ovoo,2) * size(X_ovoo,3), &
-                  t1    , size(t1,1), &
-             0d0, Y_ovov, size(Y_ovov,1) * size(Y_ovov,2) * size(Y_ovov,3))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,J1,Y_ovov) &
-  !$omp private(i,beta,a,u) &
-  !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do beta = 1, nV
-      do a = 1, nV
-        do u = 1, nO
-          J1(u,a,beta,i) = J1(u,a,beta,i) + Y_ovov(u,a,i,beta)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
-  !$omp end parallel
   deallocate(X_ovoo)
 
   double precision, allocatable :: tmp_cc(:,:,:), J1_tmp(:,:,:,:)
-  allocate(tmp_cc(cholesky_mo_num,nV,nO), J1_tmp(nV,nO,nV,nO))
-
-  call dgemm('N','T', cholesky_mo_num*nV, nO, nV, 1.d0, &
-      cc_space_v_vv_chol, cholesky_mo_num*nV, &
-      t1, nO, &
-      0.d0, tmp_cc, cholesky_mo_num*nV)
-
-  call dgemm('T','N', nV*nO, nV*nO, cholesky_mo_num, 1.d0, &
-      tmp_cc, cholesky_mo_num, cc_space_v_vo_chol, cholesky_mo_num, &
-      0.d0, J1_tmp, nV*nO)
-
-  deallocate(tmp_cc)
-
-  do i=1,nO
-    do b=1,nV
-      do a=1,nV
-        do u=1,nO
-          J1(u,a,b,i) = J1(u,a,b,i) + J1_tmp(b,u,a,i)
-        enddo
-      enddo
-    enddo
-  enddo
-
-  deallocate(J1_tmp)
 
   !- cc_space_v_vvoo(a,b,i,j) * (0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)) &
   double precision, allocatable :: X_voov(:,:,:,:), Z_ovvo(:,:,:,:)
diff --git a/devel/ccsd_gpu/gpu.c b/devel/ccsd_gpu/gpu.c
index 041b84b..346a6cd 100644
--- a/devel/ccsd_gpu/gpu.c
+++ b/devel/ccsd_gpu/gpu.c
@@ -73,9 +73,11 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                 double* cc_space_v_vvoo,
                 double* cc_space_v_oovo,
                 double* cc_space_v_ovvo,
+                double* cc_space_v_ovoo,
                 double* cc_space_f_vo,
                 double* H_vv,
                 double* g_occ,
+                double* J1,
                 double* r2)
 {
 
@@ -132,6 +134,11 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
       cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(double), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
 
+      double* d_cc_space_v_ovvo;
+      lda = nO*nV;
+      cudaMalloc((void **)&d_cc_space_v_ovvo, nO*nV*nV*nO * sizeof(double));
+      cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_ovvo, lda, d_cc_space_v_ovvo, lda);
+
       lda = nO;
       cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
       cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
@@ -499,11 +506,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
           }
           cudaFree(d_cc_space_v_oovo);
 
-          double* d_cc_space_v_ovvo;
-          lda = nO*nV;
-          cudaMalloc((void **)&d_cc_space_v_ovvo, nO*nV*nV*nO * sizeof(double));
-          cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_ovvo, lda, d_cc_space_v_ovvo, lda);
-
           double* d_X_vovo;
           lda = nV*nO;
           cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(double));
@@ -518,7 +520,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
-          cudaFree(d_cc_space_v_ovvo);
 
           double* d_Y_oovo;
           lda = nO*nO;
@@ -558,6 +559,107 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
           }
           cudaFree(d_X_oovv);
         }
+
+
+        #pragma omp section
+        {
+
+          double* d_J1;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+
+          alpha = 1.0;
+          beta =  0.0;
+          A = d_cc_space_v_ovvo; lda = nO*nV;
+          B = d_cc_space_v_ovvo; ldb = nO*nV;
+          C = d_J1; ldc = nO*nV;
+          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nV*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+
+          double* d_cc_space_v_ovoo;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_cc_space_v_ovoo, nO*nV*nO*nO * sizeof(double));
+          cublasSetMatrix(lda, nO*nO, sizeof(double), cc_space_v_ovoo, lda, d_cc_space_v_ovoo, lda);
+
+          double* d_X_ovoo;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(double));
+          for (int j=0 ; j<nO ; ++j) {
+            for (int i=0 ; i<nO ; ++i) {
+              alpha = 0.0;
+              beta =  1.0;
+              A = &(d_X_ovoo[nO*nV*(i+nO*j)]); lda = nO;
+              B = &(d_cc_space_v_ovoo[nO*nV*(j+nO*i)]); ldb = nO;
+              C = &(d_X_ovoo[nO*nV*(i+nO*j)]); ldc = nO;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+
+          double* d_Y_ovov;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV*nO;  n=nV;  k=nO;
+          A=d_X_ovoo;  lda=nO*nV*nO;
+          B=d_t1;  ldb=nO;
+          C=d_Y_ovov;  ldc=nO*nV*nO;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          for (int j=0 ; j<nO ; ++j) {
+            for (int i=0 ; i<nV ; ++i) {
+              alpha = 1.0;
+              beta =  -1.0;
+              A = &(d_J1[nO*nV*(i+nV*j)]); lda = nO;
+              B = &(d_Y_ovov[nO*nV*(j+nO*i)]); ldb = nO;
+              C = &(d_J1[nO*nV*(i+nV*j)]); ldc = nO;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          double* d_tmp_cc;
+          lda = cholesky_mo_num;
+          cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(double));
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=cholesky_mo_num*nV;  n=nO;  k=nV;
+          A=d_cc_space_v_vv_chol;  lda=cholesky_mo_num*nV;
+          B=d_t1;  ldb=nO;
+          C=d_tmp_cc;  ldc=cholesky_mo_num*nV;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          double* d_J1_tmp;
+          lda = nV*nO;
+          cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(double));
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nV*nO;  n=nV*nO;  k=cholesky_mo_num;
+          A=d_tmp_cc;  lda=cholesky_mo_num;
+          B=d_cc_space_v_vo_chol;  ldb=cholesky_mo_num;
+          C=d_J1_tmp;  ldc=nV*nO;
+          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_tmp_cc);
+
+          for (int i=0 ; i<nO ; ++i) {
+              alpha = 1.0;
+              beta =  1.0;
+              A = &(d_J1[nO*nV*nV*i]); lda = nO*nV;
+              B = &(d_J1_tmp[nV*nO*nV*i]); ldb = nV;
+              C = &(d_J1[nO*nV*nV*i]); ldc = nO*nV;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO*nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          }
+          cudaFree(d_J1_tmp);
+
+lda = nO*nV;
+cublasGetMatrix(nO*nV, nV*nO, sizeof(double), d_J1, lda, J1, lda);
+          cudaFree(d_J1);
+
+        }
       }
 
 
@@ -662,169 +764,3 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
 }
 
 
-
-
-
-
-void compute_g_vir_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
-                double* t1,
-                double* t2,
-                double* tau,
-                double* cc_space_v_ov_chol,
-                double* cc_space_v_vo_chol,
-                double* cc_space_v_vv_chol,
-                double* cc_space_v_oooo,
-                double* cc_space_v_vooo,
-                double* cc_space_v_oovv,
-                double* cc_space_v_vvoo,
-                double* cc_space_f_vo,
-                double* H_vv,
-                double* r2)
-{
-      int m,n,k, lda, ldb, ldc;
-      double alpha, beta;
-      double* A;
-      double* B;
-      double* C;
-
-      int ithread = omp_get_thread_num();
-      int igpu = ithread ;
-
-      cudaSetDevice(igpu);
-      cublasHandle_t handle;
-
-      cublasCreate(&handle);
-
-      double* d_tau;
-      double* d_r2;
-      double* d_cc_space_v_vv_chol;
-      double* d_cc_space_v_vo_chol;
-      double* d_cc_space_v_ov_chol;
-      double* d_t1;
-      double* d_t2;
-      double* d_tmp_cc;
-
-
-
-      lda = nO * nO;
-      cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
-      cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
-
-      lda = nO*nO;
-      cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
-      cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
-      lda = nO;
-      cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
-      cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
-      lda = cholesky_mo_num * nV;
-      cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
-
-      lda = cholesky_mo_num * nO;
-      cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
-
- // ---
-           double* d_g_vir;
-           cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
-           cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_g_vir, nV);
-
-           double* d_cc_space_f_vo;
-           cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
-           cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
-
-           alpha = -1.0;
-           beta =  1.0;
-           m=nV ; n=nV; k=nO;
-           A = d_cc_space_f_vo  ; lda = nV;
-           B = d_t1 ; ldb = nO;
-           C = d_g_vir; ldc = nV;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_cc_space_f_vo);
-
-           double* d_tmp_k;
-           cudaMalloc((void**)&d_tmp_k, cholesky_mo_num*sizeof(double));
-           alpha = 1.0;
-           beta =  0.0;
-           m=cholesky_mo_num ; n=1; k=nO*nV;
-           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num;
-           B = d_t1 ; ldb = nO*nV;
-           C = d_tmp_k; ldc = cholesky_mo_num;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-           alpha = 2.0;
-           beta =  1.0;
-           m=nV*nV; n=1; k=cholesky_mo_num;
-           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num;
-           B = d_tmp_k ; ldb = cholesky_mo_num;
-           C = d_g_vir; ldc = nV*nV;
-           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_tmp_k);
-
-
-           double* d_tmp_vo;
-           cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(double));
-           alpha = 1.0;
-           beta =  0.0;
-           m=cholesky_mo_num*nV ; n=nO; k=nV;
-           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num*nV;
-           B = d_t1 ; ldb = nO;
-           C = d_tmp_vo; ldc = cholesky_mo_num*nV;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-           double* d_tmp_vo2;
-           cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(double));
-           for (int i=0 ; i<nO ; ++i) {
-             alpha = -1.0;
-             beta =  0.0;
-             A = &(d_tmp_vo[cholesky_mo_num*nV*i]); lda = cholesky_mo_num;
-             B = &(d_tmp_vo[cholesky_mo_num*nV*i]); ldb = cholesky_mo_num;
-             C = &(d_tmp_vo2[cholesky_mo_num*i]); ldc = cholesky_mo_num*nO;
-             cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
-           }
-           cudaFree(d_tmp_vo);
-
-           alpha = 1.0;
-           beta =  1.0;
-           m=nV ; n=nV; k=nO*cholesky_mo_num;
-           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num*nO;
-           B = d_tmp_vo2 ; ldb = cholesky_mo_num*nO;
-           C = d_g_vir; ldc = nV;
-           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_tmp_vo2);
-
-           double* d_Y_oovv;
-           cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(double));
-           alpha = 1.0;
-           beta =  0.0;
-           m=nO*nO*nV ; n=nV; k=nV;
-           A = d_t2; lda = nO*nO*nV;
-           B = d_g_vir; ldb = nV;
-           C = d_Y_oovv; ldc = nO*nO*nV;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_g_vir);
-
-
-           alpha = 1.0;
-           beta =  1.0;
-           A = d_r2; lda = nO*nO;
-           B = d_Y_oovv; ldb = nO*nO;
-           C = d_r2; ldc = nO*nO;
-           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
-           for (int j=0 ; j<nV ; ++j) {
-             for (int i=0 ; i<nV ; ++i) {
-                alpha = 1.0;
-                beta =  1.0;
-                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
-                B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
-                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
-                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
-
-             }
-           }
-           cudaFree(d_Y_oovv);
-
-           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, nO*nO, r2, nO*nO);
-}
-
-
diff --git a/devel/ccsd_gpu/gpu_module.f90 b/devel/ccsd_gpu/gpu_module.f90
index cb6f174..014f803 100644
--- a/devel/ccsd_gpu/gpu_module.f90
+++ b/devel/ccsd_gpu/gpu_module.f90
@@ -7,8 +7,8 @@ module gpu_module
     subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
           cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
           cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
-          cc_space_v_oovo, cc_space_v_ovvo, &
-          cc_space_f_vo, H_vv, g_occ, r2) bind(C)
+          cc_space_v_oovo, cc_space_v_ovvo, cc_space_v_ovoo, &
+          cc_space_f_vo, H_vv, g_occ, J1, r2) bind(C)
         import c_int, c_double
         integer(c_int), value :: nO, nV, cholesky_mo_num
         real(c_double), intent(in)  :: t1(nO,nV)
@@ -24,9 +24,11 @@ module gpu_module
         real(c_double), intent(in)  :: cc_space_v_vvoo(nV,nV,nO,nO)
         real(c_double), intent(in)  :: cc_space_v_oovo(nO,nO,nV,nO)
         real(c_double), intent(in)  :: cc_space_v_ovvo(nO,nV,nV,nO)
+        real(c_double), intent(in)  :: cc_space_v_ovoo(nO,nV,nO,nO)
         real(c_double), intent(in)  :: cc_space_f_vo(nV,nO)
         real(c_double), intent(in)  :: H_vv(nV,nV)
         real(c_double), intent(in)  :: g_occ(nO,nO)
+        real(c_double)              :: J1(nO,nV,nV,nO)
         real(c_double), intent(out) :: r2(nO,nO,nV,nV)
     end subroutine