print conf with large coef in WF

src/info_prop.F90

@@ -35,7 +35,10 @@ Module Info_prop
   type, public :: prop_infotype
      !> @brief Number of properties
      Integer(KIND=kd_int) :: nprop
-     character*5, dimension(nprop_max) :: whichprop 
+     !> @brief Names of desired properties
+     character*5, dimension(nprop_max) :: whichprop
+     !> @brief WF printing threshold, default 0.05
+     Real (KIND=kd_dble) :: seuilcoef
   end type prop_infotype
 
 end Module Info_prop

src/lect_data_prop.F90

@@ -92,8 +92,9 @@ subroutine read_propinp(prop_info, nref0lu, iunit)
   character*5, dimension(:), allocatable :: whichprop
   integer (KIND=kd_int) :: stot, sz, vec_irrep, nvec, &
        nref0, nelact, itmp
+  Real (KIND=kd_dble) :: seuilcoef
   
-  namelist /propinp/ nprop, whichprop
+  namelist /propinp/ nprop, whichprop, seuilcoef
   namelist /vecinp/ stot, sz, vec_irrep, nvec, nref0, nelact
 
   ! vecinp
@@ -113,11 +114,12 @@ subroutine read_propinp(prop_info, nref0lu, iunit)
   nprop = 0
   allocate(whichprop(nprop_max))
   whichprop(1:nprop_max) = "xxxxx"
+  seuilcoef = 0.05d0
   read(iunit, propinp)
 
   prop_info%nprop = nprop
   prop_info%whichprop(1:nprop) =  whichprop(1:nprop)
-
+  prop_info%seuilcoef = seuilcoef
 
 end subroutine read_propinp
 

src/prop.F90

@@ -88,7 +88,7 @@ Program proprietes
 
   ! autres
   integer ::  idet, ivec, iprop, nb_thread
-  integer ::  nr, nq, iwr, i, j, nref0lu
+  integer ::  nr, nq, iwr, i, j, n, nref0lu
   real(kd_dble) :: tstart, tend, wstart, wend, wstart1, wend1,t1, t2, wt1, wt2
 
   CHARACTER(LEN=8)     :: date      ! returned values from DATE_AND_TIME()
@@ -206,9 +206,16 @@ Program proprietes
 !!$-----
 !!$ --- Projection sur le cas
 !!$-----
+  write(f_output,*)
+  write(f_output,*)"================================================================================"
+  write(f_output,*)
+  write(f_output,*)"                            Properties section                                  "
+  write(f_output,*)
+  write(f_output,*)"================================================================================"
   call gettime(t1,wt1)
   write(f_output,*)
-  write(f_output,*) " >>> Vectors projection on active space"
+  n=len("Vectors projection on active space")
+  call wrt_mess(f_output,n,"Vectors projection on active space")
   call cass_proj(psi, ener, ndet, nvec, nelact,det, det_info, o_info, &
        nblock, shtblkdet,nblkdet,deter_index,f_output)
   call gettime(t2,wt2)
@@ -230,9 +237,12 @@ Program proprietes
      do iprop = 1, nprop
         call lowercase(whichprop(iprop))
         select case (whichprop(iprop))
-        case("s-s2")
-           write(f_output,*) 
-           write(f_output,*) " >>>  Total S "
+           
+        case("s-s2") !------- Compute total S -------------------------------------
+           write(f_output,*)
+           write(f_output,*)
+           n=len("Total S")
+           call wrt_mess(f_output,n,"Total S")
            allocate(psi_S(nvec))
            psi_S(:) = 0.d0
            call s2(psi,psi_S,ndet,nvec,nelact, det, det_info, d, o_info, f_output)
@@ -253,21 +263,33 @@ Program proprietes
            write(f_output,*)
            deallocate(psi_S)
 9012       format(10(F18.8,2x))
+           
         case("wf") !------- Print WFs ---------------------------------------------
