added documentation for the two-rdm

src/casscf/densities.irp.f

@@ -42,7 +42,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
    integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
    
   if (bavard) then
-    write(6,*) ' providing density matrix P0'
+    write(6,*) ' providing the 2 body RDM on the active part'
   endif
 
   P0tuvx= 0.d0
@@ -55,11 +55,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
       uu = list_act(u)
       do t = 1, n_act_orb
        tt = list_act(t)
-       P0tuvx(t,u,v,x) =                                   &
+       P0tuvx(t,u,v,x) = act_two_rdm_spin_trace_mo(t,v,u,x)
-         state_average_weight(istate) *                    &
-         ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) +    &
-           two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) +    &
-           two_rdm_beta_beta_mo  (tt,uu,vv,xx,istate) )
       enddo
      enddo 
     enddo

src/casscf/get_energy.irp.f

@@ -1,5 +1,10 @@
 program print_2rdm
  implicit none
+ BEGIN_DOC
+ ! get the active part of the bielectronic energy on a given wave function.
+ !
+ ! useful to test the active part of the spin trace 2 rdms
+ END_DOC
  read_wf = .True.
  touch read_wf
  call routine
@@ -23,18 +28,9 @@ subroutine routine
      i = list_act(ii)
      integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
      accu(1) += act_two_rdm_spin_trace_mo(ii,jj,kk,ll) * integral
-    !if(dabs(act_two_rdm_spin_trace_mo(ii,jj,kk,ll)).gt.thr)then
-    !print*,'',ii,jj,kk,ll,act_two_rdm_spin_trace_mo(ii,jj,kk,ll)*integral
-    !print*,'accu',accu(1)
-    !endif
     enddo
    enddo
   enddo
  enddo
  print*,'accu             = ',accu(1)
- print*,'psi_energy_two_e = ',psi_energy_two_e
-!double precision :: hij
-!call i_H_j_double_alpha_beta(psi_det(1,1,1),psi_det(1,1,2),N_int,hij)
-!print*,'hij * 2',hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0
-!print*,'psi diag         = ',psi_energy_two_e - hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0
 end

src/two_body_rdm/ab_only_routines.irp.f

@@ -1,9 +1,9 @@
 
- subroutine two_rdm_dm_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze)
+ subroutine two_rdm_ab_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze)
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
+  ! Computes the alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS 
   !
   ! Assumes that the determinants are in psi_det
   !
@@ -27,7 +27,7 @@
       size(u_t, 1),                                                  &
       N_det, N_st)
 
-  call two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
+  call two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
   deallocate(u_t)
 
   do k=1,N_st
@@ -37,11 +37,11 @@
  end
 
 
- subroutine two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+ subroutine two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
+  ! Computes the alpha/beta part of the two-body density matrix
   !
   ! Default should be 1,N_det,0,1
   END_DOC
@@ -55,20 +55,20 @@
 
   select case (N_int)
     case (1)
-      call two_rdm_dm_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+      call two_rdm_ab_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
     case (2)
-      call two_rdm_dm_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+      call two_rdm_ab_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
     case (3)
-      call two_rdm_dm_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+      call two_rdm_ab_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
     case (4)
-      call two_rdm_dm_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+      call two_rdm_ab_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
     case default
-      call two_rdm_dm_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+      call two_rdm_ab_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
   end select
  end
  BEGIN_TEMPLATE
 
- subroutine two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+ subroutine two_rdm_ab_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
   use bitmasks
   implicit none
   integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep

src/two_body_rdm/all_2rdm_routines.irp.f

@@ -2,7 +2,7 @@ subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,
    use bitmasks
    implicit none
    BEGIN_DOC
-   ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
+   ! Computes the alpha/alpha, beta/beta and alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS 
    !
    ! Assumes that the determinants are in psi_det
    !

src/two_body_rdm/orb_range_2_rdm.irp.f

