Merge branch 'casscf' of github.com:QuantumPackage/qp2 into casscf

src/bitmask/core_inact_act_virt.irp.f

@@ -141,6 +141,10 @@ END_PROVIDER
  n_act_orb_tmp = 0
  n_virt_orb_tmp = 0
  n_del_orb_tmp = 0
+ core_bitmask = 0_bit_kind
+ inact_bitmask = 0_bit_kind
+ act_bitmask = 0_bit_kind
+ virt_bitmask = 0_bit_kind
  do i = 1, mo_num
   if(mo_class(i) == 'Core')then
    n_core_orb_tmp += 1

src/casscf/bavard.irp.f

@@ -1,6 +1,6 @@
 ! -*- F90 -*- 
 BEGIN_PROVIDER [logical, bavard]
-     bavard=.true.
+!     bavard=.true.
-!     bavard=.false.
+     bavard=.false.
 END_PROVIDER
 

src/casscf/casscf.irp.f

@@ -4,7 +4,8 @@ program casscf
 ! TODO : Put the documentation of the program here
   END_DOC
   no_vvvv_integrals = .True.
-  SOFT_TOUCH no_vvvv_integrals
+  pt2_max = 0.02
+  SOFT_TOUCH no_vvvv_integrals pt2_max
   call run
 end
 
@@ -19,8 +20,7 @@ subroutine run
   mo_label = "MCSCF"
   iteration = 1
   do while (.not.converged)
-    call run_cipsi
+    call run_stochastic_cipsi
-
     energy_old = energy
     energy = eone+etwo+ecore
 
@@ -30,14 +30,19 @@ subroutine run
     call write_double(6,energy_improvement, 'Predicted energy improvement')
 
     converged = dabs(energy_improvement) < thresh_scf
+    pt2_max = dabs(energy_improvement / pt2_relative_error)
 
     mo_coef = NewOrbs
     call save_mos
     call map_deinit(mo_integrals_map)
-    N_det = 1
     iteration += 1
-    FREE mo_integrals_map mo_two_e_integrals_in_map psi_det psi_coef
+    N_det = N_det/2
-    SOFT_TOUCH mo_coef N_det
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    read_wf = .True.
+    FREE mo_integrals_map mo_two_e_integrals_in_map
+    SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef
+
   enddo
 
 end

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/neworbs.irp.f

@@ -25,7 +25,7 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
   end do
   
   if (bavard) then
-    do i=2,nMonoEx+1
+    do i=2,nMonoEx
       write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
     end do
   end if
@@ -77,14 +77,14 @@ END_PROVIDER
   
   energy_improvement = SXeigenval(best_vector)
   
+  c0=SXeigenvec(1,best_vector)
+
   if (bavard) then
     write(6,*) ' SXdiag : eigenvalue for best overlap with '
     write(6,*) '  previous orbitals = ',SXeigenval(best_vector)
     write(6,*) ' weight of the 1st element ',c0
   endif
 
-  c0=SXeigenvec(1,best_vector)
-
   do i=1,nMonoEx+1
     SXvector(i)=SXeigenvec(i,best_vector)/c0
   end do

src/casscf/test_two_rdm.irp.f

@@ -1,30 +0,0 @@
-program print_two_rdm
-  implicit none
-  integer                        :: i,j,k,l
-  read_wf = .True.
-  TOUCH read_wf
-
-  double precision, parameter    :: thr = 1.d-15
-
-  double precision :: accu,twodm
-  accu = 0.d0
-  do i=1,mo_num
-    do j=1,mo_num
-      do k=1,mo_num
-        do l=1,mo_num
-         twodm = coussin_peter_two_rdm_mo(i,j,k,l,1)
-            if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then
-             print*,''
-             print*,'sum'
-             write(*,'(3X,4(I2,X),3(F16.13,X))'),  i, j, k, l, twodm,P0tuvx(i,j,k,l),dabs(twodm - P0tuvx(i,j,k,l))
-             print*,''
-            endif
-            accu += dabs(twodm - P0tuvx(i,j,k,l))
-        enddo
-      enddo
-    enddo
-  enddo
-  print*,'accu = ',accu
-  print*,'<accu> ',accu / dble(mo_num**4)
-
-end

src/cipsi/cipsi.irp.f

@@ -13,6 +13,7 @@ subroutine run_cipsi
   rss = memory_of_double(N_states)*4.d0
   call check_mem(rss,irp_here)
 
+  N_iter = 1
   allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))
 
   double precision               :: hf_energy_ref

src/cipsi/pt2_stoch_routines.irp.f

@@ -135,7 +135,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
     PROVIDE psi_occ_pattern_hii det_to_occ_pattern
   endif
 
-  if (N_det < max(4,N_states)) then
+  if (N_det <= max(4,N_states)) then
     pt2=0.d0
     variance=0.d0
     norm=0.d0
@@ -719,6 +719,15 @@ END_PROVIDER
 
   double precision :: rss
   double precision, external :: memory_of_double, memory_of_int
