src/casscf/densities.irp.f

@@ -19,15 +19,14 @@ 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 ONLY IN THE RANGE OF ACTIVE MOS
+   ! the second-order density matrix in the basis of the starting MOs
-   ! The values are state averaged
+   ! matrices are state averaged
    !
-   ! We use the spin-free generators of mono-excitations
+   ! 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>
    !
-   ! P0tuvx(p,q,r,s) = chemist notation : 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
    END_DOC
    implicit none
    integer                        :: t,u,v,x
@@ -43,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 the 2 body RDM on the active part'
+    write(6,*) ' providing density matrix P0'
   endif
 
   P0tuvx= 0.d0
@@ -56,7 +55,11 @@ 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) = act_two_rdm_spin_trace_mo(t,v,u,x)
+       P0tuvx(t,u,v,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,10 +1,5 @@
 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
@@ -28,9 +23,18 @@ 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/README.rst

@@ -3,6 +3,6 @@ two_body_rdm
 ============
 
 Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
-arrays, with pysicists notation, consistent with the two-electron integrals in the
+maps, with pysicists notation, consistent with the two-electron integrals in the
 MO basis.
 

src/two_body_rdm/ab_only_routines.irp.f

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

src/two_body_rdm/orb_range_2_rdm.irp.f

@@ -4,10 +4,6 @@
  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
@@ -21,10 +17,6 @@
  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
@@ -38,10 +30,6 @@
  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
@@ -60,14 +48,6 @@
 
  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 IN CHEMIST NOTATIONS 
+! routine that update the DIAGONAL PART of the alpha/beta two body rdm 
  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 IN CHEMIST NOTATIONS 
+! routine that update the DIAGONAL PART of ALL THREE two body rdm 
  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  IN CHEMIST NOTATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS 
  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 IN CHEMIST NOTATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS 
  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 IN CHEMIST NOTATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS 
  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 IN CHEMIST NOTATIONS 
+! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for SINGLE EXCITATIONS 
  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 IN CHEMIST NOTATIONS 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS 
  END_DOC
  implicit none
  integer, intent(in) :: dim1,dim2,dim3,dim4

src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f

@@ -28,20 +28,7 @@
  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 the two body rdms in a specific range of orbitals for a given determinant det_1
+! routine that update the DIAGONAL PART of ALL THREE two body rdm 
-! 
-! 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
@@ -167,24 +154,7 @@
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS 
-!
-! 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
@@ -249,24 +219,7 @@
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS 
-!
-! 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
@@ -367,24 +320,7 @@
 
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS 
-!
-! 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
@@ -447,24 +383,7 @@
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for SINGLE EXCITATIONS 
-!
-! 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
@@ -530,24 +449,7 @@
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS 
-!
-! 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
@@ -603,24 +505,7 @@
  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 two body rdms in a specific range of orbitals for 
+! routine that update the OFF DIAGONAL PART of the beta /beta  2RDM only for DOUBLE EXCITATIONS 
-!
-! 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