More cleaning

src/casscf/bielec.irp.f

@@ -1,104 +1,151 @@
-! -*- F90 -*- 
  BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
-&BEGIN_PROVIDER[real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
   BEGIN_DOC
   ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
-! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
   ! indices are unshifted orbital numbers
-! all integrals are read from files
   END_DOC
   implicit none
-          integer :: i,j,p,q,indx,kk
+  integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
-          real*8 :: hhh
+  double precision, allocatable  :: integrals_array(:,:)
-          logical :: lread
+  real*8                         :: mo_two_e_integral
   
-          do i=1,n_core_orb+n_act_orb
+  allocate(integrals_array(mo_num,mo_num))
-           do j=1,n_core_orb+n_act_orb
+  
+  bielec_PQxx = 0.d0
+  
+  do i=1,n_core_orb
+    ii=list_core(i)
+    do j=i,n_core_orb
+      jj=list_core(j)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
       do p=1,mo_num
         do q=1,mo_num
-              bielec_PQxx(p,q,i,j)=0.D0
+          bielec_PQxx(p,q,i,j)=integrals_array(p,q)
-              bielec_PxxQ(p,i,j,q)=0.D0
+          bielec_PQxx(p,q,j,i)=integrals_array(p,q)
+        end do
+      end do
+    end do
+    do j=1,n_act_orb
+      jj=list_act(j)
+      j3=j+n_core_orb
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do p=1,mo_num
+        do q=1,mo_num
+          bielec_PQxx(p,q,i,j3)=integrals_array(p,q)
+          bielec_PQxx(p,q,j3,i)=integrals_array(p,q)
         end do
       end do
     end do
   end do
 
-          open(unit=12,form='formatted',status='old',file='bielec_PQxx.tmp')
-          lread=.true.
-          indx=0
-          do while (lread)
-           read(12,*,iostat=kk) p,q,i,j,hhh
-           if (kk.ne.0) then
-            lread=.false.
-           else
-! stored with p.le.q, i.le.j
-            bielec_PQxx(p,q,i,j)=hhh
-            bielec_PQxx(q,p,i,j)=hhh
-            bielec_PQxx(q,p,j,i)=hhh
-            bielec_PQxx(p,q,j,i)=hhh
-            indx+=1
-           end if
-          end do
-          close(12)
-          write(6,*) ' read ',indx,' integrals PQxx into core '
 
-          open(unit=12,form='formatted',status='old',file='bielec_PxxQ.tmp')
+  ! (ij|pq)
-          lread=.true.
+  do i=1,n_act_orb
-          indx=0
+    ii=list_act(i)
-          do while (lread)
+    i3=i+n_core_orb
-           read(12,*,iostat=kk) p,i,j,q,hhh
+    do j=i,n_act_orb
-           if (kk.ne.0) then
+      jj=list_act(j)
-            lread=.false.
+      j3=j+n_core_orb
-           else
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
-! stored with (ip).le.(jq)
+      do p=1,mo_num
-            bielec_PxxQ(p,i,j,q)=hhh
+        do q=1,mo_num
-            bielec_PxxQ(q,j,i,p)=hhh
+          bielec_PQxx(p,q,i3,j3)=integrals_array(p,q)
-            indx+=1
+          bielec_PQxx(p,q,j3,i3)=integrals_array(p,q)
-           end if
         end do
-          write(6,*) ' read ',indx,' integrals PxxQ into core '
+      end do
-          close(12)
+    end do
-          write(6,*) ' provided integrals (PQ|xx) and (Px|xQ) '
+  end do
+
+  write(6,*) ' provided integrals (PQ|xx) '
 END_PROVIDER
 
+
+
+BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
+  BEGIN_DOC
+  ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
+  ! indices are unshifted orbital numbers
+  END_DOC
+  implicit none
+  integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
+  double precision, allocatable  :: integrals_array(:,:)
+  real*8                         :: mo_two_e_integral
+  
+  allocate(integrals_array(mo_num,mo_num))
+  
+  bielec_PxxQ = 0.d0
+  
+  do i=1,n_core_orb
+    ii=list_core(i)
+    do j=i,n_core_orb
+      jj=list_core(j)
+      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do p=1,mo_num
+        do q=1,mo_num
+          bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
+          bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
+        end do
+      end do
+    end do
+    do j=1,n_act_orb
+      jj=list_act(j)
+      j3=j+n_core_orb
+      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do p=1,mo_num
+        do q=1,mo_num
+          bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
+          bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
+        end do
+      end do
+    end do
+  end do
+
+
+  ! (ip|qj)
+  do i=1,n_act_orb
+    ii=list_act(i)
+    i3=i+n_core_orb
+    do j=i,n_act_orb
+      jj=list_act(j)
+      j3=j+n_core_orb
+      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do p=1,mo_num
+        do q=1,mo_num
+          bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
+          bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
+        end do
+      end do
+    end do
+  end do
+  write(6,*) ' provided integrals (Px|xQ) '
+END_PROVIDER
+
+
 BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   BEGIN_DOC
   ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
   ! index p runs over the whole basis, t,u,v only over the active orbitals
-! inegrals read from file
   END_DOC
   implicit none
-          integer :: i,j,k,p,t,u,v,kk,indx
+  integer                        :: i,j,k,p,t,u,v
-          real*8 :: hhh
+  double precision, allocatable  :: integrals_array(:)
-          logical :: lread
+  real*8                         :: mo_two_e_integral
   
-          write(6,*) ' reading integrals bielecCI '
+  allocate(integrals_array(mo_num))
   
   do i=1,n_act_orb
+    t=list_act(i)
     do j=1,n_act_orb
+      u=list_act(j)
       do k=1,n_act_orb
+        v=list_act(k)
+        ! (tu|vp)
+        call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map)
         do p=1,mo_num
-              bielecCI(i,k,j,p)=0.D0
+          bielecCI(i,k,j,p)=integrals_array(p)
         end do
       end do
     end do
   end do
-
-          open(unit=12,form='formatted',status='old',file='bielecCI.tmp')
-          lread=.true.
-          indx=0
-          do while (lread)
-           read(12,*,iostat=kk) i,j,k,p,hhh
-           if (kk.ne.0) then
-            lread=.false.
-           else
-            bielecCI(i,j,k,p)=hhh
-            bielecCI(j,i,k,p)=hhh
-            indx+=1
-           end if
-          end do
-          write(6,*) ' read ',indx,' integrals (aa|aP) into core '
-          close(12)
   write(6,*) ' provided integrals (tu|xP) '
 END_PROVIDER
 

src/casscf/bielec_create.irp.f