+  if (N_det_generators == 1) then
+    pt2_w = 1.d0
+    pt2_cw = 1.d0
+    pt2_W_T = 1.d0
+    pt2_u_0 = 1.d0
+    pt2_n_0 = 1
+    return
+  endif
+
   rss = memory_of_double(2*N_det_generators+1)
   call check_mem(rss,irp_here)
 
@@ -754,7 +763,7 @@ END_PROVIDER
     end if
     pt2_n_0(1) += 1
     if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then
-      stop "teeth building failed"
+      print *, "teeth building failed"
     end if
   end do
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

src/cipsi/stochastic_cipsi.irp.f

@@ -12,6 +12,7 @@ subroutine run_stochastic_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated N_generators_bitmask
 
+  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 

src/generators_cas/generators.irp.f

@@ -55,6 +55,7 @@ END_PROVIDER
       nongen(inongen) = i
     endif
   enddo
+  ASSERT (m == N_det_generators)
 
   psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators)
   psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :)

src/two_body_rdm/README.rst

@@ -2,5 +2,7 @@
 two_body_rdm
 ============
 
-Contains the two rdms (aa,bb,ab) stored as plain arrays 
+Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
+maps, with pysicists notation, consistent with the two-electron integrals in the
+MO basis.
 

src/two_body_rdm/all_2rdm_routines.irp.f

