diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f
index 0a5cfb49..2a595fe7 100644
--- a/src/casscf/get_energy.irp.f
+++ b/src/casscf/get_energy.irp.f
@@ -27,7 +27,7 @@ subroutine routine
     do ii = 1, n_act_orb
      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
+     accu(1) += state_av_act_two_rdm_spin_trace_mo(ii,jj,kk,ll) * integral
     enddo
    enddo
   enddo
diff --git a/src/two_body_rdm/all_states_2_rdm.irp.f b/src/two_body_rdm/all_states_2_rdm.irp.f
index b168da56..cd74758f 100644
--- a/src/two_body_rdm/all_states_2_rdm.irp.f
+++ b/src/two_body_rdm/all_states_2_rdm.irp.f
@@ -14,7 +14,7 @@
  ! condition for alpha/beta spin
  ispin = 1 
  all_states_act_two_rdm_alpha_alpha_mo = 0.D0
- call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  END_PROVIDER 
 
@@ -31,7 +31,7 @@
  ! condition for alpha/beta spin
  ispin = 2
  all_states_act_two_rdm_beta_beta_mo = 0.d0
- call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  END_PROVIDER 
 
@@ -53,7 +53,7 @@
  ispin = 3 
  print*,'ispin = ',ispin
  all_states_act_two_rdm_alpha_beta_mo = 0.d0
- call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  END_PROVIDER 
 
@@ -77,7 +77,7 @@
  all_states_act_two_rdm_spin_trace_mo = 0.d0
  integer :: i
 
- call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ call orb_range_all_states_two_rdm_openmp(all_states_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  END_PROVIDER 
 
diff --git a/src/two_body_rdm/routines_compute_2rdm_all_states.irp.f b/src/two_body_rdm/routines_compute_2rdm_all_states.irp.f
index 27b2dfe3..7606e353 100644
--- a/src/two_body_rdm/routines_compute_2rdm_all_states.irp.f
+++ b/src/two_body_rdm/routines_compute_2rdm_all_states.irp.f
@@ -5,7 +5,7 @@
 ! routine that update the DIAGONAL PART of the alpha/beta two body rdm in a specific range of orbitals 
  END_DOC
  implicit none
- integer, intent(in) :: dim1
+ integer, intent(in) :: dim1,N_st
  double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
  integer(bit_kind), intent(in)  :: det_1(N_int,2)
  integer(bit_kind), intent(in)  :: orb_bitmask(N_int)
@@ -123,8 +123,8 @@
      i2 = occ(j,2)
      h1 = list_orb_reverse(i1)
      h2 = list_orb_reverse(i2)
-     big_array(h1,h2,h1,h2,istate) += 0.5d0 * (c_1 )
-     big_array(h2,h1,h2,h1,istate) += 0.5d0 * (c_1 )
+     big_array(h1,h2,h1,h2,istate) += 0.5d0 * c_1(istate)
+     big_array(h2,h1,h2,h1,istate) += 0.5d0 * c_1(istate)
     enddo 
    enddo
    do i = 1, n_occ_ab(1)
@@ -334,6 +334,7 @@
    p1 = exc(1,2,2)
    if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
    p1 = list_orb_reverse(p1)
+   do istate = 1, N_st
     do i = 1, n_occ_ab(1)
      h2 = occ(i,1)
      if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
@@ -341,6 +342,7 @@
      big_array(h1,h2,p1,h2,istate) +=  0.5d0 * c_1(istate) * phase
      big_array(h2,h1,h2,p1,istate) +=  0.5d0 * c_1(istate) * phase
     enddo 
+   enddo
   endif
  endif
  end