Fast 4idx in CASSCF

src/casscf_cipsi/gradient.irp.f

@@ -163,14 +163,21 @@ real*8 function gradvec_it(i,t)
   tt=list_act(t)
   gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
   gradvec_it-=occnum(tt)*Fipq(ii,tt)
+  do y=1,n_act_orb ! active
+!   y3=y+n_core_inact_orb ! list_act(y)
+    do x=1,n_act_orb ! active
+!      x3=x+n_core_inact_orb ! list_act(x)
       do v=1,n_act_orb ! active
         vv=list_act(v)
-    do x=1,n_act_orb ! active 
-      x3=x+n_core_inact_orb ! list_act(x) 
-      do y=1,n_act_orb ! active
-        y3=y+n_core_inact_orb ! list_act(y)
         !                Gamma(2)  a a a a      1/r12      i  a  a  a
-        gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
+!        gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
+         integer :: ichol
+         double precision :: tmp
+         tmp = 0.d0
+         do ichol=1,cholesky_mo_num
+           tmp = tmp + cholesky_no_total_transp(ichol,vv,ii) * cholesky_no_total_transp(ichol,list_act(x),list_act(y))
+         enddo
+         gradvec_it = gradvec_it - 2.D0*P0tuvx_no(t,v,x,y)*tmp
       end do
     end do
   end do

src/casscf_cipsi/hessian.irp.f

@@ -17,7 +17,8 @@ real*8 function hessmat_itju(i,t,j,u)
   if (i.eq.j) then
     if (t.eq.u) then
       ! diagonal element
-      term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt))     &
+      term = occnum(tt)*Fipq(ii,ii) +      &
+           2.D0*(Fipq(tt,tt)+Fapq(tt,tt))     &
           -2.D0*(Fipq(ii,ii)+Fapq(ii,ii))
       term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i))
       term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt)             &

src/casscf_cipsi/natorb.irp.f

@@ -72,84 +72,27 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
   BEGIN_DOC
   ! 4-index transformation of 2part matrices
   END_DOC
-  integer                        :: i,j,k,l,p,q
-  real*8                         :: d(n_act_orb)
 
-  ! index per index
+  double precision, allocatable :: tmp(:,:,:,:)
-  ! first quarter
+  allocate(tmp(n_act_orb,n_act_orb,n_act_orb,n_act_orb))
-  P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:)
 
-  do j=1,n_act_orb
+  call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
-    do k=1,n_act_orb
+          P0tuvx, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
-      do l=1,n_act_orb
+          tmp, (n_act_orb*n_act_orb*n_act_orb))
-        do p=1,n_act_orb
+
-          d(p)=0.D0
+  call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
-        end do
+          tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
-        do p=1,n_act_orb
+          P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb))
-          do q=1,n_act_orb
+
-            d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
+  call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
-          end do
+          P0tuvx_no, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
-        end do
+          tmp, (n_act_orb*n_act_orb*n_act_orb))
-        do p=1,n_act_orb
+
-          P0tuvx_no(p,j,k,l)=d(p)
+  call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
-        end do
+          tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
-      end do
+          P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb))
-    end do
+
-  end do
+  deallocate(tmp)
-  ! 2nd quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,p,k,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 3rd quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,k,p,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          P0tuvx_no(j,k,l,p)=d(p)
-        end do
-      end do
-    end do
-  end do
 
 END_PROVIDER
 
@@ -160,6 +103,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
   BEGIN_DOC
   ! Transformed one-e integrals
   END_DOC
+
   integer :: i,j, p, q
   real*8 :: d(n_act_orb)
   one_ints_no(:,:)=mo_one_e_integrals(:,:)
@@ -168,10 +112,8 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
   do j=1,mo_num
     do p=1,n_act_orb
       d(p)=0.D0
-    end do
-    do p=1,n_act_orb
       do q=1,n_act_orb
-        d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
+        d(p) = d(p) + one_ints_no(list_act(q),j)*natorbsCI(q,p)
       end do
     end do
     do p=1,n_act_orb
@@ -183,8 +125,6 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
   do j=1,mo_num
     do p=1,n_act_orb
       d(p)=0.D0
-    end do
-    do p=1,n_act_orb
       do q=1,n_act_orb
         d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
       end do

src/mo_two_e_ints/map_integrals.irp.f

@@ -86,10 +86,8 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_s
 
    call set_multiple_levels_omp(.False.)
 
-
    !$OMP PARALLEL DO PRIVATE(k,l,ii) SCHEDULE(dynamic)
    do l=mo_integrals_cache_min,mo_integrals_cache_max
-     print *, l
      do k=mo_integrals_cache_min,mo_integrals_cache_max
          ii = int(l-mo_integrals_cache_min,8)
          ii = ior( shiftl(ii,mo_integrals_cache_shift), int(k-mo_integrals_cache_min,8))