@@ -1,118 +0,0 @@
-! -*- F90 -*- 
- BEGIN_PROVIDER[real*8, bielec_PQxxtmp, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
-&BEGIN_PROVIDER[real*8, bielec_PxxQtmp, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
-&BEGIN_PROVIDER[integer, num_PQxx_written]
-&BEGIN_PROVIDER[integer, num_PxxQ_written]
-BEGIN_DOC
-! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
-! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
-! indices are unshifted orbital numbers
-END_DOC
-        implicit none
-          integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
-          double precision, allocatable :: integrals_array1(:,:)
-          double precision, allocatable :: integrals_array2(:,:)
-          real*8 :: mo_two_e_integral
-
-          allocate(integrals_array1(mo_num,mo_num))
-          allocate(integrals_array2(mo_num,mo_num))
-
-          do i=1,n_core_orb+n_act_orb
-           do j=1,n_core_orb+n_act_orb
-            do p=1,mo_num
-             do q=1,mo_num
-              bielec_PQxxtmp(p,q,i,j)=0.D0
-              bielec_PxxQtmp(p,i,j,q)=0.D0
-             end do
-            end do
-           end do
-          end do
-
-          do i=1,n_core_orb
-           ii=list_core(i)
-           do j=i,n_core_orb
-            jj=list_core(j)
-! (ij|pq)
-            call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
-! (ip|qj)
-            call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array2,mo_integrals_map)
-            do p=1,mo_num
-             do q=1,mo_num
-              bielec_PQxxtmp(p,q,i,j)=integrals_array1(p,q)
-              bielec_PQxxtmp(p,q,j,i)=integrals_array1(p,q)
-              bielec_PxxQtmp(p,i,j,q)=integrals_array2(p,q)
-              bielec_PxxQtmp(p,j,i,q)=integrals_array2(q,p)
-             end do
-            end do
-           end do
-           do j=1,n_act_orb
-            jj=list_act(j)
-            j3=j+n_core_orb
-! (ij|pq)
-            call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
-! (ip|qj)
-            call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array2,mo_integrals_map)
-            do p=1,mo_num
-             do q=1,mo_num
-              bielec_PQxxtmp(p,q,i,j3)=integrals_array1(p,q)
-              bielec_PQxxtmp(p,q,j3,i)=integrals_array1(p,q)
-              bielec_PxxQtmp(p,i,j3,q)=integrals_array2(p,q)
-              bielec_PxxQtmp(p,j3,i,q)=integrals_array2(q,p)
-             end do
-            end do
-           end do
-          end do
-          do i=1,n_act_orb
-           ii=list_act(i)
-           i3=i+n_core_orb
-           do j=i,n_act_orb
-            jj=list_act(j)
-            j3=j+n_core_orb
-! (ij|pq)
-            call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
-! (ip|qj)
-            call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array2,mo_integrals_map)
-            do p=1,mo_num
-             do q=1,mo_num
-              bielec_PQxxtmp(p,q,i3,j3)=integrals_array1(p,q)
-              bielec_PQxxtmp(p,q,j3,i3)=integrals_array1(p,q)
-              bielec_PxxQtmp(p,i3,j3,q)=integrals_array2(p,q)
-              bielec_PxxQtmp(p,j3,i3,q)=integrals_array2(q,p)
-             end do
-            end do
-           end do
-          end do
-          write(6,*) ' provided integrals (PQ|xx) '
-          write(6,*) ' provided integrals (Px|xQ) '
-!!$          write(6,*) ' 1 1 1 2 = ',bielec_PQxxtmp(2,2,2,3),bielec_PxxQtmp(2,2,2,3)
-END_PROVIDER
-
-BEGIN_PROVIDER[real*8, bielecCItmp, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
-BEGIN_DOC
-! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
-! index p runs over the whole basis, t,u,v only over the active orbitals
-END_DOC
-        implicit none
-          integer :: i,j,k,p,t,u,v
-          double precision, allocatable :: integrals_array1(:)
-          real*8 :: mo_two_e_integral
-
-          allocate(integrals_array1(mo_num))
- 
-          do i=1,n_act_orb
-           t=list_act(i)
-           do j=1,n_act_orb
-            u=list_act(j)
-            do k=1,n_act_orb
-             v=list_act(k)
-! (tu|vp)
-             call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array1,mo_integrals_map)
-             do p=1,mo_num
-              bielecCItmp(i,k,j,p)=integrals_array1(p)
-             end do
-            end do
-           end do
-          end do
-          write(6,*) ' provided integrals (tu|xP) '
-END_PROVIDER
-

src/casscf/bielec_natorb.irp.f

@@ -0,0 +1,273 @@
+ BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
+  BEGIN_DOC
+  ! integral (pq|xx) in the basis of natural MOs
+  ! indices are unshifted orbital numbers
+  END_DOC
+  implicit none
+  integer                        :: i,j,k,l,t,u,p,q,pp
+  real*8                         :: d(n_act_orb)
+
+  bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:)
+
+  do j=1,mo_num
+    do k=1,n_core_orb+n_act_orb
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PQxx_no(list_act(p),j,k,l)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 2nd quarter
+  do j=1,mo_num
+    do k=1,n_core_orb+n_act_orb
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PQxx_no(j,list_act(p),k,l)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 3rd quarter
+  do j=1,mo_num
+    do k=1,mo_num
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PQxx_no(j,k,n_core_orb+q,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PQxx_no(j,k,n_core_orb+p,l)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 4th quarter
+  do j=1,mo_num
+    do k=1,mo_num
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PQxx_no(j,k,l,n_core_orb+q)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PQxx_no(j,k,l,n_core_orb+p)=d(p)
+        end do
+      end do
+    end do
+  end do
+  write(6,*) ' transformed PQxx'
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
+  BEGIN_DOC
+  ! integral (px|xq) in the basis of natural MOs
+  ! indices are unshifted orbital numbers
+  END_DOC
+  implicit none
+  integer                        :: i,j,k,l,t,u,p,q,pp
+  real*8                         :: d(n_act_orb)
+
+  bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:)
+
+  do j=1,mo_num
+    do k=1,n_core_orb+n_act_orb
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0             
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PxxQ_no(list_act(p),k,l,j)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 2nd quarter
+  do j=1,mo_num
+    do k=1,n_core_orb+n_act_orb
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PxxQ_no(j,k,l,list_act(p))=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 3rd quarter
+  do j=1,mo_num
+    do k=1,mo_num
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PxxQ_no(j,n_core_orb+q,l,k)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PxxQ_no(j,n_core_orb+p,l,k)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 4th quarter
+  do j=1,mo_num
+    do k=1,mo_num
+      do l=1,n_core_orb+n_act_orb
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielec_PxxQ_no(j,l,n_core_orb+q,k)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielec_PxxQ_no(j,l,n_core_orb+p,k)=d(p)
+        end do
+      end do
+    end do
+  end do
+  write(6,*) ' transformed PxxQ '
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
+  BEGIN_DOC
+  ! integrals (tu|vp) in the basis of natural MOs
+  ! index p runs over the whole basis, t,u,v only over the active orbitals
+  END_DOC
+  implicit none
+  integer                        :: i,j,k,l,t,u,p,q,pp
+  real*8                         :: d(n_act_orb)
+  
+  bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:)
+
+  do j=1,n_act_orb
+    do k=1,n_act_orb
+      do l=1,mo_num
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielecCI_no(p,j,k,l)=d(p)
+        end do
+      end do
+    end do
+  end do
+  ! 2nd quarter
+  do j=1,n_act_orb
+    do k=1,n_act_orb
+      do l=1,mo_num
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielecCI_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,mo_num
+        do p=1,n_act_orb
+          d(p)=0.D0
+        end do
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielecCI_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
+          pp=n_act_orb-p+1
+          do q=1,n_act_orb
+            d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
+          end do
+        end do
+        do p=1,n_act_orb
+          bielecCI_no(j,k,l,list_act(p))=d(p)
+        end do
+      end do 
+    end do  
+  end do  
+  write(6,*) ' transformed tuvP '
+
+END_PROVIDER
+

src/casscf/casscf.irp.f

@@ -19,7 +19,13 @@ subroutine run
     N_det = 1
     TOUCH N_det psi_det psi_coef 
     call run_cipsi
-    call driver_wdens
+
+    write(6,*) ' total energy = ',eone+etwo+ecore
+    mo_label = "MCSCF"
+    mo_label = "Natural"
+    mo_coef(:,:) = NatOrbsFCI(:,:)
+    call save_mos
+
     call driver_optorb
     energy_old = energy
     energy = eone+etwo+ecore

src/casscf/densities.irp.f

@@ -1,56 +1,30 @@
-! -*- F90 -*-
+use bitmasks
-use bitmasks ! you need to include the bitmasks_module.f90 features
 
 BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
-&BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
   BEGIN_DOC
   ! the first-order density matrix in the basis of the starting MOs
-! the second-order density matrix in the basis of the starting MOs
   ! matrices are state averaged
   !
   ! we use the spin-free generators of mono-excitations
   ! E_pq destroys q and creates p
   ! D_pq   =     <0|E_pq|0> = D_qp
-! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
   !
   END_DOC
   implicit none
   integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
   integer                        :: ierr
-      integer(bit_kind), allocatable :: det_mu(:,:)
+  integer(bit_kind)              :: det_mu(N_int,2)
-      integer(bit_kind), allocatable :: det_mu_ex(:,:)
+  integer(bit_kind)              :: det_mu_ex(N_int,2)
-      integer(bit_kind), allocatable :: det_mu_ex1(:,:)
+  integer(bit_kind)              :: det_mu_ex1(N_int,2)
-      integer(bit_kind), allocatable :: det_mu_ex11(:,:)
+  integer(bit_kind)              :: det_mu_ex2(N_int,2)
-      integer(bit_kind), allocatable :: det_mu_ex12(:,:)
+  real*8                         :: phase1,phase2,term
-      integer(bit_kind), allocatable :: det_mu_ex2(:,:)
+  integer                        :: nu1,nu2
-      integer(bit_kind), allocatable :: det_mu_ex21(:,:)
+  integer                        :: ierr1,ierr2
-      integer(bit_kind), allocatable :: det_mu_ex22(:,:)
+  real*8                         :: cI_mu(N_states)
-      real*8 :: phase1,phase11,phase12,phase2,phase21,phase22
-      integer :: nu1,nu2,nu11,nu12,nu21,nu22
-      integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22
-      real*8 :: cI_mu(N_states),term
-      allocate(det_mu(N_int,2))
-      allocate(det_mu_ex(N_int,2))
-      allocate(det_mu_ex1(N_int,2))
-      allocate(det_mu_ex11(N_int,2))
-      allocate(det_mu_ex12(N_int,2))
-      allocate(det_mu_ex2(N_int,2))
-      allocate(det_mu_ex21(N_int,2))
-      allocate(det_mu_ex22(N_int,2))
   
   write(6,*) ' providing density matrices D0 and P0 '
   
-! set all to zero
+  D0tu = 0.d0
-       do t=1,n_act_orb
-        do u=1,n_act_orb
-         D0tu(u,t)=0.D0
-         do v=1,n_act_orb
-          do x=1,n_act_orb
-           P0tuvx(x,v,u,t)=0.D0
-          end do
-         end do
-        end do
-       end do
   
   ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
   do mu=1,n_det
@@ -72,6 +46,73 @@ END_DOC
           do istate=1,n_states
             term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
             D0tu(t,u)+=term
+          end do
+        end if
+        ! det_mu_ex2 is in the list
+        if (nu2.ne.-1) then
+          do istate=1,n_states
+            term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
+            D0tu(t,u)+=term
+          end do
+        end if
+      end do
+    end do
+  end do
+  
+  ! we average by just dividing by the number of states
+  do x=1,n_act_orb
+    do v=1,n_act_orb
+      D0tu(v,x)*=1.0D0/dble(N_states)
+    end do
+  end do
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
+   BEGIN_DOC
+   ! the second-order density matrix in the basis of the starting MOs
+   ! matrices are state averaged
+   !
+   ! we use the spin-free generators of mono-excitations
+   ! E_pq destroys q and creates p
+   ! D_pq   =     <0|E_pq|0> = D_qp
+   ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
+   !
+   END_DOC
+   implicit none
+   integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
+   integer                        :: ierr
+   real*8                         :: phase1,phase11,phase12,phase2,phase21,phase22
+   integer                        :: nu1,nu2,nu11,nu12,nu21,nu22
+   integer                        :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22
+   real*8                         :: cI_mu(N_states),term
+   integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex
+   integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12
+   integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
+   
+   write(6,*) ' providing density matrices D0 and P0 '
+
+   P0tuvx = 0.d0
+   
+   ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
+   do mu=1,n_det
+     call det_extract(det_mu,mu,N_int)
+     do istate=1,n_states
+       cI_mu(istate)=psi_coef(mu,istate)
+     end do
+     do t=1,n_act_orb
+       ipart=list_act(t)
+       do u=1,n_act_orb
+         ihole=list_act(u)
+         ! apply E_tu
+         call det_copy(det_mu,det_mu_ex1,N_int)
+         call det_copy(det_mu,det_mu_ex2,N_int)
+         call do_spinfree_mono_excitation(det_mu,det_mu_ex1          &
+             ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
+         ! det_mu_ex1 is in the list
+         if (nu1.ne.-1) then
+           do istate=1,n_states
+             term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
              ! and we fill P0_tvvu
              do v=1,n_act_orb
                P0tuvx(t,v,v,u)-=term
@@ -82,7 +123,6 @@ END_DOC
          if (nu2.ne.-1) then
            do istate=1,n_states
              term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
-           D0tu(t,u)+=term
              do v=1,n_act_orb
                P0tuvx(t,v,v,u)-=term
              end do
@@ -165,7 +205,6 @@ END_DOC
    ! we average by just dividing by the number of states
    do x=1,n_act_orb
      do v=1,n_act_orb
-        D0tu(v,x)*=1.0D0/dble(N_states)
        do u=1,n_act_orb
          do t=1,n_act_orb
            P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)

src/casscf/det_manip.irp.f

@@ -1,5 +1,4 @@
-! -*- F90 -*-
+use bitmasks
-use bitmasks ! you need to include the bitmasks_module.f90 features
 
 subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
       ispin,phase,ierr)

src/casscf/driver_wdens.irp.f