@@ -1,443 +1,442 @@
-
+subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze)
- subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze)
+   use bitmasks
-  use bitmasks
+   implicit none
-  implicit none
+   BEGIN_DOC
-  BEGIN_DOC
+   ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
-  ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
+   !
-  !
+   ! Assumes that the determinants are in psi_det
-  ! Assumes that the determinants are in psi_det
+   !
-  !
+   ! istart, iend, ishift, istep are used in ZMQ parallelization.
-  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+   END_DOC
-  END_DOC
+   integer, intent(in)            :: N_st,sze
-  integer, intent(in)            :: N_st,sze
+   integer, intent(in)            :: dim1,dim2,dim3,dim4
-  integer, intent(in) :: dim1,dim2,dim3,dim4
+   double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: u_0(sze,N_st)
-  double precision, intent(inout) :: u_0(sze,N_st)
+   integer                        :: k
-  integer :: k
+   double precision, allocatable  :: u_t(:,:)
-  double precision, allocatable  :: u_t(:,:)
+   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
+   allocate(u_t(N_st,N_det))
-  allocate(u_t(N_st,N_det))
+   do k=1,N_st
-  do k=1,N_st
+     call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
-    call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
+   enddo
-  enddo
+   call dtranspose(                                                  &
-  call dtranspose(                                                   &
+       u_0,                                                          &
-      u_0,                                                           &
+       size(u_0, 1),                                                 &
-      size(u_0, 1),                                                  &
+       u_t,                                                          &
-      u_t,                                                           &
+       size(u_t, 1),                                                 &
-      size(u_t, 1),                                                  &
+       N_det, N_st)
-      N_det, N_st)
+   
-
+   call all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
-  call all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
+   deallocate(u_t)
-  deallocate(u_t)
+   
-
+   do k=1,N_st
-  do k=1,N_st
+     call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
-    call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
+   enddo
-  enddo
+   
-
+end
- end
+ 
-
+ 
-
+subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
- subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+   use bitmasks
-  use bitmasks
+   implicit none
-  implicit none
+   BEGIN_DOC
-  BEGIN_DOC
+   ! Computes two-rdm
-  ! Computes two-rdm
+   !
-  !
+   ! Default should be 1,N_det,0,1
-  ! Default should be 1,N_det,0,1
+   END_DOC
-  END_DOC
+   integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
-  integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
+   integer, intent(in)            :: dim1,dim2,dim3,dim4
-  integer, intent(in) :: dim1,dim2,dim3,dim4
+   double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(in)   :: u_t(N_st,N_det)
-  double precision, intent(in)   :: u_t(N_st,N_det)
+   
-
+   
-  
+   PROVIDE N_int
-  PROVIDE N_int 
+   
-
+   select case (N_int)
-  select case (N_int)
+     case (1)
-    case (1)
+       call all_two_rdm_dm_nstates_openmp_work_1(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-      call all_two_rdm_dm_nstates_openmp_work_1(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+     case (2)
-    case (2)
+       call all_two_rdm_dm_nstates_openmp_work_2(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-      call all_two_rdm_dm_nstates_openmp_work_2(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+     case (3)
-    case (3)
+       call all_two_rdm_dm_nstates_openmp_work_3(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-      call all_two_rdm_dm_nstates_openmp_work_3(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+     case (4)
-    case (4)
+       call all_two_rdm_dm_nstates_openmp_work_4(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-      call all_two_rdm_dm_nstates_openmp_work_4(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+       case default
-    case default
+       call all_two_rdm_dm_nstates_openmp_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-      call all_two_rdm_dm_nstates_openmp_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+   end select
-  end select
+end
- end
+ 
-
+ BEGIN_TEMPLATE
-BEGIN_TEMPLATE
+ 
-
 subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
-  use bitmasks
+   use bitmasks
-  implicit none
+   implicit none
-  BEGIN_DOC
+   BEGIN_DOC
-  ! Computes $v_t = H | u_t \\rangle$ and $s_t = S^2  | u_t \\rangle$
+   ! Computes $v_t = H | u_t \\rangle$ and $s_t = S^2  | u_t \\rangle$
-  !
+   !
-  ! Default should be 1,N_det,0,1
+   ! Default should be 1,N_det,0,1
-  END_DOC
+   END_DOC
-  integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
+   integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
-  double precision, intent(in)   :: u_t(N_st,N_det)
+   double precision, intent(in)   :: u_t(N_st,N_det)
-  integer, intent(in) :: dim1,dim2,dim3,dim4
+   integer, intent(in)            :: dim1,dim2,dim3,dim4
-  double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
-  double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
+   double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
-
+   
-  integer                        :: i,j,k,l
+   integer                        :: i,j,k,l
-  integer                        :: k_a, k_b, l_a, l_b, m_a, m_b
+   integer                        :: k_a, k_b, l_a, l_b, m_a, m_b
-  integer                        :: istate
+   integer                        :: istate
-  integer                        :: krow, kcol, krow_b, kcol_b
+   integer                        :: krow, kcol, krow_b, kcol_b
-  integer                        :: lrow, lcol
+   integer                        :: lrow, lcol
-  integer                        :: mrow, mcol
+   integer                        :: mrow, mcol
-  integer(bit_kind)              :: spindet($N_int)
+   integer(bit_kind)              :: spindet($N_int)
-  integer(bit_kind)              :: tmp_det($N_int,2)
+   integer(bit_kind)              :: tmp_det($N_int,2)
-  integer(bit_kind)              :: tmp_det2($N_int,2)
+   integer(bit_kind)              :: tmp_det2($N_int,2)
-  integer(bit_kind)              :: tmp_det3($N_int,2)
+   integer(bit_kind)              :: tmp_det3($N_int,2)
-  integer(bit_kind), allocatable :: buffer(:,:)
+   integer(bit_kind), allocatable :: buffer(:,:)
-  integer                        :: n_doubles
+   integer                        :: n_doubles
-  integer, allocatable           :: doubles(:)
+   integer, allocatable           :: doubles(:)
-  integer, allocatable           :: singles_a(:)
+   integer, allocatable           :: singles_a(:)
-  integer, allocatable           :: singles_b(:)
+   integer, allocatable           :: singles_b(:)
-  integer, allocatable           :: idx(:), idx0(:)
+   integer, allocatable           :: idx(:), idx0(:)
-  integer                        :: maxab, n_singles_a, n_singles_b, kcol_prev
+   integer                        :: maxab, n_singles_a, n_singles_b, kcol_prev
-  integer*8                      :: k8
+   integer*8                      :: k8
-
+   
-  maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
+   maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
-  allocate(idx0(maxab))
+   allocate(idx0(maxab))
-
+   
-  do i=1,maxab
+   do i=1,maxab
-    idx0(i) = i
+     idx0(i) = i
-  enddo
+   enddo
-
+   
-  ! Prepare the array of all alpha single excitations
+   ! Prepare the array of all alpha single excitations
-  ! -------------------------------------------------
+   ! -------------------------------------------------
-
+   
-  PROVIDE N_int nthreads_davidson
+   PROVIDE N_int nthreads_davidson
- !!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)        &
+   !!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
- !    !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,                &
+       !    !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,          &
- !    !$OMP          psi_bilinear_matrix_columns,                    &
+       !    !$OMP          psi_bilinear_matrix_columns,              &
- !    !$OMP          psi_det_alpha_unique, psi_det_beta_unique,      &
+       !    !$OMP          psi_det_alpha_unique, psi_det_beta_unique,&
- !    !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,   &
+       !    !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,&
- !    !$OMP          psi_bilinear_matrix_transp_rows,                &
+       !    !$OMP          psi_bilinear_matrix_transp_rows,          &
- !    !$OMP          psi_bilinear_matrix_transp_columns,             &
+       !    !$OMP          psi_bilinear_matrix_transp_columns,       &
- !    !$OMP          psi_bilinear_matrix_transp_order, N_st,         &
+       !    !$OMP          psi_bilinear_matrix_transp_order, N_st,   &
- !    !$OMP          psi_bilinear_matrix_order_transp_reverse,       &
+       !    !$OMP          psi_bilinear_matrix_order_transp_reverse, &
- !    !$OMP          psi_bilinear_matrix_columns_loc,                &
+       !    !$OMP          psi_bilinear_matrix_columns_loc,          &
- !    !$OMP          psi_bilinear_matrix_transp_rows_loc,            &
+       !    !$OMP          psi_bilinear_matrix_transp_rows_loc,      &
- !    !$OMP          istart, iend, istep, irp_here, v_t, s_t,        &
+       !    !$OMP          istart, iend, istep, irp_here, v_t, s_t,  &
- !    !$OMP          ishift, idx0, u_t, maxab)                       &
+       !    !$OMP          ishift, idx0, u_t, maxab)                 &
- !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
+       !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,&
- !    !$OMP          lcol, lrow, l_a, l_b,                           &
+       !    !$OMP          lcol, lrow, l_a, l_b,                     &
- !    !$OMP          buffer, doubles, n_doubles,                     &
+       !    !$OMP          buffer, doubles, n_doubles,               &
- !    !$OMP          tmp_det2, idx, l, kcol_prev,          &
+       !    !$OMP          tmp_det2, idx, l, kcol_prev,              &
- !    !$OMP          singles_a, n_singles_a, singles_b,              &
+       !    !$OMP          singles_a, n_singles_a, singles_b,        &
- !    !$OMP          n_singles_b, k8)
+       !    !$OMP          n_singles_b, k8)
-
+   
-  ! Alpha/Beta double excitations
+   ! Alpha/Beta double excitations
-  ! =============================
+   ! =============================
-
+   
-  allocate( buffer($N_int,maxab),                                     &
+   allocate( buffer($N_int,maxab),                                   &
-      singles_a(maxab),                                              &
+       singles_a(maxab),                                             &
-      singles_b(maxab),                                              &
+       singles_b(maxab),                                             &
-      doubles(maxab),                                                &
+       doubles(maxab),                                               &
-      idx(maxab))
+       idx(maxab))
-
+   
-  kcol_prev=-1
+   kcol_prev=-1
-
+   
-  ASSERT (iend <= N_det)
+   ASSERT (iend <= N_det)
-  ASSERT (istart > 0)
+   ASSERT (istart > 0)
-  ASSERT (istep  > 0)
+   ASSERT (istep  > 0)
-
+   
- !!$OMP DO SCHEDULE(dynamic,64)
+   !!$OMP DO SCHEDULE(dynamic,64)
-  do k_a=istart+ishift,iend,istep
+   do k_a=istart+ishift,iend,istep
-
+     
-    krow = psi_bilinear_matrix_rows(k_a)
+     krow = psi_bilinear_matrix_rows(k_a)
-    ASSERT (krow <= N_det_alpha_unique)
+     ASSERT (krow <= N_det_alpha_unique)
-
+     
-    kcol = psi_bilinear_matrix_columns(k_a)
+     kcol = psi_bilinear_matrix_columns(k_a)
-    ASSERT (kcol <= N_det_beta_unique)
+     ASSERT (kcol <= N_det_beta_unique)
-
+     
-    tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
+     tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-    tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+     tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-
+     
-    if (kcol /= kcol_prev) then
+     if (kcol /= kcol_prev) then
-      call get_all_spin_singles_$N_int(                              &
+       call get_all_spin_singles_$N_int(                             &
-          psi_det_beta_unique, idx0,                                 &
+           psi_det_beta_unique, idx0,                                &
-          tmp_det(1,2), N_det_beta_unique,                           &
+           tmp_det(1,2), N_det_beta_unique,                          &
-          singles_b, n_singles_b)
+           singles_b, n_singles_b)
-    endif
+     endif
-    kcol_prev = kcol
+     kcol_prev = kcol
-
+     
-    ! Loop over singly excited beta columns
+     ! Loop over singly excited beta columns
-    ! -------------------------------------
+     ! -------------------------------------
-
+     
-    do i=1,n_singles_b
+     do i=1,n_singles_b
-      lcol = singles_b(i)
+       lcol = singles_b(i)
-
+       
-      tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
+       tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
-
+       
-      l_a = psi_bilinear_matrix_columns_loc(lcol)
+       l_a = psi_bilinear_matrix_columns_loc(lcol)
-      ASSERT (l_a <= N_det)
+       ASSERT (l_a <= N_det)
-
+       
-      do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
+       do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
-        lrow = psi_bilinear_matrix_rows(l_a)
+         lrow = psi_bilinear_matrix_rows(l_a)
-        ASSERT (lrow <= N_det_alpha_unique)
+         ASSERT (lrow <= N_det_alpha_unique)
-
+         
-        buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
+         buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
-
+         
-        ASSERT (l_a <= N_det)
+         ASSERT (l_a <= N_det)
-        idx(j) = l_a
+         idx(j) = l_a
-        l_a = l_a+1
+         l_a = l_a+1
-      enddo
+       enddo
-      j = j-1
+       j = j-1
-
+       
-      call get_all_spin_singles_$N_int(                              &
+       call get_all_spin_singles_$N_int(                             &
-          buffer, idx, tmp_det(1,1), j,                              &
+           buffer, idx, tmp_det(1,1), j,                             &
-          singles_a, n_singles_a )
+           singles_a, n_singles_a )
-
+       
-      ! Loop over alpha singles
+       ! Loop over alpha singles
-      ! -----------------------
+       ! -----------------------
-
+       
-      do k = 1,n_singles_a
+       do k = 1,n_singles_a
-        l_a = singles_a(k)
+         l_a = singles_a(k)
-        ASSERT (l_a <= N_det)
+         ASSERT (l_a <= N_det)
-
+         
-        lrow = psi_bilinear_matrix_rows(l_a)
+         lrow = psi_bilinear_matrix_rows(l_a)
-        ASSERT (lrow <= N_det_alpha_unique)
+         ASSERT (lrow <= N_det_alpha_unique)
-
+         
-        tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
+         tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
-       !call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij)
+         !call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij)
-        do l= 1, N_states
+         do l= 1, N_states
+           c_1(l) = u_t(l,l_a)
+           c_2(l) = u_t(l,k_a)
+         enddo
+         call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
+       enddo
+       
+     enddo
+     
+   enddo
+   ! !$OMP END DO
+   
+   ! !$OMP DO SCHEDULE(dynamic,64)
+   do k_a=istart+ishift,iend,istep
+     
+     
+     ! Single and double alpha exitations
+     ! ===================================
+     
+     
+     ! Initial determinant is at k_a in alpha-major representation
+     ! -----------------------------------------------------------------------
+     
+     krow = psi_bilinear_matrix_rows(k_a)
+     ASSERT (krow <= N_det_alpha_unique)
+     
+     kcol = psi_bilinear_matrix_columns(k_a)
+     ASSERT (kcol <= N_det_beta_unique)
+     
+     tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
+     tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+     
+     ! Initial determinant is at k_b in beta-major representation
+     ! ----------------------------------------------------------------------
+     
+     k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
+     ASSERT (k_b <= N_det)
+     
+     spindet(1:$N_int) = tmp_det(1:$N_int,1)
+     
+     ! Loop inside the beta column to gather all the connected alphas
+     lcol = psi_bilinear_matrix_columns(k_a)
+     l_a = psi_bilinear_matrix_columns_loc(lcol)
+     do i=1,N_det_alpha_unique
+       if (l_a > N_det) exit
+       lcol = psi_bilinear_matrix_columns(l_a)
+       if (lcol /= kcol) exit
+       lrow = psi_bilinear_matrix_rows(l_a)
+       ASSERT (lrow <= N_det_alpha_unique)
+       
+       buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
+       idx(i) = l_a
+       l_a = l_a+1
+     enddo
+     i = i-1
+     
+     call get_all_spin_singles_and_doubles_$N_int(                   &
+         buffer, idx, spindet, i,                                    &
+         singles_a, doubles, n_singles_a, n_doubles )
+     
+     ! Compute Hij for all alpha singles
+     ! ----------------------------------
+     
+     tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+     do i=1,n_singles_a
+       l_a = singles_a(i)
+       ASSERT (l_a <= N_det)
+       
+       lrow = psi_bilinear_matrix_rows(l_a)
+       ASSERT (lrow <= N_det_alpha_unique)
+       
+       tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
+       do l= 1, N_states
          c_1(l) = u_t(l,l_a)
          c_2(l) = u_t(l,k_a)
-        enddo
+       enddo
-        call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
+       ! increment the alpha/beta part for single excitations
-      enddo
+       call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
-
+       ! increment the alpha/alpha part for single excitations
-    enddo
+       call off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
-
+       
-  enddo
+     enddo
-! !$OMP END DO
+     
-
+     
-! !$OMP DO SCHEDULE(dynamic,64)
+     ! Compute Hij for all alpha doubles
-  do k_a=istart+ishift,iend,istep
+     ! ----------------------------------
-
+     
-
+     do i=1,n_doubles
-    ! Single and double alpha exitations
+       l_a = doubles(i)
-    ! ===================================
+       ASSERT (l_a <= N_det)
-
+       
-
+       lrow = psi_bilinear_matrix_rows(l_a)
-    ! Initial determinant is at k_a in alpha-major representation
+       ASSERT (lrow <= N_det_alpha_unique)
-    ! -----------------------------------------------------------------------
+       
-
+       do l= 1, N_states
-    krow = psi_bilinear_matrix_rows(k_a)
+         c_1(l) = u_t(l,l_a)
-    ASSERT (krow <= N_det_alpha_unique)
+         c_2(l) = u_t(l,k_a)
-
+       enddo
-    kcol = psi_bilinear_matrix_columns(k_a)
+       call off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
-    ASSERT (kcol <= N_det_beta_unique)
+     enddo
-
+     
-    tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
+     
-    tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+     ! Single and double beta excitations
-
+     ! ==================================
-    ! Initial determinant is at k_b in beta-major representation
+     
-    ! ----------------------------------------------------------------------
+     
-
+     ! Initial determinant is at k_a in alpha-major representation
-    k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
+     ! -----------------------------------------------------------------------
-    ASSERT (k_b <= N_det)
+     
-
+     krow = psi_bilinear_matrix_rows(k_a)
-    spindet(1:$N_int) = tmp_det(1:$N_int,1)
+     kcol = psi_bilinear_matrix_columns(k_a)
-
+     
-    ! Loop inside the beta column to gather all the connected alphas
+     tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-    lcol = psi_bilinear_matrix_columns(k_a)
+     tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-    l_a = psi_bilinear_matrix_columns_loc(lcol)
+     
-    do i=1,N_det_alpha_unique
+     spindet(1:$N_int) = tmp_det(1:$N_int,2)
-      if (l_a > N_det) exit
+     
-      lcol = psi_bilinear_matrix_columns(l_a)
+     ! Initial determinant is at k_b in beta-major representation
-      if (lcol /= kcol) exit
+     ! -----------------------------------------------------------------------
-      lrow = psi_bilinear_matrix_rows(l_a)
+     
-      ASSERT (lrow <= N_det_alpha_unique)
+     k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
-
+     ASSERT (k_b <= N_det)
-      buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
+     
-      idx(i) = l_a
+     ! Loop inside the alpha row to gather all the connected betas
-      l_a = l_a+1
+     lrow = psi_bilinear_matrix_transp_rows(k_b)
-    enddo
+     l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
-    i = i-1
+     do i=1,N_det_beta_unique
-
+       if (l_b > N_det) exit
-    call get_all_spin_singles_and_doubles_$N_int(                    &
+       lrow = psi_bilinear_matrix_transp_rows(l_b)
-        buffer, idx, spindet, i,                                     &
+       if (lrow /= krow) exit
-        singles_a, doubles, n_singles_a, n_doubles )
+       lcol = psi_bilinear_matrix_transp_columns(l_b)
-
+       ASSERT (lcol <= N_det_beta_unique)
-    ! Compute Hij for all alpha singles
+       
-    ! ----------------------------------
+       buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
-
+       idx(i) = l_b
-    tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+       l_b = l_b+1
-    do i=1,n_singles_a
+     enddo
-      l_a = singles_a(i)
+     i = i-1
-      ASSERT (l_a <= N_det)
+     
-
+     call get_all_spin_singles_and_doubles_$N_int(                   &
-      lrow = psi_bilinear_matrix_rows(l_a)
+         buffer, idx, spindet, i,                                    &
-      ASSERT (lrow <= N_det_alpha_unique)
+         singles_b, doubles, n_singles_b, n_doubles )
-
+     
-      tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
+     ! Compute Hij for all beta singles
-      do l= 1, N_states
+     ! ----------------------------------
-       c_1(l) = u_t(l,l_a)
+     
-       c_2(l) = u_t(l,k_a)
+     tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-      enddo
+     do i=1,n_singles_b
-      ! increment the alpha/beta part for single excitations 
+       l_b = singles_b(i)
-      call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
+       ASSERT (l_b <= N_det)
-      ! increment the alpha/alpha part for single excitations 
+       
-      call off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
+       lcol = psi_bilinear_matrix_transp_columns(l_b)
-
+       ASSERT (lcol <= N_det_beta_unique)
-    enddo
+       
-
+       tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
-
+       l_a = psi_bilinear_matrix_transp_order(l_b)
-    ! Compute Hij for all alpha doubles
+       do l= 1, N_states
-    ! ----------------------------------
+         c_1(l) = u_t(l,l_a)
-
+         c_2(l) = u_t(l,k_a)
-    do i=1,n_doubles
+       enddo
-      l_a = doubles(i)
+       ! increment the alpha/beta  part for single excitations
-      ASSERT (l_a <= N_det)
+       call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
-
+       ! increment the beta /beta  part for single excitations
-      lrow = psi_bilinear_matrix_rows(l_a)
+       call off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
-      ASSERT (lrow <= N_det_alpha_unique)
+     enddo
-
+     
-      do l= 1, N_states
+     ! Compute Hij for all beta doubles
-       c_1(l) = u_t(l,l_a)
+     ! ----------------------------------
-       c_2(l) = u_t(l,k_a)
+     
-      enddo
+     do i=1,n_doubles
-      call off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
+       l_b = doubles(i)
-    enddo
+       ASSERT (l_b <= N_det)
-
+       
-
+       lcol = psi_bilinear_matrix_transp_columns(l_b)
-    ! Single and double beta excitations
+       ASSERT (lcol <= N_det_beta_unique)
-    ! ==================================
+       
-
+       l_a = psi_bilinear_matrix_transp_order(l_b)
-
+       do l= 1, N_states
-    ! Initial determinant is at k_a in alpha-major representation
+         c_1(l) = u_t(l,l_a)
-    ! -----------------------------------------------------------------------
+         c_2(l) = u_t(l,k_a)
-
+       enddo
-    krow = psi_bilinear_matrix_rows(k_a)
+       call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
-    kcol = psi_bilinear_matrix_columns(k_a)
+       ASSERT (l_a <= N_det)
-
+       
-    tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
+     enddo
-    tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+     
-
+     
-    spindet(1:$N_int) = tmp_det(1:$N_int,2)
+     ! Diagonal contribution
-
+     ! =====================
-    ! Initial determinant is at k_b in beta-major representation
+     
-    ! -----------------------------------------------------------------------
+     
-
+     ! Initial determinant is at k_a in alpha-major representation
-    k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
+     ! -----------------------------------------------------------------------
-    ASSERT (k_b <= N_det)
+     
-
+     krow = psi_bilinear_matrix_rows(k_a)
-    ! Loop inside the alpha row to gather all the connected betas
+     ASSERT (krow <= N_det_alpha_unique)
-    lrow = psi_bilinear_matrix_transp_rows(k_b)
+     
-    l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
+     kcol = psi_bilinear_matrix_columns(k_a)
-    do i=1,N_det_beta_unique
+     ASSERT (kcol <= N_det_beta_unique)
-      if (l_b > N_det) exit
+     
-      lrow = psi_bilinear_matrix_transp_rows(l_b)
+     tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-      if (lrow /= krow) exit
+     tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-      lcol = psi_bilinear_matrix_transp_columns(l_b)
+     
-      ASSERT (lcol <= N_det_beta_unique)
+     double precision, external     :: diag_wee_mat_elem, diag_S_mat_elem
-
+     
-      buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
+     double precision               :: c_1(N_states),c_2(N_states)
-      idx(i) = l_b
+     do l = 1, N_states
-      l_b = l_b+1
+       c_1(l) = u_t(l,k_a)
-    enddo
+     enddo
-    i = i-1
+     
-
+     call diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4)
-    call get_all_spin_singles_and_doubles_$N_int(                    &
+     
-        buffer, idx, spindet, i,                                     &
+   end do
-        singles_b, doubles, n_singles_b, n_doubles )
+   !!$OMP END DO
-
+   deallocate(buffer, singles_a, singles_b, doubles, idx)
-    ! Compute Hij for all beta singles
+   !!$OMP END PARALLEL
-    ! ----------------------------------
+   
-
-    tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-    do i=1,n_singles_b
-      l_b = singles_b(i)
-      ASSERT (l_b <= N_det)
-
-      lcol = psi_bilinear_matrix_transp_columns(l_b)
-      ASSERT (lcol <= N_det_beta_unique)
-
-      tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
-      l_a = psi_bilinear_matrix_transp_order(l_b)
-      do l= 1, N_states
-       c_1(l) = u_t(l,l_a)
-       c_2(l) = u_t(l,k_a)
-      enddo
-      ! increment the alpha/beta  part for single excitations 
-      call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
-      ! increment the beta /beta  part for single excitations 
-      call off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
-    enddo
-
-    ! Compute Hij for all beta doubles
-    ! ----------------------------------
-
-    do i=1,n_doubles
-      l_b = doubles(i)
-      ASSERT (l_b <= N_det)
-
-      lcol = psi_bilinear_matrix_transp_columns(l_b)
-      ASSERT (lcol <= N_det_beta_unique)
-
-      l_a = psi_bilinear_matrix_transp_order(l_b)
-      do l= 1, N_states
-       c_1(l) = u_t(l,l_a)
-       c_2(l) = u_t(l,k_a)
-      enddo
-      call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
-      ASSERT (l_a <= N_det)
-
-    enddo
-
-
-    ! Diagonal contribution
-    ! =====================
-
-
-    ! Initial determinant is at k_a in alpha-major representation
-    ! -----------------------------------------------------------------------
-
-    krow = psi_bilinear_matrix_rows(k_a)
-    ASSERT (krow <= N_det_alpha_unique)
-
-    kcol = psi_bilinear_matrix_columns(k_a)
-    ASSERT (kcol <= N_det_beta_unique)
-
-    tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-    tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-
-    double precision, external :: diag_wee_mat_elem, diag_S_mat_elem
-
-    double precision :: c_1(N_states),c_2(N_states)
-    do l = 1, N_states
-     c_1(l) = u_t(l,k_a)
-    enddo
-
-    call diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4)
-
-  end do
- !!$OMP END DO
-  deallocate(buffer, singles_a, singles_b, doubles, idx)
- !!$OMP END PARALLEL
-
 end
