Swap loop indices in CC

devel/cc/form_cW_nc.irp.f

@@ -17,6 +17,7 @@ subroutine form_cW_nc(nO,nV,t1,t2,tau,cWoooo,cWovvo,cWvvvv)
   integer                       :: i,j,m,n
   integer                       :: a,b,e,f
   double precision,external     :: Kronecker_Delta
+  double precision :: x
 
 ! Output variables
 
@@ -24,25 +25,33 @@ subroutine form_cW_nc(nO,nV,t1,t2,tau,cWoooo,cWovvo,cWvvvv)
   double precision,intent(out)  :: cWovvo(nO,nV,nV,nO)
   double precision,intent(out)  :: cWvvvv(nV,nV,nV,nV)
 
+
 ! OOOO block  
 
   cWoooo(:,:,:,:) = OOOO(:,:,:,:)
 
-  do m=1,nO
+  do e=1,nV
-    do n=1,nO
+    do j=1,nO
       do i=1,nO
-        do j=1,nO
+        do n=1,nO
-
+          do m=1,nO
-          do e=1,nV
             cWoooo(m,n,i,j) = cWoooo(m,n,i,j) + t1(j,e)*OOOV(m,n,i,e) - t1(i,e)*OOOV(m,n,j,e)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do e=1,nV
+  do f=1,nV
-            do f=1,nV
+    do e=1,nV
-            cWoooo(m,n,i,j) = cWoooo(m,n,i,j) + 0.25d0*tau(i,j,e,f)*OOVV(m,n,e,f)
+      do j=1,nO
+        do i=1,nO
+          x = 0.25d0*tau(i,j,e,f)
+          do n=1,nO
+            do m=1,nO
+              cWoooo(m,n,i,j) = cWoooo(m,n,i,j) + x*OOVV(m,n,e,f)
             end do
           end do
-
         end do
       end do
     end do
@@ -52,26 +61,40 @@ subroutine form_cW_nc(nO,nV,t1,t2,tau,cWoooo,cWovvo,cWvvvv)
 
   cWovvo(:,:,:,:) = OVVO(:,:,:,:)
 
-  do m=1,nO
+  do f=1,nV
-    do b=1,nV
+    do j=1,nO
       do e=1,nV
-        do j=1,nO
+        do b=1,nV
-         
+          do m=1,nO
-          do f=1,nV
             cWovvo(m,b,e,j) = cWovvo(m,b,e,j) + t1(j,f)*OVVV(m,b,e,f)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do n=1,nO
+  do j=1,nO
+    do e=1,nV
+      do b=1,nV
+        do n=1,nO
+          do m=1,nO
             cWovvo(m,b,e,j) = cWovvo(m,b,e,j) - t1(n,b)*OOVO(m,n,e,j)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do n=1,nO
+  do j=1,nO
-            do f=1,nV
+    do f=1,nV
-              cWovvo(m,b,e,j) = cWovvo(m,b,e,j) &
+      do b=1,nV
-                              - ( 0.5d0*t2(j,n,f,b) + t1(j,f)*t1(n,b) )*OOVV(m,n,e,f)
+        do n=1,nO
+          x = 0.5d0*t2(j,n,f,b)  + t1(j,f)*t1(n,b)
+          do e=1,nV
+            do m=1,nO
+              cWovvo(m,b,e,j) = cWovvo(m,b,e,j) - x *OOVV(m,n,e,f)
             end do
           end do
-
         end do
       end do
     end do
@@ -81,21 +104,29 @@ subroutine form_cW_nc(nO,nV,t1,t2,tau,cWoooo,cWovvo,cWvvvv)
 
   cWvvvv(:,:,:,:) = VVVV(:,:,:,:)
 
-  do a=1,nV
+  do f=1,nV
-    do b=1,nV
+    do e=1,nV
-      do e=1,nV
+      do b=1,nV
-        do f=1,nV
+        do a=1,nV
-
           do m=1,nO
             cWvvvv(a,b,e,f) = cWvvvv(a,b,e,f) - t1(m,b)*VOVV(a,m,e,f) + t1(m,a)*VOVV(b,m,e,f)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do m=1,nO
+  do f=1,nV
-            do n=1,nO
+    do e=1,nV
-              cWvvvv(a,b,e,f) = cWvvvv(a,b,e,f) + 0.25d0*tau(m,n,a,b)*OOVV(m,n,e,f)
+      do b=1,nV
+        do a=1,nV
+          x = 0.d0
+          do n=1,nO
+            do m=1,nO
+              x = x + tau(m,n,a,b)*OOVV(m,n,e,f)
             end do
           end do