-           write(f_output,*) 
-           write(f_output,*) " >>>  Total WFs"
+           write(f_output,*)
+           write(f_output,*)
+           n=len("Total WFs")
+           call wrt_mess(f_output,n,"Total WFs")
            call wrt_WF(ndet,nvec,1,nvec,det,psi,o_info,f_output)
+           
         case("pref1") !---- Project WF on Ref1 ------------------------------------
-           write(f_output,*) 
-           write(f_output,*) " >>>  Projection on Ref1"
+           write(f_output,*)
+           write(f_output,*)
+           n=len("Projection on Ref1")
+           call wrt_mess(f_output,n,"Projection on Ref1")
            write(f_output,*) "      Not yet implemented"
+           
         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"
+           write(f_output,*)
+           write(f_output,*)
+           n=len("Projection on determinants with coefficients largest than a threshold")
+           call wrt_mess(f_output,n,"Projection on determinants with coefficients largest than a threshold")
+           write(f_output,'("  Threshold = ",f8.5)')  prop_info%seuilcoef
+           call largecoef(prop_info%seuilcoef,ndet,nvec,det,psi,o_info,f_output)
+           
         case("dens") !----- 1e-density matrix --------------------------------------
-           write(f_output,*) 
-           write(f_output,*) ">>> One particule density matrix"
+           write(f_output,*)
+           n=len("One particule density matrices")
+           call wrt_mess(f_output,n,"One particule density matrices")
 #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)
@@ -283,7 +305,11 @@ Program proprietes
   end if
 
   write(f_output,*)
-  write(f_output,*) ' <<< End of prop code >>> '
+  write(f_output,*)"================================================================================"
+  write(f_output,*)
+  write(f_output,*)"                        End of properties section                               "
+  write(f_output,*)
+  write(f_output,*)"================================================================================"
   write(f_output,*)
 
   deallocate(shtblkdet, nblkdet, deter_index)
@@ -306,4 +332,27 @@ Program proprietes
   call spindetact_all_free(rspin)
   call deter_blocklist_free(d)
 
+contains
+  subroutine wrt_mess(iunit,n,mess)
+    implicit none
+    character(len=*), intent(in) :: mess
+    Integer(KIND=kd_int), intent(in) :: iunit,n 
+    Integer :: i,d
+    
+    d= (80-n-2)/2
+    do i=1,d
+       write(iunit,9001,advance="no") ">"
+    end do
+    write(iunit,9001,advance="no") " "
+    do i=1,n
+       write(iunit,9001,advance="no") mess(i:i)
+    end do
+    write(iunit,9001,advance="no") " "
+    do i=1, 80-n-2-d
+       write(iunit,9001,advance="no") "<"
+    end do
+    write(iunit,*)
+9001 format(a1)
+    
+  end subroutine wrt_mess
 end Program proprietes

src/proprietes.F90

@@ -373,6 +373,220 @@ subroutine s2(psi,psi_S,ndet,nvec, nelact, det, det_info, d, o_info, iunit)
   end do
   !if (ldebug)
   ! write(f_output,9002) psi_S2(1:nvec)
-9002 format(" S2   ",10(2x,F18.8))
+  !9002 format(" S2   ",10(2x,F18.8))
  
 end subroutine s2
