termes diagonaux semblent OK - prog verif

cmake/RelaxSE_dev.cmake

@@ -6,4 +6,5 @@ set(RELAXSE_dev
   src/generator/dev/dens_exdiag_tp.F90
   src/generator/dev/dens_exdiag_ap.F90
   src/generator/dev/dens_exdiag_ta.F90
+  src/generator/dev/verif_dens.F90
   )

src/gencode/aaaa_000_000_generated.F90

@@ -147,7 +147,7 @@ contains
                 end do
                 deallocate(andIJ)
                 deallocate(spinandIJ)
-                elm = sign_act * psign * (1 - 2*modulo((Nh_I*Ndiffact_I + Nh_J*Ndiffact_J),2)) * &
+                elm = sign_act * psign * (1 - 2*modulo((Nh_I**Ndiffact_I + Nh_J*Ndiffact_J),2)) * &
                      tsign * JK
 
                   do m = 1, nvec

src/prop.F90

@@ -45,6 +45,7 @@ Program proprietes
   use utils_wrt
 #ifdef VAR_DEV
   use densite
+!  use density_matrix
 #endif
   
 !!$ -------- Donnes locales -----------------------------------
@@ -78,7 +79,7 @@ Program proprietes
   character*5, dimension(:), allocatable :: whichprop
 
   ! proprietes
-  real(kd_dble), dimension(:,:,:),  allocatable :: rho_tot,rho_spin
+  !  real(kd_dble), dimension(:,:,:),  allocatable :: rho_tot,rho_spin
   
   ! impressions
   integer, parameter :: pas = 10
@@ -92,9 +93,13 @@ Program proprietes
   CHARACTER(LEN=8)     :: date      ! returned values from DATE_AND_TIME()
   CHARACTER(LEN=10)    :: time
   CHARACTER(LEN=5)     :: zone
-  INTEGER,DIMENSION(8) :: values 
-  
-  
+  INTEGER,DIMENSION(8) :: values
+
+#ifdef VAR_DEV
+  integer :: iorb1, iorb2, norb
+  real(kind=kd_dble), allocatable :: rho(:, :)
+  real(kind=kd_dble), allocatable :: WF(:)
+#endif
   
 !!$============================================================
 !!$ -------- Code ---------------------------------------------
@@ -211,9 +216,9 @@ Program proprietes
        'CPU time ',t2-t1,'s, Wall time ',wt2-wt1,' s'
   flush(f_output)
   
-!!$-----
+!!$-----------------------------------------------------------------------------
 !!$ --- Autres proprits
-!!$-----
+!!$-----------------------------------------------------------------------------
   nprop = prop_info%nprop
   if (nprop.eq.0) goto 9999
   
@@ -247,26 +252,41 @@ Program proprietes
         write(f_output,*)
         deallocate(psi_S)
 9012    format(10(F18.8,2x))
-     case("wf")
+     case("wf") !------- Print WFs ---------------------------------------------
         write(f_output,*) 
         write(f_output,*) " >>>  Total WFs"
         call wrt_WF(ndet,nvec,1,nvec,det,psi,o_info,f_output)
-     case("pref1")
-        ! Projection sur les ref1
+     case("pref1") !---- Project WF on Ref1 ------------------------------------
         write(f_output,*) 
         write(f_output,*) " >>>  Projection on Ref1"
         write(f_output,*) "      Not yet implemented"
-     case("lcoef")
-        ! Projection sur les ref1
+     case("lcoef") !---- Print largest determinants ----------------------------
         write(f_output,*) 
         write(f_output,*) " >>>  Projection on determinants with coefficients largest than 0.05"
         write(f_output,*) "      Not yet implemented"
-     case("dens")
-        ! Calcul de la matrice densite
+     case("dens") !----- 1e-density matrix --------------------------------------
         write(f_output,*) 
         write(f_output,*) ">>> One particule density matrix"
 #ifdef VAR_DEV
         call  dens(ndet,nvec, psi,det, d,rspin, o_info,g_info,prog_info)
+        call  verif_dens(ndet,nvec, psi,det, d,rspin, o_info,g_info,prog_info)
+!         allocate(WF(ndet))
+!         Wf(:) = psi(:,2)
+!         norb = o_info%ntot + o_info%ngel + o_info%ndel
+!         allocate(rho(norb,norb))
+!         rho(:,:) = 0.d0
+!         rho = calc_density_matrix(WF, o_info, det)
+        
+!         write(f_output,*) "    calc_density_matrix : non zero terms"
+!         do iorb2=1,norb
+!            !           do iorb1=1,norb
+!            iorb1 = iorb2
+!               if (abs(rho(iorb1,iorb2)).gt.1.d-10) write(f_output,9983) iorb1,iorb2, rho(iorb1,iorb2)
+! !           end do
+!         end do
+        write(f_output,*)
+        write(f_output,*)
+9983    format (5x,2(I4,1x),2x,10(F22.16,1x))
 #endif
      end select
   end do

src/spindetact.F90

@@ -441,7 +441,7 @@ contains
     integer,intent(in)                  :: iout, nact
 
     integer :: i
-    logical :: debug = .true.
+    logical :: debug = .false.
     integer :: nCM1, nCM2
     integer, pointer :: detCM1(:), detCM2(:)
 

src/utils_wrt.F90

@@ -490,7 +490,6 @@ contains
 
     ! 10  if (mod(last_vec - first_vec + 1,pas) .eq. 0) return
 
-    write(iunit,*) "***** numeros particules faux"
     
     if (npas.eq.0) goto 10
     kp : do kpas = 1,npas