-
+          cWvvvv(a,b,e,f) = cWvvvv(a,b,e,f) + 0.25d0*x
         end do
       end do
     end do

devel/cc/form_r2_nc.irp.f

@@ -6,90 +6,343 @@ subroutine form_r2_nc(nO,nV,t1,t2,tau,cFoo,cFov,cFvv,cWoooo,cWvvvv,cWovvo,r2)
 
 ! Input variables
 
-  integer,intent(in)            :: nO,nV
+  integer,intent(in)             :: nO,nV
   
-  double precision,intent(in)   :: cFoo(nO,nO)
+  double precision,intent(in)    :: cFoo(nO,nO)
-  double precision,intent(in)   :: cFov(nO,nV)
+  double precision,intent(in)    :: cFov(nO,nV)
-  double precision,intent(in)   :: cFvv(nV,nV)
+  double precision,intent(in)    :: cFvv(nV,nV)
   
-  double precision,intent(in)   :: cWoooo(nO,nO,nO,nO)
+  double precision,intent(in)    :: cWoooo(nO,nO,nO,nO)
-  double precision,intent(in)   :: cWvvvv(nV,nV,nV,nV)
+  double precision,intent(in)    :: cWvvvv(nV,nV,nV,nV)
-  double precision,intent(in)   :: cWovvo(nO,nV,nV,nO)
+  double precision,intent(in)    :: cWovvo(nO,nV,nV,nO)
   
-  double precision,intent(in)   :: t1(nO,nV)
+  double precision,intent(in)    :: t1(nO,nV)
-  double precision,intent(in)   :: t2(nO,nO,nV,nV)
+  double precision,intent(in)    :: t2(nO,nO,nV,nV)
-  double precision,intent(in)   :: tau(nO,nO,nV,nV)
+  double precision,intent(in)    :: tau(nO,nO,nV,nV)
+  double precision               :: x
 
 ! Local variables
 
-  integer                       :: i,j,m,n
+  integer                        :: i,j,m,n
-  integer                       :: a,b,e,f
+  integer                        :: a,b,e,f
 
 ! Output variables
 
-  double precision,intent(out)  :: r2(nO,nO,nV,nV)
+  double precision,intent(out)   :: r2(nO,nO,nV,nV)
 
   r2(:,:,:,:) = OOVV(:,:,:,:)
 
-  do i=1,nO
+  do e=1,nV
-    do j=1,nO
+    do b=1,nV
       do a=1,nV
-        do b=1,nV
+        do j=1,nO
-
+          do i=1,nO
-          do e=1,nV
             r2(i,j,a,b) = r2(i,j,a,b) + t2(i,j,a,e)*cFvv(b,e)
+          end do
+        end do
+      end do
+    end do
+  end do
+
+
+  do e=1,nV
+    do b=1,nV
+      do a=1,nV
+        do j=1,nO
+          do i=1,nO
             r2(i,j,a,b) = r2(i,j,a,b) - t2(i,j,b,e)*cFvv(a,e)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do m=1,nO
+
-            do e=1,nV
+ do b=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b) - 0.5d0*t2(i,j,a,e)*t1(m,b)*cFov(m,e)
+   do e=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(i,j,b,e)*t1(m,a)*cFov(m,e)
+     do a=1,nV
-            end do
+       do j=1,nO
+         do i=1,nO
+           x = 0.5d0*t2(i,j,a,e)
+           if (x /= 0.d0) then
+              do m=1,nO
+                r2(i,j,a,b) = r2(i,j,a,b) - x*t1(m,b)*cFov(m,e)
+              end do
+            endif
           end do
+        end do
+      end do
+    end do
+  end do
 
+
+  do e=1,nV
+    do b=1,nV
+      do a=1,nV
+        do j=1,nO
+          do i=1,nO
+            x = 0.5d0*t2(i,j,b,e)
+            if (x /= 0.d0) then
+              do m=1,nO
+                r2(i,j,a,b) = r2(i,j,a,b) + x*t1(m,a)*cFov(m,e)
+              end do
+            endif
+          end do
+        end do
+      end do
+    end do
+  end do
+  
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do m=1,nO
+          x = cFoo(m,j)
+          do i=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,a,b)*x
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
           do m=1,nO