@@ -1,154 +0,0 @@
-       subroutine driver_wdens
-         implicit none
-         integer :: istate,p,q,r,s,indx,i,j
-
-
-         write(6,*) ' total energy = ',eone+etwo+ecore
-         write(6,*) ' generating natural orbitals '
-         write(6,*)
-         write(6,*)
-         call trf_to_natorb
-
-         write(6,*) ' all data available ! '
-         write(6,*) '    writing out files '
-
-         open(unit=12,file='D0tu.dat',form='formatted',status='unknown')
-         do p=1,n_act_orb
-          do q=1,n_act_orb
-           if (abs(D0tu(p,q)).gt.1.D-12) then
-            write(12,'(2i8,E20.12)') p,q,D0tu(p,q)
-           end if
-          end do
-         end do 
-         close(12)
-         
-real*8 :: approx,np,nq,nr,ns
-logical :: lpq,lrs,lps,lqr
-         open(unit=12,file='P0tuvx.dat',form='formatted',status='unknown')
-         do p=1,n_act_orb
-          np=D0tu(p,p)
-          do q=1,n_act_orb
-           lpq=p.eq.q
-           nq=D0tu(q,q)
-           do r=1,n_act_orb
-            lqr=q.eq.r
-            nr=D0tu(r,r)
-            do s=1,n_act_orb
-             lrs=r.eq.s
-             lps=p.eq.s
-             approx=0.D0
-             if (lpq.and.lrs) then
-              if (lqr) then
-! pppp
-               approx=0.5D0*np*(np-1.D0)
-              else
-! pprr
-               approx=0.5D0*np*nr
-              end if
-             else
-              if (lps.and.lqr.and..not.lpq) then
-! pqqp
-               approx=-0.25D0*np*nq
-              end if
-             end if
-             if (abs(P0tuvx(p,q,r,s)).gt.1.D-12) then
-              write(12,'(4i4,2E20.12)') p,q,r,s,P0tuvx(p,q,r,s),approx
-             end if
-            end do
-           end do
-          end do 
-         end do 
-         close(12)
-         
-         open(unit=12,form='formatted',status='unknown',file='onetrf.tmp')
-         indx=0
-         do q=1,mo_num
-          do p=q,mo_num
-            if (abs(onetrf(p,q)).gt.1.D-12) then
-             write(12,'(2i6,E20.12)') p,q,onetrf(p,q)
-             indx+=1
-            end if
-          end do 
-         end do 
-         write(6,*) ' wrote ',indx,' mono-electronic integrals'
-         close(12)
- 
-         
-         open(unit=12,form='formatted',status='unknown',file='bielec_PQxx.tmp')
-         indx=0
-         do p=1,mo_num
-          do q=p,mo_num
-           do r=1,n_core_orb+n_act_orb
-            do s=r,n_core_orb+n_act_orb
-             if (abs(bielec_PQxxtmp(p,q,r,s)).gt.1.D-12) then
-              write(12,'(4i8,E20.12)') p,q,r,s,bielec_PQxxtmp(p,q,r,s)
-              indx+=1
-             end if
-            end do
-           end do
-          end do 
-         end do 
-         write(6,*) ' wrote ',indx,' integrals (PQ|xx)'
-         close(12)
- 
-         open(unit=12,form='formatted',status='unknown',file='bielec_PxxQ.tmp')
-         indx=0
-         do p=1,mo_num
-          do q=1,n_core_orb+n_act_orb
-           do r=q,n_core_orb+n_act_orb
-integer ::s_start
-            if (q.eq.r) then
-             s_start=p
-            else
-             s_start=1
-            end if
-            do s=s_start,mo_num
-             if (abs(bielec_PxxQtmp(p,q,r,s)).gt.1.D-12) then
-              write(12,'(4i8,E20.12)') p,q,r,s,bielec_PxxQtmp(p,q,r,s)
-              indx+=1
-             end if
-            end do
-           end do
-          end do 
-         end do 
-         write(6,*) ' wrote ',indx,' integrals (Px|xQ)'
-         close(12)
- 
-         open(unit=12,form='formatted',status='unknown',file='bielecCI.tmp')
-         indx=0
-         do p=1,n_act_orb
-          do q=p,n_act_orb
-           do r=1,n_act_orb
-            do s=1,mo_num
-             if (abs(bielecCItmp(p,q,r,s)).gt.1.D-12) then
-              write(12,'(4i8,E20.12)') p,q,r,s,bielecCItmp(p,q,r,s)
-              indx+=1
-             end if
-            end do
-           end do
-          end do 
-         end do 
-         write(6,*) ' wrote ',indx,' integrals (tu|xP)'
-         close(12)
-
-         write(6,*)
-         write(6,*) '   creating new orbitals '
-         do i=1,mo_num
-          write(6,*) ' Orbital No ',i
-          write(6,'(5F14.6)') (NatOrbsFCI(j,i),j=1,mo_num)
-          write(6,*)
-         end do
-
-         mo_label = "MCSCF"
-         mo_label = "Natural"
-         do i=1,mo_num
-          do j=1,ao_num
-           mo_coef(j,i)=NatOrbsFCI(j,i)
-          end do
-         end do
-         call save_mos
- 
-         write(6,*) '   ... done '
-
-       end
-

src/casscf/gradient.irp.f

@@ -1,10 +1,8 @@
-! -*- F90 -*-
+use bitmasks
-
-use bitmasks ! you need to include the bitmasks_module.f90 features
 
 BEGIN_PROVIDER [ integer, nMonoEx ]
   BEGIN_DOC
-! 
+  ! Number of single excitations
   END_DOC
   implicit none
   nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
@@ -205,7 +203,7 @@ END_DOC
       x3=x+n_core_orb
       do y=1,n_act_orb
         y3=y+n_core_orb
-            gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx(ii,vv,x3,y3)
+        gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
       end do
     end do
   end do
@@ -241,7 +239,7 @@ END_DOC
   do v=1,n_act_orb
     do x=1,n_act_orb
       do y=1,n_act_orb
-            gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa)
+        gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa)
       end do
     end do
   end do

src/casscf/hessian.irp.f

@@ -1,6 +1,4 @@
-! -*- F90 -*-
+use bitmasks
-
-use bitmasks ! you need to include the bitmasks_module.f90 features
 
 BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
   BEGIN_DOC
@@ -304,19 +302,19 @@ END_DOC
       ! diagonal element
       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(tt,i,i,tt)-bielec_pqxx(tt,tt,i,i))
+      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(tt,i,i,tt) &
+      term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt)             &
-              -bielec_pqxx(tt,tt,i,i))
+          -bielec_pqxx_no(tt,tt,i,i))
       term-=occnum(tt)*Fipq(tt,tt)
       do v=1,n_act_orb
         vv=list_act(v)
         do x=1,n_act_orb
           xx=list_act(x)
-           term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(vv,xx,i,i) &
+          term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(vv,xx,i,i)      &
               +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*              &
-                bielec_pxxq(vv,i,i,xx))
+              bielec_pxxq_no(vv,i,i,xx))
           do y=1,n_act_orb
-            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx)
+            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
           end do
         end do
       end do
@@ -324,21 +322,21 @@ END_DOC
       ! it/iu, t != u
       uu=list_act(u)
       term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu))
-         term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) &
+      term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) &
-            -bielec_PQxx(tt,uu,i,j)) 
+          -bielec_PQxx_no(tt,uu,i,j))
       term-=occnum(tt)*Fipq(uu,tt)
       term-=(occnum(tt)+occnum(uu))                                  &
-           *(3.D0*bielec_PxxQ(tt,i,i,uu)-bielec_PQxx(uu,tt,i,i))
+          *(3.D0*bielec_PxxQ_no(tt,i,i,uu)-bielec_PQxx_no(uu,tt,i,i))
       do v=1,n_act_orb
         vv=list_act(v)
         !         term-=D0tu(u,v)*Fipq(tt,vv)   ! published, but inverting t and u seems more correct
         do x=1,n_act_orb
           xx=list_act(x)
-           term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,i) &
+          term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,i)      &
               +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v))               &
-                *bielec_pxxq(vv,i,i,xx))
+              *bielec_pxxq_no(vv,i,i,xx))
           do y=1,n_act_orb
-            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(u,v,y,xx)
+            term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(u,v,y,xx)
           end do
         end do
       end do
@@ -355,18 +353,18 @@ END_DOC
     else
       term=0.D0
     end if
-        term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) &
+    term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt)   &
-            -bielec_PQxx(tt,uu,i,j)) 
+        -bielec_PQxx_no(tt,uu,i,j))
     term-=(occnum(tt)+occnum(uu))*                                   &
-           (4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) &
+        (4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt)          &
-           -bielec_PQxx(uu,tt,i,j))
+        -bielec_PQxx_no(uu,tt,i,j))
     do v=1,n_act_orb
       vv=list_act(v)
       do x=1,n_act_orb
         xx=list_act(x)
-          term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,j) &
+        term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,j)        &
             +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v))                 &
-               *bielec_pxxq(vv,i,j,xx))
+            *bielec_pxxq_no(vv,i,j,xx))
       end do
     end do
   end if
@@ -390,17 +388,17 @@ END_DOC
   tt=list_act(t)
   jj=list_core(j)
   aa=list_virt(a)
-       term=2.D0*(4.D0*bielec_pxxq(aa,j,i,tt) &
+  term=2.D0*(4.D0*bielec_pxxq_no(aa,j,i,tt)                             &
-          -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt))
+      -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt))
-       term-=occnum(tt)*(4.D0*bielec_pxxq(aa,j,i,tt) &
+  term-=occnum(tt)*(4.D0*bielec_pxxq_no(aa,j,i,tt)                      &
-          -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt))
+      -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt))
   if (i.eq.j) then
     term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt))
     term-=0.5D0*occnum(tt)*Fipq(aa,tt)
     do v=1,n_act_orb
       do x=1,n_act_orb
         do y=1,n_act_orb
-           term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa)
+          term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa)
         end do
       end do
     end do
@@ -430,8 +428,8 @@ END_DOC
   else
     term=0.D0
   end if