@@ -4,6 +4,10 @@
  BEGIN_PROVIDER [double precision, act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
  implicit none
  double precision, allocatable :: state_weights(:) 
+ BEGIN_DOC
+! act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs
+!                                     = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
+ END_DOC 
  allocate(state_weights(N_states))
  state_weights = 1.d0/dble(N_states)
  integer :: ispin
@@ -17,6 +21,10 @@
  BEGIN_PROVIDER [double precision, act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
  implicit none
  double precision, allocatable :: state_weights(:) 
+ BEGIN_DOC
+! act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs
+!                                     = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
+ END_DOC 
  allocate(state_weights(N_states))
  state_weights = 1.d0/dble(N_states)
  integer :: ispin
@@ -30,6 +38,10 @@
  BEGIN_PROVIDER [double precision, act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
  implicit none
  double precision, allocatable :: state_weights(:) 
+ BEGIN_DOC
+! act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
+!                                     = <Psi| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi>
+ END_DOC 
  allocate(state_weights(N_states))
  state_weights = 1.d0/dble(N_states)
  integer :: ispin
@@ -48,6 +60,14 @@
 
  BEGIN_PROVIDER [double precision, act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
  implicit none
+ BEGIN_DOC
+! act_two_rdm_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices
+! The active part of the two-electron energy can be computed as:
+!
+!   \sum_{i,j,k,l = 1, n_act_orb} act_two_rdm_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > 
+!
+!   with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l)
+ END_DOC
  double precision, allocatable :: state_weights(:) 
  allocate(state_weights(N_states))
  state_weights = 1.d0/dble(N_states)

src/two_body_rdm/routines_compute_2rdm.irp.f

@@ -3,7 +3,7 @@
  subroutine diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the DIAGONAL PART of the alpha/beta two body rdm 
+! routine that update the DIAGONAL PART of the alpha/beta two body rdm IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4
@@ -31,7 +31,7 @@
  subroutine diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the DIAGONAL PART of ALL THREE two body rdm 
+! routine that update the DIAGONAL PART of ALL THREE two body rdm IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4
@@ -77,7 +77,7 @@
  subroutine off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS  IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4
@@ -101,7 +101,7 @@
  subroutine off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4
@@ -140,7 +140,7 @@
 
  subroutine off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS 
  END_DOC
  use bitmasks
  implicit none
@@ -177,7 +177,7 @@
  subroutine off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4
@@ -214,7 +214,7 @@
  subroutine off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4

src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f

@@ -28,7 +28,20 @@
  subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the DIAGONAL PART of ALL THREE two body rdm 
+! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,ispin
@@ -154,7 +167,24 @@
  subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 3 or 4 will do something
  END_DOC
  implicit none
  integer, intent(in) :: dim1,ispin
@@ -219,7 +249,24 @@
  subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 3 or 4 will do something
  END_DOC
  implicit none
  integer, intent(in) :: dim1,ispin
@@ -320,7 +367,24 @@
 
  subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 1 or 4 will do something
  END_DOC
  use bitmasks
  implicit none
@@ -383,7 +447,24 @@
  subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for SINGLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a BETA  SINGLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 2 or 4 will do something
  END_DOC
  implicit none
  integer, intent(in) :: dim1,ispin
@@ -449,7 +530,24 @@
  subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 1 or 4 will do something
  END_DOC
  implicit none
  integer, intent(in) :: dim1,ispin
@@ -505,7 +603,24 @@
  subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin)
  use bitmasks
  BEGIN_DOC
-! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for DOUBLE EXCITATIONS 
+! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for 
+!
+! a given couple of determinant det_1, det_2 being a BETA /BETA  DOUBLE excitation with respect to one another
+! 
+! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
+! 
+! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation 
+!
+! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
+!
+! ispin determines which spin-spin component of the two-rdm you will update 
+!
+! ispin   == 1 :: alpha/ alpha
+! ispin   == 2 :: beta / beta
+! ispin   == 3 :: alpha/ beta
+! ispin   == 4 :: spin traced <=> total two-rdm 
+!
+! here, only ispin == 2 or 4 will do something
  END_DOC
  implicit none