Using fast 2RDM s

src/casscf/densities.irp.f

@@ -29,7 +29,9 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
    !
    END_DOC
    implicit none
-   integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
+   integer                        :: t,u,v,x
+   integer                        :: tt,uu,vv,xx
+   integer                        :: 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
@@ -43,125 +45,25 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
     write(6,*) ' providing density matrix P0'
   endif
 
-   P0tuvx = 0.d0
+  P0tuvx= 0.d0
-   
+  do istate=1,N_states
-   ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
+   do x = 1, n_act_orb
-   do mu=1,n_det
+    xx = list_act(x)
-     call det_extract(det_mu,mu,N_int)
+    do v = 1, n_act_orb
-     do istate=1,n_states
+     vv = list_act(v)
-       cI_mu(istate)=psi_coef(mu,istate)
+     do u = 1, n_act_orb
-     end do
+      uu = list_act(u)
-     do t=1,n_act_orb
+      do t = 1, n_act_orb
-       ipart=list_act(t)
+       tt = list_act(t)
-       do u=1,n_act_orb
+       P0tuvx(t,u,v,x) =                                   &
-         ihole=list_act(u)
+         state_average_weight(istate) *                    &
-         ! apply E_tu
+         ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) +    &
-         call det_copy(det_mu,det_mu_ex1,N_int)
+           two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) +    &
-         call det_copy(det_mu,det_mu_ex2,N_int)
+           two_rdm_beta_beta_mo  (tt,uu,vv,xx,istate) )
-         call do_spinfree_mono_excitation(det_mu,det_mu_ex1          &
+      enddo
-             ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
+     enddo 
-         ! det_mu_ex1 is in the list
+    enddo
-         if (nu1.ne.-1) then
+   enddo
-           do istate=1,n_states
+  enddo
-             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
-             end do
-           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
-             do v=1,n_act_orb
-               P0tuvx(t,v,v,u)-=term
-             end do
-           end do
-         end if
-       end do
-     end do
-   end do
-   ! now we do the double excitation E_tu E_vx |0>
-   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 v=1,n_act_orb
-       ipart=list_act(v)
-       do x=1,n_act_orb
-         ihole=list_act(x)
-         ! apply E_vx
-         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)
-         ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0>
-         if (ierr1.eq.1) then
-           do t=1,n_act_orb
-             jpart=list_act(t)
-             do u=1,n_act_orb
-               jhole=list_act(u)
-               call det_copy(det_mu_ex1,det_mu_ex11,N_int)
-               call det_copy(det_mu_ex1,det_mu_ex12,N_int)
-               call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11&
-                   ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12)
-               if (nu11.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)&
-                       *phase11*phase1
-                 end do
-               end if
-               if (nu12.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)&
-                       *phase12*phase1
-                 end do
-               end if
-             end do
-           end do
-         end if
-         
-         ! we apply E_tu to the second resultant determinant
-         if (ierr2.eq.1) then
-           do t=1,n_act_orb
-             jpart=list_act(t)
-             do u=1,n_act_orb
-               jhole=list_act(u)
-               call det_copy(det_mu_ex2,det_mu_ex21,N_int)
-               call det_copy(det_mu_ex2,det_mu_ex22,N_int)
-               call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21&
-                   ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22)
-               if (nu21.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)&
-                       *phase21*phase2
-                 end do
-               end if
-               if (nu22.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)&
-                       *phase22*phase2
-                 end do
-               end if
-             end do
-           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
-       do u=1,n_act_orb
-         do t=1,n_act_orb
-           P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)
-         end do
-       end do
-     end do
-   end do
 
 END_PROVIDER

src/casscf/test_two_rdm.irp.f

@@ -8,11 +8,11 @@ program print_two_rdm
 
   double precision :: accu,twodm
   accu = 0.d0
-  do i=1,mo_num
+  do i=1,n_act_orb
-    do j=1,mo_num
+    do j=1,n_act_orb
-      do k=1,mo_num
+      do k=1,n_act_orb
-        do l=1,mo_num
+        do l=1,n_act_orb
-         twodm = coussin_peter_two_rdm_mo(i,j,k,l,1)
+         twodm = coussin_peter_two_rdm_mo(list_act(i),list_act(j),list_act(k),list_act(l))
             if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then
              print*,''
              print*,'sum'

src/two_body_rdm/two_rdm.irp.f

@@ -1,23 +1,27 @@
-
+BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num)]
- BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
  implicit none
  BEGIN_DOC
  !  coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF 
  END_DOC
- integer :: i,j,k,l
+ integer :: i,j,k,l, istate
+ coussin_peter_two_rdm_mo = 0.d0
+ do istate=1,N_states
   do l = 1, mo_num
    do k = 1, mo_num
     do j = 1, mo_num
      do i = 1, mo_num
-     coussin_peter_two_rdm_mo(i,j,k,l,:) = 0.5d0 * (two_rdm_alpha_beta_mo(i,j,k,l,:) + two_rdm_alpha_beta_mo(i,j,k,l,:)) & 
+      coussin_peter_two_rdm_mo(i,j,k,l) =            &
-                                         +          two_rdm_alpha_alpha_mo(i,j,k,l,:)                                  &
+        state_average_weight(istate) *               &
-                                         +          two_rdm_beta_beta_mo(i,j,k,l,:)         
+        ( two_rdm_alpha_beta_mo(i,j,k,l,istate) +    &
+          two_rdm_alpha_alpha_mo(i,j,k,l,istate)+    &
+          two_rdm_beta_beta_mo(i,j,k,l,istate) )
+     enddo
     enddo 
    enddo
   enddo
  enddo
   
- END_PROVIDER 
+END_PROVIDER 
 
 
  BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]