-       term-=occnum(uu)*(bielec_pqxx(aa,ii,t3,u3)-4.D0*bielec_pqxx(aa,uu,t3,i) &
+  term-=occnum(uu)*(bielec_pqxx_no(aa,ii,t3,u3)-4.D0*bielec_pqxx_no(aa,uu,t3,i)&
-                +bielec_pxxq(aa,t3,u3,ii))
+      +bielec_pxxq_no(aa,t3,u3,ii))
   do v=1,n_act_orb
     vv=list_act(v)
     v3=v+n_core_orb
@@ -439,9 +437,9 @@ END_DOC
       integer                        :: x3
       xx=list_act(x)
       x3=x+n_core_orb
-         term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,ii,v3,x3) &
+      term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,ii,v3,x3)        &
           +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u))                   &
-              *bielec_pqxx(aa,xx,v3,i))
+          *bielec_pqxx_no(aa,xx,v3,i))
     end do
   end do
   if (t.eq.u) then
@@ -467,19 +465,19 @@ END_DOC
     if (a.eq.b) then
       ! ia/ia
       term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii))
-         term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,aa)-bielec_pqxx(aa,aa,i,i))
+      term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,aa)-bielec_pqxx_no(aa,aa,i,i))
     else
       bb=list_virt(b)
       ! ia/ib
       term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb))
-         term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,bb)-bielec_pqxx(aa,bb,i,i))
+      term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,bb)-bielec_pqxx_no(aa,bb,i,i))
     end if
   else
     ! ia/jb
     jj=list_core(j)
     bb=list_virt(b)
-        term=2.D0*(4.D0*bielec_pxxq(aa,i,j,bb)-bielec_pqxx(aa,bb,i,j) &
+    term=2.D0*(4.D0*bielec_pxxq_no(aa,i,j,bb)-bielec_pqxx_no(aa,bb,i,j)    &
-         -bielec_pxxq(aa,j,i,bb))
+        -bielec_pxxq_no(aa,j,i,bb))
     if (a.eq.b) then
       term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj))
     end if
@@ -503,15 +501,15 @@ END_DOC
   tt=list_act(t)
   bb=list_virt(b)
   t3=t+n_core_orb
-       term=occnum(tt)*(4.D0*bielec_pxxq(aa,i,t3,bb)-bielec_pxxq(aa,t3,i,bb) &
+  term=occnum(tt)*(4.D0*bielec_pxxq_no(aa,i,t3,bb)-bielec_pxxq_no(aa,t3,i,bb)&
-             -bielec_pqxx(aa,bb,i,t3))
+      -bielec_pqxx_no(aa,bb,i,t3))
   if (a.eq.b) then
     term-=Fipq(tt,ii)+Fapq(tt,ii)
     term-=0.5D0*occnum(tt)*Fipq(tt,ii)
     do v=1,n_act_orb
       do x=1,n_act_orb
         do y=1,n_act_orb
-           term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,ii)
+          term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,ii)
         end do
       end do
     end do
@@ -545,11 +543,11 @@ END_DOC
         do x=1,n_act_orb
           xx=list_act(x)
           x3=x+n_core_orb
-           t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,aa,v3,x3) &
+          t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3)      &
               +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*              &
-                bielec_pxxq(aa,x3,v3,aa))
+              bielec_pxxq_no(aa,x3,v3,aa))
           do y=1,n_act_orb
-            t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx)
+            t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
           end do
         end do
       end do
@@ -565,9 +563,9 @@ END_DOC
         do x=1,n_act_orb
           xx=list_act(x)
           x3=x+n_core_orb
-           term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,bb,v3,x3) &
+          term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3)    &
               +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))               &
-                *bielec_pxxq(aa,x3,v3,bb))
+              *bielec_pxxq_no(aa,x3,v3,bb))
         end do
       end do
     end if
@@ -582,9 +580,9 @@ END_DOC
       do x=1,n_act_orb
         xx=list_act(x)
         x3=x+n_core_orb
-          term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,bb,v3,x3) &
+        term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3)      &
             +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v))                 &
-               *bielec_pxxq(aa,x3,v3,bb))
+            *bielec_pxxq_no(aa,x3,v3,bb))
       end do
     end do
     if (a.eq.b) then
@@ -592,8 +590,8 @@ END_DOC
       do v=1,n_act_orb
         do x=1,n_act_orb
           do y=1,n_act_orb
-            term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,uu)
+            term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu)
-            term-=P0tuvx_no(u,v,x,y)*bielecCI(x,y,v,tt)
+            term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt)
           end do
         end do
       end do

src/casscf/mcscf_fock.irp.f

@@ -1,9 +1,41 @@
-! -*- F90 -*- 
 BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
-&BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
    BEGIN_DOC
-! the inactive and the active Fock matrices, in molecular 
+   ! the inactive Fock matrix, in molecular orbitals
-! orbitals
+   END_DOC
+   implicit none
+   integer                        :: p,q,k,kk,t,tt,u,uu
+   
+   do q=1,mo_num
+     do p=1,mo_num
+       Fipq(p,q)=one_ints_no(p,q)
+     end do
+   end do
+   
+   ! the inactive Fock matrix
+   do k=1,n_core_orb
+     kk=list_core(k)
+     do q=1,mo_num
+       do p=1,mo_num
+         Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q)
+       end do
+     end do
+   end do
+   
+   if (bavard) then
+     integer                        :: i
+     write(6,*)
+     write(6,*) ' the diagonal of the inactive effective Fock matrix '
+     write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
+     write(6,*)
+   end if
+   
+   
+END_PROVIDER
+ 
+ 
+BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
+   BEGIN_DOC
+   ! the active active Fock matrix, in molecular orbitals
    ! we create them in MOs, quite expensive
    !
    ! for an implementation in AOs we need first the natural orbitals
@@ -11,36 +43,17 @@ BEGIN_DOC
    !
    END_DOC
    implicit none
-       double precision, allocatable :: integrals_array1(:,:)
-       double precision, allocatable :: integrals_array2(:,:)
    integer                        :: p,q,k,kk,t,tt,u,uu
-       allocate(integrals_array1(mo_num,mo_num))
-       allocate(integrals_array2(mo_num,mo_num))
    
-       do p=1,mo_num
+   Fapq = 0.d0
-        do q=1,mo_num
-         Fipq(p,q)=one_ints(p,q)
-         Fapq(p,q)=0.D0
-        end do
-       end do
-
-! the inactive Fock matrix
-       do k=1,n_core_orb
-        kk=list_core(k)
-        do p=1,mo_num
-         do q=1,mo_num
-          Fipq(p,q)+=2.D0*bielec_pqxx(p,q,k,k) -bielec_pxxq(p,k,k,q)
-         end do
-        end do
-       end do
    
    ! the active Fock matrix, D0tu is diagonal
    do t=1,n_act_orb
      tt=list_act(t)
-        do p=1,mo_num
      do q=1,mo_num
+       do p=1,mo_num
          Fapq(p,q)+=occnum(tt)                                       &
-              *(bielec_pqxx(p,q,tt,tt)-0.5D0*bielec_pxxq(p,tt,tt,q))
+             *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q))
        end do
      end do
    end do

src/casscf/natorb.irp.f

@@ -1,21 +1,40 @@
-! -*- F90 -*-
+ BEGIN_PROVIDER [real*8, occnum, (mo_num)]
-! diagonalize D0tu
-! save the diagonal somewhere, in inverse order
-! 4-index-transform the 2-particle density matrix over active orbitals
-! correct the bielectronic integrals
-! correct the monoelectronic integrals
-! put integrals on file, as well orbitals, and the density matrices
-!  
-     subroutine trf_to_natorb
   implicit none
