Merge pull request #253 from QuantumPackage/dev-stable-rdm

modified the factor two in rdm

src/basis_correction/pbe_on_top.irp.f

@@ -41,11 +41,9 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top =  on_top_cas_mu_r(ipoint,istate)
-    on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
-    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
 !  We take the extrapolated on-top pair density 
@@ -105,11 +103,9 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top = on_top_cas_mu_r(ipoint,istate)
-    on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
-    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
 !  We take the extrapolated on-top pair density 
@@ -169,11 +165,9 @@
    
    if(mu_of_r_potential == "cas_ful")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
-    ! factor 2 because convention  N(N-1)/ 2 in provider on_top_cas_mu_r
+    on_top = on_top_cas_mu_r(ipoint,istate)
-    on_top = 1.d0 * on_top_cas_mu_r(ipoint,istate)
    else
     ! You take the on-top of the CAS wave function computed separately
-    ! No factor 2 because convention  N(N-1) in provider total_cas_on_top_density
     on_top = total_cas_on_top_density(ipoint,istate)
    endif
 !  We DO NOT take the extrapolated on-top pair density

src/mu_of_r/f_hf_utils.irp.f

@@ -86,6 +86,9 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
    enddo
   enddo
  enddo
+ ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
+ f_HF_val_ab *= 2.d0 
+ two_bod_dens *= 2.d0
 end
 
 
@@ -136,4 +139,6 @@ subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab)
    enddo
   enddo
  enddo
+ ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
+ int_f_HF_val_ab *= 2.d0
 end

src/mu_of_r/f_psi_i_a_v_utils.irp.f

@@ -49,6 +49,9 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
    enddo
   enddo
  enddo
+ ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
+ f_ii_val_ab  *= 2.d0 
+ two_bod_dens *= 2.d0
 end
 
 
@@ -142,6 +145,9 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
    f_ia_val_ab += v_tilde(i,a) * rho_tilde(i,a)
   enddo
  enddo
+ ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
+ f_ia_val_ab *= 2.d0 
+ two_bod_dens *= 2.d0
 end
 
 
@@ -194,7 +200,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
   do b = 1, n_act_orb  ! 2
    do c = 1, n_act_orb  ! 1
     do d = 1, n_act_orb  ! 2
-     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * 0.5d0 * act_2_rdm_ab_mo(d,c,b,a,istate)
+     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
      rho_tilde(b,a) += rho 
      two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) 
     enddo
@@ -222,6 +228,8 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
    f_aa_val_ab += v_tilde(b,a) * rho_tilde(b,a)
   enddo
  enddo
+
+ ! DO NOT multiply by two as in give_f_ii_val_ab and give_f_ia_val_ab because the N(N-1) normalization condition of the active two-rdm
 end
 
 BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)]