-
+ 
-SUBST [ N_int ]
+ SUBST [ N_int ]
-
+ 
-1;;
+ 1;;
-2;;
+ 2;;
-3;;
+ 3;;
-4;;
+ 4;;
-N_int;;
+ N_int;;
-
+ 
-END_TEMPLATE
+ END_TEMPLATE
-
+ 

src/two_body_rdm/two_rdm.irp.f

@@ -1,84 +1,62 @@
-
- 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
- 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,:)) & 
-                                         +          two_rdm_alpha_alpha_mo(i,j,k,l,:)                                  &
-                                         +          two_rdm_beta_beta_mo(i,j,k,l,:)         
-    enddo
-   enddo 
-  enddo
- enddo
-  
- END_PROVIDER 
-
-
  BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
- implicit none
+  implicit none
- BEGIN_DOC
+  BEGIN_DOC
- !  two_rdm_alpha_beta(i,j,k,l) = <Psi| a^{dagger}_{j,alpha} a^{dagger}_{l,beta} a_{k,beta} a_{i,alpha} | Psi>
+  !  two_rdm_alpha_beta(i,j,k,l) = <Psi| a^{dagger}_{j,alpha} a^{dagger}_{l,beta} a_{k,beta} a_{i,alpha} | Psi>
- !                     1 1 2 2  = chemist notations 
+  !                     1 1 2 2  = chemist notations
- !  note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
+  !  note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
- !  
+  !
- END_DOC
+  END_DOC
- integer :: dim1,dim2,dim3,dim4
+  integer                        :: dim1,dim2,dim3,dim4
- double precision :: cpu_0,cpu_1
+  double precision               :: cpu_0,cpu_1
- dim1 = mo_num
+  dim1 = mo_num
- dim2 = mo_num
+  dim2 = mo_num
- dim3 = mo_num
+  dim3 = mo_num
- dim4 = mo_num
+  dim4 = mo_num
- two_rdm_alpha_beta_mo = 0.d0
+  two_rdm_alpha_beta_mo = 0.d0
- two_rdm_alpha_alpha_mo= 0.d0
+  two_rdm_alpha_alpha_mo= 0.d0
- two_rdm_beta_beta_mo  = 0.d0
+  two_rdm_beta_beta_mo  = 0.d0
- print*,'providing two_rdm_alpha_beta ...'
+  print*,'providing two_rdm_alpha_beta ...'
- call wall_time(cpu_0)
+  call wall_time(cpu_0)
- call all_two_rdm_dm_nstates_openmp(two_rdm_alpha_alpha_mo,two_rdm_beta_beta_mo,two_rdm_alpha_beta_mo,dim1,dim2,dim3,dim4,psi_coef,size(psi_coef,2),size(psi_coef,1))
+  call all_two_rdm_dm_nstates_openmp(two_rdm_alpha_alpha_mo,two_rdm_beta_beta_mo,two_rdm_alpha_beta_mo,dim1,dim2,dim3,dim4,psi_coef,size(psi_coef,2),size(psi_coef,1))
- call wall_time(cpu_1)
+  call wall_time(cpu_1)
- print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)
+  print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)
-
+  
- END_PROVIDER 
+END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
- implicit none
+  implicit none
- BEGIN_DOC
+  BEGIN_DOC
- !  two_rdm_alpha_beta_mo_physicist,(i,j,k,l) = <Psi| a^{dagger}_{k,alpha} a^{dagger}_{l,beta} a_{j,beta} a_{i,alpha} | Psi>
+  !  two_rdm_alpha_beta_mo_physicist,(i,j,k,l) = <Psi| a^{dagger}_{k,alpha} a^{dagger}_{l,beta} a_{j,beta} a_{i,alpha} | Psi>
- !                                   1 2 1 2  = physicist notations 
+  !                                   1 2 1 2  = physicist notations
- !  note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
+  !  note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
- !  
+  !
- END_DOC
+  END_DOC
- integer :: i,j,k,l,istate
+  integer                        :: i,j,k,l,istate
- double precision :: cpu_0,cpu_1
+  double precision               :: cpu_0,cpu_1
- two_rdm_alpha_beta_mo_physicist = 0.d0
+  two_rdm_alpha_beta_mo_physicist = 0.d0
- print*,'providing two_rdm_alpha_beta_mo_physicist ...'
+  print*,'providing two_rdm_alpha_beta_mo_physicist ...'
- call wall_time(cpu_0)
+  call wall_time(cpu_0)
- do istate = 1, N_states 
+  do istate = 1, N_states
-  do i = 1, mo_num
+    do i = 1, mo_num
-   do j = 1, mo_num
+      do j = 1, mo_num
-    do k = 1, mo_num
+        do k = 1, mo_num
-     do l = 1, mo_num
+          do l = 1, mo_num
-      !                               1 2 1 2                                 1 1 2 2 
+            !                               1 2 1 2                                 1 1 2 2
-      two_rdm_alpha_beta_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_beta_mo(i,l,j,k,istate)
+            two_rdm_alpha_beta_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_beta_mo(i,l,j,k,istate)
-      two_rdm_alpha_alpha_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_alpha_mo(i,l,j,k,istate)
+            two_rdm_alpha_alpha_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_alpha_mo(i,l,j,k,istate)
-      two_rdm_beta_beta_mo_physicist(l,k,i,j,istate) = two_rdm_beta_beta_mo(i,l,j,k,istate)
+            two_rdm_beta_beta_mo_physicist(l,k,i,j,istate) = two_rdm_beta_beta_mo(i,l,j,k,istate)
-     enddo
+          enddo
+        enddo
+      enddo
     enddo
-   enddo
   enddo
- enddo
+  call wall_time(cpu_1)
- call wall_time(cpu_1)
+  print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
- print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
+  
-
+END_PROVIDER
- END_PROVIDER