-      integer :: i,j,k,l,t,u,p,q,pp
+  BEGIN_DOC
-      real*8 :: eigval(n_act_orb),natorbsCI(n_act_orb,n_act_orb)
+  ! MO occupation numbers
-      real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
+  END_DOC
+  
+  integer :: i
+  occnum=0.D0
+  do i=1,n_core_orb
+    occnum(list_core(i))=2.D0
+  end do
+  
+  do i=1,n_act_orb
+    occnum(list_act(i))=occ_act(n_act_orb-i+1)
+  end do
+
+  write(6,*) ' occupation numbers '
+  do i=1,mo_num
+    write(6,*) i,occnum(i)
+  end do
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ]
+&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ]
+ implicit none
+ BEGIN_DOC
+ ! Natural orbitals of CI
+ END_DOC
+ integer                        :: i, j
+ 
+ call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
  
-      call lapack_diag(eigval,natorbsCI,D0tu,n_act_orb,n_act_orb)
  write(6,*) ' found occupation numbers as '
  do i=1,n_act_orb
-       write(6,*) i,eigval(i)
+   write(6,*) i,occ_act(i)
  end do
  
  if (bavard) then
@@ -38,24 +57,26 @@ real*8 :: xmx
        natOrbsCI(j,i)*=xmx
      end do
      
-         
      write(6,*) ' Eigenvector No ',i
      write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb)
    end do
  end if
   
-      do i=1,n_act_orb
+END_PROVIDER
-       do j=1,n_act_orb
+
-        D0tu(i,j)=0.D0
+
-       end do
+BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
-! fill occupation numbers in descending order
+  implicit none
-       D0tu(i,i)=eigval(n_act_orb-i+1)
+  BEGIN_DOC
-      end do
-!
   ! 4-index transformation of 2part matrices
-!
+  END_DOC
+  integer                        :: i,j,k,l,p,q,pp
+  real*8                         :: d(n_act_orb)
+
   ! index per index
   ! first quarter
+  P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:)
+
   do j=1,n_act_orb
     do k=1,n_act_orb
       do l=1,n_act_orb
@@ -65,11 +86,11 @@ real*8 :: xmx
         do p=1,n_act_orb
           pp=n_act_orb-p+1
           do q=1,n_act_orb
-            d(pp)+=P0tuvx(q,j,k,l)*natorbsCI(q,p)
+            d(pp)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
           end do
         end do
         do p=1,n_act_orb
-           P0tuvx(p,j,k,l)=d(p)
+          P0tuvx_no(p,j,k,l)=d(p)
         end do
       end do
     end do
@@ -84,11 +105,11 @@ real*8 :: xmx
         do p=1,n_act_orb
           pp=n_act_orb-p+1
           do q=1,n_act_orb
-            d(pp)+=P0tuvx(j,q,k,l)*natorbsCI(q,p)
+            d(pp)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
           end do
         end do
         do p=1,n_act_orb
-           P0tuvx(j,p,k,l)=d(p)
+          P0tuvx_no(j,p,k,l)=d(p)
         end do
       end do
     end do
@@ -103,11 +124,11 @@ real*8 :: xmx
         do p=1,n_act_orb
           pp=n_act_orb-p+1
           do q=1,n_act_orb
-            d(pp)+=P0tuvx(j,k,q,l)*natorbsCI(q,p)
+            d(pp)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
           end do
         end do
         do p=1,n_act_orb
-           P0tuvx(j,k,p,l)=d(p)
+          P0tuvx_no(j,k,p,l)=d(p)
         end do
       end do
     end do
@@ -122,24 +143,30 @@ real*8 :: xmx
         do p=1,n_act_orb
           pp=n_act_orb-p+1
           do q=1,n_act_orb
-            d(pp)+=P0tuvx(j,k,l,q)*natorbsCI(q,p)
+            d(pp)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
           end do
         end do
         do p=1,n_act_orb
-           P0tuvx(j,k,l,p)=d(p)
+          P0tuvx_no(j,k,l,p)=d(p)
         end do
       end do
     end do
   end do
   write(6,*) ' transformed P0tuvx '
-!
+
-! one-electron integrals 
+END_PROVIDER
-!
+
-       do i=1,mo_num
+
-        do j=1,mo_num
+
-         onetrf(i,j)=mo_one_e_integrals(i,j)
+BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
-        end do
+  implicit none
-       end do
+  BEGIN_DOC
+  ! Transformed one-e integrals
+  END_DOC
+  integer :: i,j, p, pp, q
+  real*8 :: d(n_act_orb)
+  one_ints_no(:,:)=mo_one_e_integrals(:,:)
+
   ! 1st half-trf
   do j=1,mo_num
     do p=1,n_act_orb
@@ -148,13 +175,14 @@ real*8 :: xmx
     do p=1,n_act_orb
       pp=n_act_orb-p+1
       do q=1,n_act_orb
-          d(pp)+=onetrf(list_act(q),j)*natorbsCI(q,p)
+        d(pp)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
       end do
     end do
     do p=1,n_act_orb
-          onetrf(list_act(p),j)=d(p)
+      one_ints_no(list_act(p),j)=d(p)
     end do
   end do
+
   ! 2nd half-trf
   do j=1,mo_num
     do p=1,n_act_orb
@@ -163,22 +191,26 @@ real*8 :: xmx
     do p=1,n_act_orb
       pp=n_act_orb-p+1
       do q=1,n_act_orb
-          d(pp)+=onetrf(j,list_act(q))*natorbsCI(q,p)
+        d(pp)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
       end do
     end do
     do p=1,n_act_orb
-         onetrf(j,list_act(p))=d(p)
+      one_ints_no(j,list_act(p))=d(p)
-        end do
-       end do
-       write(6,*) ' transformed onetrf '
-! 
-! Orbitals
-!
-       do j=1,ao_num
-        do i=1,mo_num
-         NatOrbsFCI(j,i)=mo_coef(j,i)
     end do
   end do
+  write(6,*) ' transformed one_ints '
+END_PROVIDER
+
+
+BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
+  implicit none
+  BEGIN_DOC
+! FCI natural orbitals
+  END_DOC
+  integer :: i,j, p, pp, q
+  real*8 :: d(n_act_orb)
+
+  NatOrbsFCI(:,:)=mo_coef(:,:)
   
   do j=1,ao_num
     do p=1,n_act_orb
@@ -195,229 +227,25 @@ real*8 :: xmx
     end do
   end do
   write(6,*) ' transformed orbitals '
-!
+END_PROVIDER
-! now the bielectronic integrals
-!
-!!$       write(6,*) ' before the transformation '
-!!$integer :: kk,ll,ii,jj
-!!$real*8 :: h1,h2,h3
-!!$       do i=1,n_act_orb
-!!$          ii=list_act(i)
-!!$          do j=1,n_act_orb
-!!$             jj=list_act(j)
-!!$             do k=1,n_act_orb
-!!$                kk=list_act(k)
-!!$                do l=1,n_act_orb
-!!$                   ll=list_act(l)
-!!$                   h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb)
-!!$                   h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll)
-!!$                   h3=bielecCItmp(i,j,k,ll)
-!!$                   if ((h1.ne.h2).or.(h1.ne.h3)) then
-!!$                      write(6,9901) i,j,k,l,h1,h2,h3
-!!$9901                  format(' aie   ',4i4,3E20.12)
-!!$9902                  format('correct',4i4,3E20.12)
-!!$                   else
-!!$                      write(6,9902) i,j,k,l,h1,h2,h3
-!!$                   end if
-!!$                end do
-!!$             end do
-!!$          end do
-!!$       end do
 