-            r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,a,b)*cFoo(m,j)
             r2(i,j,a,b) = r2(i,j,a,b) + t2(j,m,a,b)*cFoo(m,i)
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do m=1,nO
+
-            do e=1,nV
+  do b=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b) - 0.5d0*t2(i,m,a,b)*t1(j,e)*cFov(m,e)
+    do a=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(j,m,a,b)*t1(i,e)*cFov(m,e)
+      do e=1,nV
+        do m=1,nO
+          x = 0.5d0*cFov(m,e)
+          if (x == 0.d0) cycle
+          do j=1,nO
+            do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - x*t2(i,m,a,b)*t1(j,e)
             end do
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do m=1,nO
+
-            do n=1,nO
+  do b=1,nV
+    do a=1,nV
+      do e=1,nV
+        do m=1,nO
+          x = 0.5d0*cFov(m,e)
+          do j=1,nO
+            do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + x*t2(j,m,a,b)*t1(i,e)
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
+          do n=1,nO
+            do m=1,nO
               r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*tau(m,n,a,b)*cWoooo(m,n,i,j)
             end do
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do e=1,nV
+
-            do f=1,nV
+  do f=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*tau(i,j,e,f)*cWvvvv(a,b,e,f)
+    do e=1,nV
+      do b=1,nV
+        do a=1,nV
+          x = 0.5d0*cWvvvv(a,b,e,f) 
+          do j=1,nO
+            do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + x*tau(i,j,e,f)
             end do
           end do
+        end do
+      end do
+    end do
+  end do
 
-          do m=1,nO
-            do e=1,nV
-              r2(i,j,a,b) = r2(i,j,a,b)                                                 & 
-                          + t2(i,m,a,e)*cWovvo(m,b,e,j) - t1(i,e)*t1(m,a)*OVVO(m,b,e,j) &
-                          - t2(j,m,a,e)*cWovvo(m,b,e,i) + t1(j,e)*t1(m,a)*OVVO(m,b,e,i) &
-                          - t2(i,m,b,e)*cWovvo(m,a,e,j) + t1(i,e)*t1(m,b)*OVVO(m,a,e,j) &
-                          + t2(j,m,b,e)*cWovvo(m,a,e,i) - t1(j,e)*t1(m,b)*OVVO(m,a,e,i)
-            end do
-          end do
 
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
         do e=1,nV
-          r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*VVVO(a,b,e,j) - t1(j,e)*VVVO(a,b,e,i)
+          do m=1,nO
+            do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + t2(i,m,a,e)*cWovvo(m,b,e,j) 
+            end do
+          end do
         end do
+      end do
+    end do
+  end do
 
-        do m=1,nO
+  do b=1,nV
-          r2(i,j,a,b) = r2(i,j,a,b) - t1(m,a)*OVOO(m,b,i,j) + t1(m,b)*OVOO(m,a,i,j)
+    do a=1,nV
+      do j=1,nO
+        do e=1,nV
+          do m=1,nO
+            x = t1(m,a)*OVVO(m,b,e,j)
+            do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - t1(i,e)*x
+            end do
+          end do
         end do
+      end do
+    end do
+  end do
 
+  do b=1,nV
+    do e=1,nV
+      do a=1,nV
+        do j=1,nO
+          do i=1,nO
+            do m=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - t2(j,m,a,e)*cWovvo(m,b,e,i) 
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do b=1,nV
+    do a=1,nV
+      do e=1,nV
+        do j=1,nO
+          do i=1,nO
+            do m=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + t1(j,e)*t1(m,a)*OVVO(m,b,e,i)
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+ do e=1,nV
+   do b=1,nV
+     do a=1,nV
+       do j=1,nO
+         do m=1,nO
+           do i=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,b,e)*cWovvo(m,a,e,j) 
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do e=1,nV
+          do i=1,nO
+            do m=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*t1(m,b)*OVVO(m,a,e,j) 
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do e=1,nV
+    do b=1,nV
+      do a=1,nV
+        do j=1,nO
+          do i=1,nO
+            do m=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) + t2(j,m,b,e)*cWovvo(m,a,e,i) 
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do b=1,nV
+    do e=1,nV
+      do a=1,nV
+        do j=1,nO
+          do i=1,nO
+            do m=1,nO
+              r2(i,j,a,b) = r2(i,j,a,b) - t1(j,e)*t1(m,b)*OVVO(m,a,e,i)
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do e=1,nV
+          do i=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*VVVO(a,b,e,j) 
+          end do
+        end do
+      end do
+    end do
+  end do
+
+
+  do b=1,nV
+    do a=1,nV
+      do e=1,nV
+        do i=1,nO
+          do j=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) - t1(j,e)*VVVO(a,b,e,i)
+          end do
+        end do
+      end do
+    end do
+  end do
+
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
+          do m=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) - t1(m,a)*OVOO(m,b,i,j) 
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  do b=1,nV
+    do a=1,nV
+      do j=1,nO
+        do i=1,nO
+          do m=1,nO
+            r2(i,j,a,b) = r2(i,j,a,b) + t1(m,b)*OVOO(m,a,i,j)
+          end do
         end do
       end do
     end do