+
+
+!$====================================================================
+!> @brief Print WF contributions larger than seuilcoef
+!> @author MBL
+!> @date April 2022
+!
+!> @param[in] vect(ndet,nvec) = states matrix
+!> @param[in] det(ndet) = determinants
+!> @param[in] ndet = number of determinants
+!> @param[in] nvec = number of vectors
+!> @param[in] nvec = seuilcoef printing threshold 
+!> @param[in] iunit : output file
+!$====================================================================
+  Subroutine  largecoef(seuilcoef,ndet,nvec,det,vect,o_info,iunit)
+    ! -------- Donnes globales -----------------------------------
+    use dimensions
+    use info
+    use files
+    use utils_char
+    use typedet
+    use spindetact
+    Implicit none
+    ! -------- Donnes transmises -----------------------------------
+    Integer (KIND=kd_int), intent(in) :: ndet,nvec, iunit
+    Real (KIND=kd_dble),   intent(in) :: seuilcoef
+    Real (KIND=kd_dble), dimension(ndet,nvec), intent(in) :: vect
+    type(deter), dimension(ndet), intent(in) :: det
+    type(o_infotype), intent(in)             :: o_info
+    
+    !   -------- Donnes locales -----------------------------------
+    Integer, parameter :: pas = 5
+    Integer (KIND=kd_int) :: npas, kpas
+    Integer (KIND=kd_int) :: shtrou, shpart, shtu, shtd, shpu, shpd
+    Integer (KIND=kd_int) :: idet, ivec, ibit, j
+    Logical :: OK
+    Character*1, dimension(2) :: spin
+    Integer(KIND=kd_int), dimension(2) :: sh
+
+    !============================================================
+    ! -------- Code ---------------------------------------------
+    !-----
+    !----- Initialisations 
+    !-----
+    npas = nvec/pas
+
+    shtrou  = 2*o_info%ngel +  o_info%nocc + o_info%nligo
+    shtu    = o_info%ngel
+    shtd    = o_info%ngel +  o_info%nocc + o_info%nligo
+    shpart  = 2*(o_info%ngel +  o_info%nocc + o_info%nligo + o_info%nact) + o_info%nligv + o_info%nvirt
+    shpu    = o_info%ngel + o_info%nocc + o_info%nligo + o_info%nact
+    shpd    = o_info%ngel + o_info%nocc + o_info%nligo + o_info%nact + o_info%nligv + o_info%nvirt
+
+    if (npas.eq.0) goto 10
+    kp : do kpas = 1,npas
+       do idet = 1, ndet
+          ! test seuil
+          OK = .false.
+          do ivec =1,nvec
+             if (abs(vect(idet,ivec)).gt.seuilcoef) OK=.true.
+          end do
+          if (.not.OK) cycle
+          ! OK = .true. => write det
+          write(iunit,'(i9,2x)',advance='no') idet
+          ! actives
+          do ibit = 0,o_info%nact-1
+             write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+          end do
+          write(iunit,'(" ")',advance='no')
+          do ibit = o_info%nact, 2*o_info%nact-1
+             write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+          end do
+          write(iunit,'(" ")',advance='no')
+          ! trous
+          if (det(idet)%dettr(1).eq.0) then ! 0 trou
+             write(iunit,'(11x," -> ")',advance='no')
+          else if (det(idet)%dettr(2).eq.0) then ! 1 trou
+             if (det(idet)%dettr(1).le.shtrou) then
+                spin(1)="u"
+                sh(1)  = shtu
+             else
+                spin(1)="d"
+                sh(1)  = shtd
+             end if
+             write(iunit,'(I4,a1,6x," -> ")',advance='no') det(idet)%dettr(1)-sh(1), spin(1)
+          else ! 2 trou
+             if (det(idet)%dettr(1).le.shtrou) then
+                spin(:)="u"
+                sh(:)  = shtu
+             else if (det(idet)%dettr(2).le.shtrou) then
+                spin(1)="d"
+                sh(1)  = shtd
+                spin(2)="u"
+                sh(2)  = shtu
+             else
+                spin(:)="d"
+                sh(:)  = shtd
+             end if
+             write(iunit,'(I4,a1,I4,a1,1x," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
+          end if
+          ! particules
+          if (det(idet)%detprt(1).eq.0) then ! 0 particules
+             write(iunit,'(11x,2x)',advance='no')
+          else if (det(idet)%detprt(2).eq.0) then ! 1 particules
+             if (det(idet)%detprt(1).le.shpart) then
+                spin(1)="u"
+                sh(1)  = shpu
+             else
+                spin(1)="d"
+                sh(1)  = shpd
+             end if
+             write(iunit,'(I4,a1,6x,2x)',advance='no') det(idet)%detprt(1)-sh(1), spin(1)
+          else ! 2 particules
+             if (det(idet)%detprt(2).le.shpart) then
+                spin(:)="u"
+                sh(:)  = shpu
+             else if (det(idet)%detprt(1).le.shpart) then
+                spin(1)="u"
+                sh(1)  = shpu
+                spin(2)="d"
+                sh(2)  = shpd
+             else
+                spin(:)="d"
+                sh(:)  = shpd
+             end if
+             write(iunit,'(I4,a1,I4,a1,3x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
+          end if
+          ! vecteur
+          write(iunit,9990)  (vect(idet,ivec), ivec = 1 + (kpas-1)*pas, 1 + kpas*pas)
+       end do
+       write(iunit,*)
+    end do kp
+    
+10  if (mod(nvec,pas) .eq. 0) return
+    do idet = 1, ndet
+       ! test seuil
+       OK = .false.
+       do ivec =1,nvec
+          if (abs(vect(idet,ivec)).gt.seuilcoef) OK=.true.
+       end do
+       if (.not.OK) cycle
+       ! OK = .true. => write det
+       write(iunit,'(i9,2x)',advance='no') idet
+       ! actives
+       do ibit = 0,o_info%nact-1
+          write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+       end do
+       write(iunit,'(" ")',advance='no')
+       do ibit = o_info%nact, 2*o_info%nact-1
+          write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+       end do
+       write(iunit,'(" ")',advance='no')
+       ! trous
+       if (det(idet)%dettr(1).eq.0) then ! 0 trou
+          write(iunit,'(11x," -> ")',advance='no')
+       else if (det(idet)%dettr(2).eq.0) then ! 1 trou
+          if (det(idet)%dettr(1).le.shtrou) then
+             spin(1)="u"
+             sh(1)  = shtu
+          else
+             spin(1)="d"
+             sh(1)  = shtd
+          end if
+          write(iunit,'(I4,a1,6x," -> ")',advance='no') det(idet)%dettr(1)-sh(1), spin(1)
+       else ! 2 trou
+          if (det(idet)%dettr(1).le.shtrou) then
+             spin(:)="u"
+             sh(:)  = shtu
+          else if (det(idet)%dettr(2).le.shtrou) then
+             spin(1)="d"
+             sh(1)  = shtd
+             spin(2)="u"
+             sh(2)  = shtu
+          else
+             spin(:)="d"
+             sh(:)  = shtd
+          end if
+          write(iunit,'(I4,a1,I4,a1,1x," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
+       end if
+       ! particules
+       if (det(idet)%detprt(1).eq.0) then ! 0 particules
+          write(iunit,'(11x,2x)',advance='no')
+       else if (det(idet)%detprt(2).eq.0) then ! 1 particules
+          if (det(idet)%detprt(1).le.shpart) then
+             spin(1)="u"
+             sh(1)  = shpu
+          else
+             spin(1)="d"
+             sh(1)  = shpd
+          end if
+          write(iunit,'(I4,a1,6x,2x)',advance='no') det(idet)%detprt(1)-sh(1), spin(1)
+       else ! 2 particules
+          if (det(idet)%detprt(2).le.shpart) then
+             spin(:)="u"
+             sh(:)  = shpu
+          else if (det(idet)%detprt(1).le.shpart) then
+             spin(1)="u"
+             sh(1)  = shpu
+             spin(2)="d"
+             sh(2)  = shpd
+          else
+             spin(:)="d"
+             sh(:)  = shpd
+          end if
+          write(iunit,'(I4,a1,I4,a1,3x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
+       end if
+       ! vecteur
+       write(iunit,9990)  (vect(idet,ivec), ivec = 1 + npas*pas, nvec )
+    end do
+    write(iunit,*)
+    
+9990 format(5(2x,F20.16))
+
+  end Subroutine largecoef

src/utils_wrt.F90

@@ -495,31 +495,32 @@ contains
     kp : do kpas = 1,npas
        do idet = 1, ndet
           write(iunit,'(i9,2x)',advance='no') idet
-       ! actives
-       do ibit = 0,o_info%nact-1
-          write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
-       end do
-       write(iunit,'(" ")',advance='no')
-       do ibit = o_info%nact, 2*o_info%nact-1
-          write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
-       end do
-       write(iunit,'(" ")',advance='no')
-       ! trous
-       if (det(idet)%dettr(1).eq.0) then ! 0 trou
-          write(iunit,'(11x," -> ")',advance='no')
-       else if (det(idet)%dettr(2).eq.0) then ! 1 trou
-          if (det(idet)%dettr(1).le.shtrou) then
-             spin(1)="u"
-             sh(1)  = shtu
-          else
-             spin(1)="d"
-             sh(1)  = shtd
-          end if
-          write(iunit,'(I4,a1,6x," -> ")',advance='no') det(idet)%dettr(1)-sh(1), spin(1)
-       else ! 2 trou
-          if (det(idet)%dettr(1).le.shtrou) then
-             spin(:)="u"
-             sh(:)  = shtu
+          ! actives
+          do ibit = 0,o_info%nact-1
+             write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+          end do
+          write(iunit,'(" ")',advance='no')
+          do ibit = o_info%nact, 2*o_info%nact-1
+             write(iunit,'(i1)',advance='no') ibits(det(idet)%detact, ibit, 1)
+          end do
+          write(iunit,'(" ")',advance='no')
+          ! trous
+          if (det(idet)%dettr(1).eq.0) then ! 0 trou
+             write(iunit,'(11x," -> ")',advance='no')
+          else if (det(idet)%dettr(2).eq.0) then ! 1 trou
+             if (det(idet)%dettr(1).le.shtrou) then
+                spin(1)="u"
+                sh(1)  = shtu
+             else
+                spin(1)="d"
+                sh(1)  = shtd
+             end if
+             write(iunit,'(I4,a1,6x," -> ")',advance='no') det(idet)%dettr(1)-sh(1), spin(1)
+          else ! 2 trou
+             if (det(idet)%dettr(1).le.shtrou) then
+                spin(:)="u"
+                sh(:)  = shtu
+                
           else if (det(idet)%dettr(2).le.shtrou) then
              spin(1)="d"
              sh(1)  = shtd
@@ -529,7 +530,7 @@ contains
              spin(:)="d"
              sh(:)  = shtd
           end if
-          write(iunit,'(I4,a1,I4,a1," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
+          write(iunit,'(I4,a1,I4,a1,1x," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
        end if
        ! particules
        if (det(idet)%detprt(1).eq.0) then ! 0 particules
@@ -556,7 +557,7 @@ contains
              spin(:)="d"
              sh(:)  = shpd
           end if
-          write(iunit,'(I4,a1,I4,a1,2x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
+          write(iunit,'(I4,a1,I4,a1,3x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
        end if
        ! vecteur
        write(iunit,9990)  (vect(idet,ivec), ivec = first_vec + (kpas-1)*pas, first_vec-1 + kpas*pas)
@@ -601,7 +602,7 @@ contains
              spin(:)="d"
              sh(:)  = shtd
           end if
-          write(iunit,'(I4,a1,I4,a1," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
+          write(iunit,'(I4,a1,I4,a1,1x," -> ")',advance='no') (det(idet)%dettr(j)-sh(j), spin(j), j=2,1,-1)
        end if
        ! particules
        if (det(idet)%detprt(1).eq.0) then ! 0 particules
@@ -628,7 +629,7 @@ contains
              spin(:)="d"
              sh(:)  = shpd
           end if
-          write(iunit,'(I4,a1,I4,a1,2x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
+          write(iunit,'(I4,a1,I4,a1,3x)',advance='no') (det(idet)%detprt(j)-sh(j), spin(j), j=1,2)
        end if
        ! vecteur
        write(iunit,9990)  (vect(idet,ivec), ivec = first_vec + npas*pas, last_vec )