-       do j=1,mo_num
+
-        do k=1,n_core_orb+n_act_orb
+
-         do l=1,n_core_orb+n_act_orb
+
-          do p=1,n_act_orb
+
-           d1(p)=0.D0
+
-           d2(p)=0.D0
+subroutine trf_to_natorb()
-          end do
+  implicit none
-          do p=1,n_act_orb
+  BEGIN_DOC
-           pp=n_act_orb-p+1
+  ! save the diagonal somewhere, in inverse order
-           do q=1,n_act_orb
+  ! 4-index-transform the 2-particle density matrix over active orbitals
-            d1(pp)+=bielec_PQxxtmp(list_act(q),j,k,l)*natorbsCI(q,p)
+  ! correct the bielectronic integrals
-            d2(pp)+=bielec_PxxQtmp(list_act(q),k,l,j)*natorbsCI(q,p)
+  ! correct the monoelectronic integrals
-           end do
+  ! put integrals on file, as well orbitals, and the density matrices
-          end do
-          do p=1,n_act_orb
-           bielec_PQxxtmp(list_act(p),j,k,l)=d1(p)
-           bielec_PxxQtmp(list_act(p),k,l,j)=d2(p)
-          end do
-         end do
-        end do
-       end do
-! 2nd quarter
-       do j=1,mo_num
-        do k=1,n_core_orb+n_act_orb
-         do l=1,n_core_orb+n_act_orb
-          do p=1,n_act_orb
-           d1(p)=0.D0
-           d2(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d1(pp)+=bielec_PQxxtmp(j,list_act(q),k,l)*natorbsCI(q,p)
-            d2(pp)+=bielec_PxxQtmp(j,k,l,list_act(q))*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielec_PQxxtmp(j,list_act(p),k,l)=d1(p)
-           bielec_PxxQtmp(j,k,l,list_act(p))=d2(p)
-          end do
-         end do
-        end do
-       end do
-! 3rd quarter
-       do j=1,mo_num
-        do k=1,mo_num
-         do l=1,n_core_orb+n_act_orb
-          do p=1,n_act_orb
-           d1(p)=0.D0
-           d2(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d1(pp)+=bielec_PQxxtmp(j,k,n_core_orb+q,l)*natorbsCI(q,p)
-            d2(pp)+=bielec_PxxQtmp(j,n_core_orb+q,l,k)*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielec_PQxxtmp(j,k,n_core_orb+p,l)=d1(p)
-           bielec_PxxQtmp(j,n_core_orb+p,l,k)=d2(p)
-          end do
-         end do
-        end do
-       end do
-! 4th quarter
-       do j=1,mo_num
-        do k=1,mo_num
-         do l=1,n_core_orb+n_act_orb
-          do p=1,n_act_orb
-           d1(p)=0.D0
-           d2(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d1(pp)+=bielec_PQxxtmp(j,k,l,n_core_orb+q)*natorbsCI(q,p)
-            d2(pp)+=bielec_PxxQtmp(j,l,n_core_orb+q,k)*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielec_PQxxtmp(j,k,l,n_core_orb+p)=d1(p)
-           bielec_PxxQtmp(j,l,n_core_orb+p,k)=d2(p)
-          end do
-         end do
-        end do
-       end do
-       write(6,*) ' transformed PQxx and PxxQ '
   !
-! and finally the bielecCI integrals
+  END_DOC
-!
+  integer                        :: i,j,k,l,t,u,p,q,pp
-       do j=1,n_act_orb
+  real*8                         :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
-        do k=1,n_act_orb
-         do l=1,mo_num
-          do p=1,n_act_orb
-           d(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d(pp)+=bielecCItmp(q,j,k,l)*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielecCItmp(p,j,k,l)=d(p)
-          end do
-         end do
-        end do
-       end do
-! 2nd quarter
-       do j=1,n_act_orb
-        do k=1,n_act_orb
-         do l=1,mo_num
-          do p=1,n_act_orb
-           d(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d(pp)+=bielecCItmp(j,q,k,l)*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielecCItmp(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,mo_num
-          do p=1,n_act_orb
-           d(p)=0.D0
-          end do
-          do p=1,n_act_orb
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d(pp)+=bielecCItmp(j,k,q,l)*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielecCItmp(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
-           pp=n_act_orb-p+1
-           do q=1,n_act_orb
-            d(pp)+=bielecCItmp(j,k,l,list_act(q))*natorbsCI(q,p)
-           end do
-          end do
-          do p=1,n_act_orb
-           bielecCItmp(j,k,l,list_act(p))=d(p)
-          end do
-         end do
-        end do
-       end do
-       write(6,*) ' transformed tuvP '
-!
-! that's all
-!
-!!$
-!!$! test coherence of the bielectronic integals
-!!$!   PQxx = PxxQ = tuvP for some of the indices
-!!$       write(6,*) ' after the transformation '
-!!$       do i=1,n_act_orb
-!!$          ii=list_act(i)
-!!$          do j=1,n_act_orb
-!!$             jj=list_act(j)
-!!$             do k=1,n_act_orb
-!!$                kk=list_act(k)
-!!$                do l=1,n_act_orb
-!!$                   ll=list_act(l)
-!!$                   h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb)
-!!$                   h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll)
-!!$                   h3=bielecCItmp(i,j,k,ll)
-!!$                   if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
-!!$                      write(6,9901) i,j,k,l,h1,h1-h2,h1-h3
-!!$                   else
-!!$                      write(6,9902) i,j,k,l,h1,h2,h3
-!!$                   end if
-!!$                end do
-!!$             end do
-!!$          end do
-!!$       end do
   
   ! we recalculate total energies
   write(6,*)
@@ -443,32 +271,31 @@ real*8 :: xmx
   e_two_all=0.D0
   do i=1,n_core_orb
     ii=list_core(i)
-        e_one_all+=2.D0*onetrf(ii,ii)
+    e_one_all+=2.D0*one_ints_no(ii,ii)
     do j=1,n_core_orb
       jj=list_core(j)
-         e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
+      e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
     end do
     do t=1,n_act_orb
       tt=list_act(t)
       t3=t+n_core_orb
-         do u=1,n_act_orb
+      e_two_all += occnum(list_act(t)) *                            &
-          uu=list_act(u)
+          (2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3))
-          u3=u+n_core_orb
-          e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
-              -bielec_PQxxtmp(tt,ii,i,u3))
-         end do
     end do
   end do
+
+    
+
   do t=1,n_act_orb
     tt=list_act(t)
+    e_one_all += occnum(list_act(t))*one_ints_no(tt,tt)
     do u=1,n_act_orb
       uu=list_act(u)
-         e_one_all+=D0tu(t,u)*onetrf(tt,uu)
       do v=1,n_act_orb
         v3=v+n_core_orb
         do x=1,n_act_orb
           x3=x+n_core_orb
-           e_two_all  +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3)
+          e_two_all  +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3)
         end do
       end do
     end do
@@ -479,12 +306,12 @@ real*8 :: xmx
   ecore_bis=nuclear_repulsion
   do i=1,n_core_orb
     ii=list_core(i)
-        ecore    +=2.D0*onetrf(ii,ii)
+    ecore    +=2.D0*one_ints_no(ii,ii)
-        ecore_bis+=2.D0*onetrf(ii,ii)
+    ecore_bis+=2.D0*one_ints_no(ii,ii)
     do j=1,n_core_orb
       jj=list_core(j)
-         ecore    +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
+      ecore    +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
-         ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii)
+      ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii)
     end do
   end do
   eone    =0.D0
@@ -495,18 +322,18 @@ real*8 :: xmx
   do t=1,n_act_orb
     tt=list_act(t)
     t3=t+n_core_orb
+    eone     += occnum(list_act(t))*one_ints_no(tt,tt)
+    eone_bis += occnum(list_act(t))*one_ints_no(tt,tt)
+    do i=1,n_core_orb
+      ii=list_core(i)
+      eone     += occnum(list_act(t)) * &
+        (2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3))
+      eone_bis += occnum(list_act(t)) * &
+        (2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt))
+    end do
     do u=1,n_act_orb
       uu=list_act(u)
       u3=u+n_core_orb
-         eone    +=D0tu(t,u)*onetrf(tt,uu)
-         eone_bis+=D0tu(t,u)*onetrf(tt,uu)
-         do i=1,n_core_orb
-          ii=list_core(i)
-          eone    +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
-            -bielec_PQxxtmp(tt,ii,i,u3))
-          eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) &
-            -bielec_PxxQtmp(tt,i,i,uu))
-         end do
       do v=1,n_act_orb
         vv=list_act(v)
         v3=v+n_core_orb
@@ -514,12 +341,12 @@ real*8 :: xmx
           xx=list_act(x)
           x3=x+n_core_orb
           real*8                         :: h1,h2,h3
-           h1=bielec_PQxxtmp(tt,uu,v3,x3)
+          h1=bielec_PQxx_no(tt,uu,v3,x3)
-           h2=bielec_PxxQtmp(tt,u3,v3,xx)
+          h2=bielec_PxxQ_no(tt,u3,v3,xx)
-           h3=bielecCItmp(t,u,v,xx)
+          h3=bielecCI_no(t,u,v,xx)
-           etwo    +=P0tuvx(t,u,v,x)*h1
+          etwo    +=P0tuvx_no(t,u,v,x)*h1
-           etwo_bis+=P0tuvx(t,u,v,x)*h2
+          etwo_bis+=P0tuvx_no(t,u,v,x)*h2
-           etwo_ter+=P0tuvx(t,u,v,x)*h3
+          etwo_ter+=P0tuvx_no(t,u,v,x)*h3
           if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
             write(6,9901) t,u,v,x,h1,h2,h3
             9901 format('aie: ',4I4,3E20.12)
@@ -540,9 +367,7 @@ real*8 :: h1,h2,h3
   write(6,*) ' ----------------------------------------- '
   write(6,*) '     sum of all        = ',eone+etwo+ecore
   write(6,*)
+  SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter 
   
 end subroutine trf_to_natorb
 
- BEGIN_PROVIDER [real*8, onetrf, (mo_num,mo_num)]
-&BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
-END_PROVIDER

src/casscf/natorb_casscf.irp.f

@@ -1,65 +0,0 @@
-! -*- F90 -*-
-BEGIN_PROVIDER [real*8, occnum, (mo_num)]
-    implicit none
-    integer :: i,kk,j
-    logical :: lread
-    real*8 :: rdum
-    do i=1,mo_num
-     occnum(i)=0.D0
-    end do
-    do i=1,n_core_orb
-     occnum(list_core(i))=2.D0
-    end do
-
-    open(unit=12,file='D0tu.dat',form='formatted',status='old')
-    lread=.true.
-    do while (lread)
-     read(12,*,iostat=kk) i,j,rdum
-     if (kk.ne.0) then
-      lread=.false.
-     else
-      if (i.eq.j) then
-       occnum(list_act(i))=rdum
-      else
-       write(6,*) ' WARNING - no natural orbitals !'
-       write(6,*) i,j,rdum
-      end if
-     end if
-    end do
-    close(12)
-    write(6,*) ' read occupation numbers '
-    do i=1,mo_num
-     write(6,*) i,occnum(i)
-    end do
-
-END_PROVIDER
-
-BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
-      implicit none
-      integer :: i,j,k,l,kk
-      real*8 :: rdum
-      logical :: lread
-
-      do i=1,n_act_orb
-       do j=1,n_act_orb
-        do k=1,n_act_orb
-         do l=1,n_act_orb
-          P0tuvx_no(l,k,j,i)=0.D0
-         end do
-        end do
-       end do
-      end do
-
-      open(unit=12,file='P0tuvx.dat',form='formatted',status='old')
-      lread=.true.
-      do while (lread)
-       read(12,*,iostat=kk) i,j,k,l,rdum
-       if (kk.ne.0) then
-        lread=.false.
-       else
-        P0tuvx_no(i,j,k,l)=rdum
-       end if
-      end do
-      close(12)
-      write(6,*) ' read the 2-particle density matrix '
-END_PROVIDER

src/casscf/one_ints.irp.f

@@ -1,26 +0,0 @@
-! -*- F90 -*-
-BEGIN_PROVIDER [real*8, one_ints, (mo_num,mo_num)]
-      implicit none 
-      integer :: i,j,kk
-      logical :: lread
-      real*8 :: rdum
-      do i=1,mo_num
-       do j=1,mo_num
-        one_ints(i,j)=0.D0
-       end do 
-      end do
-      open(unit=12,file='onetrf.tmp',status='old',form='formatted')
-      lread=.true.
-      do while (lread)
-       read(12,*,iostat=kk) i,j,rdum
-       if (kk.ne.0) then
-        lread=.false.
-       else
-        one_ints(i,j)=rdum
-        one_ints(j,i)=rdum
-       end if
-      end do
-      close(12)
-      write(6,*) ' read MCSCF natural one-electron integrals '
-END_PROVIDER
-

src/casscf/tot_en.irp.f

@@ -1,4 +1,3 @@
-! -*- F90 -*-
  BEGIN_PROVIDER [real*8, etwo]
 &BEGIN_PROVIDER [real*8, eone]
 &BEGIN_PROVIDER [real*8, eone_bis]
@@ -16,7 +15,7 @@ real*8 :: e_one_all,e_two_all
      e_one_all+=2.D0*mo_one_e_integrals(ii,ii)
      do j=1,n_core_orb
        jj=list_core(j)
-         e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
+       e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
      end do
      do t=1,n_act_orb
        tt=list_act(t)
@@ -24,8 +23,8 @@ real*8 :: e_one_all,e_two_all
        do u=1,n_act_orb
          uu=list_act(u)
          u3=u+n_core_orb
-          e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
+         e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i)           &
-              -bielec_PQxxtmp(tt,ii,i,u3))
+             -bielec_PQxx(tt,ii,i,u3))
        end do
      end do
    end do
@@ -38,7 +37,7 @@ real*8 :: e_one_all,e_two_all
          v3=v+n_core_orb
          do x=1,n_act_orb
            x3=x+n_core_orb
-           e_two_all  +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3)
+           e_two_all  +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3)
          end do
        end do
      end do
@@ -53,8 +52,8 @@ real*8 :: e_one_all,e_two_all
      ecore_bis+=2.D0*mo_one_e_integrals(ii,ii)
      do j=1,n_core_orb
        jj=list_core(j)
-         ecore    +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
+       ecore    +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
-         ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii)
+       ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii)
      end do
    end do
    eone    =0.D0
@@ -72,10 +71,10 @@ real*8 :: e_one_all,e_two_all
        eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
        do i=1,n_core_orb
          ii=list_core(i)
-          eone    +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
+         eone    +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i)            &
-            -bielec_PQxxtmp(tt,ii,i,u3))
+             -bielec_PQxx(tt,ii,i,u3))
-          eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) &
+         eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii)           &
-            -bielec_PxxQtmp(tt,i,i,uu))
+             -bielec_PxxQ(tt,i,i,uu))
        end do
        do v=1,n_act_orb
          vv=list_act(v)
@@ -84,9 +83,9 @@ real*8 :: e_one_all,e_two_all
            xx=list_act(x)
            x3=x+n_core_orb
            real*8                         :: h1,h2,h3
-           h1=bielec_PQxxtmp(tt,uu,v3,x3)
+           h1=bielec_PQxx(tt,uu,v3,x3)
-           h2=bielec_PxxQtmp(tt,u3,v3,xx)
+           h2=bielec_PxxQ(tt,u3,v3,xx)
-           h3=bielecCItmp(t,u,v,xx)
+           h3=bielecCI(t,u,v,xx)
            etwo    +=P0tuvx(t,u,v,x)*h1
            etwo_bis+=P0tuvx(t,u,v,x)*h2
            etwo_ter+=P0tuvx(t,u,v,x)*h3