Possibility to increase the number of vectors to be computed

• possibility to change nvec after a restart to ask for more vectors
• possibility to specify the number of converged vectors to increase the convergence threshold on some vectors

src/RelaxSE_diag.F90

@@ -39,6 +39,7 @@ module SASS_diag
   use utils_wrt
   use utils_ortho
   use utils_intcase
+  use utils_char
   use typetargetvec
   use compute_hv
   use sort
@@ -117,6 +118,7 @@ contains
     
     integer :: ivec, jvec, kvec, iter, ish, jsh, i, j
     integer :: ncol, icol, jcol, idet
+    Integer :: n
 
     integer :: nguess
     integer :: iconv, nbre_col
@@ -134,7 +136,15 @@ contains
     real(kd_dble), dimension(:,:), allocatable :: Vm_tmp, Wm_tmp
 
     debugdav = bdav_info%debug_dav
-    open(666)
+    if (debugdav) then
+       if (prog_info%id_cpu.eq.0) then       
+          write(f_output,*)
+          write(f_output,*)">>>>> Debug_dav on <<<<<<<"
+          write(f_output,*)
+       end if
+       call noblancs(prog_info%prefix,n)
+       open(fdav,file=prog_info%prefix(1:n)//".debug",form="formatted")
+    end if
     open(667)
     
     ngel     = o_info%ngel
@@ -145,7 +155,7 @@ contains
     nvirt    = o_info%nvirt
     nd00     = det%detblock(1)%p%ndet
     sz       = v_info%sz
-    nvectot  = v_info%nvec
+    nvectot  = v_info%nvec !total nb of vectors to be computed (from INPUT)
     Sizeheff = bdav_info%Sizeheff
     NitDavid = bdav_info%NitDavid
     tol_orth = bdav_info%tol_orth 
@@ -195,23 +205,20 @@ contains
 
     if (prog_info%id_cpu.eq.0) then
        write(f_output,*)
-       write(f_output,*)&
-            "================================================================="
-       write(f_output,*) ">>> Davidson Diagonalisation with low memory"
+       write(f_output,*) ">>> Davidson Diagonalisation" ! with low memory"
        write(f_output,'(A,I0,A,ES8.1)') "   Search for ",nvectot,&
             " vectors with a convergence threshold of ",sqrt(tol_conv)
        write(f_output,*)
        call gettime(t1,ostart)
        write(f_output,'(A,F12.4,A)') 'Total Elapsed tWall ', ostart-prog_info%wstart,'s'
-
        write(f_output, '(A,I0,A)') 'Particles are divided in batches of ', sizebatch
+       write(f_output,*)&
+            "================================================================="
     endif
 
     if (debugdav) then 
        write(fdav,*) "Guess vectors"
-       do idet = 1, nd00
-          write(fdav,*) (psi_0_guess(idet,ivec), ivec = 1,nvectot)
-       enddo
+       call wrt_vect(fdav,nd00,ndet,nvectot,psi_0_guess)
        write(fdav,*)
     end if
 
@@ -226,14 +233,15 @@ contains
     ! aux nvectot vecteurs d'essai donnés  sur les nd00 premiers déterminants
     if (prog_info%restart) then
        nguess = ndet
+       iter0 = iter0-1
     else
        nguess = nd00
     endif
     Vm(1:nguess,1:nvectot,1)     = psi_0_guess(1:nguess,1:nvectot)
     psi_SASS(1:nguess,1:nvectot) = psi_0_guess(1:nguess,1:nvectot) 
 
-!!$ Orthogonormalite des vecteurs d'essai
-    call orthonorm(Vm(1,1,1),ndet,nvectot,tol_orth,tol_norm,.false.,prog_info)
+!!$ Orthogonormalite des vecteurs d'essai 
+    if (.not.prog_info%restart) call orthonorm(Vm(1,1,1),ndet,nvectot,tol_orth,tol_norm,.false.,prog_info)
     
   
 !!$ Davidson iterations ------------------------------------------------
@@ -255,10 +263,27 @@ contains
     etmp(:)     = 0.d0   ! energies temporaires   
     vconv(:)    = 0
 
-    nconv = 0
+    nconv = 0 !nb of converged vectors
     ncol  = 1
-    nvec  = nvectot
-
+    nvec  = nvectot !nb of vectors still to be converged
+    ! si restart on met les initialisations specifiques
+    if (prog_info%restart) then 
+       nconv = v_info%nconv !nb of converged vectors from the restart file
+       ener(:) = ener_info%ener(:)
+       if (nconv.ge.nvectot) then
+          if (prog_info%id_cpu.eq.0) then
+             write(f_output,'(" Number of converged vectors ",i4,&
+                  " equal to number of required vectors",i4," No diagonalisation needed." )') nconv, nvectot
+             call SASS_quit('>>> No need to diagonalise', f_output)
+          end if
+       end if
+       if (nconv.ne.0) then
+          nvec = nvectot - nconv
+          Vm(:,:,:) = 0.d0
+          Vm(1:ndet,1:nvec,1) = psi_SASS(1:ndet, 1+nconv:nvec+nconv)
+       end if
+    end if
+    
     !Fock and Integrals
     call intkind_H_all_init(intkindlist, int_info)
        
@@ -268,6 +293,7 @@ contains
        write(f_output,'(X,A6,X,A12,X,A12,X,A12,2X,A)') &
             '#Iter', 'tallCPU (s)', 'tCPU0 (s)', 'tWall0 (s)', "Energies"
        write(*,*) '#tCPU (s) : tWall (s) : tCPU/tWall : intkind : DblockI : DblockJ : spin'
+       write(*,*) 
 #else
        write(f_output,'(X,A6,X,A12,X,A12,2X,A)') '#Iter', 'tCPU (s)', 'tWall (s)', "Energies"
        write(667,*) '#tCPU (s) : tWall (s) : tCPU/tWall : intkind : DblockI : DblockJ : spin'
@@ -279,11 +305,15 @@ contains
     !----------------------------------------------------------------------------------------------
     daviter: do iter = 1, NitDavid !---------------------------------------------------------------
 
+       !call recouvrement(Vm(1,1,1), psi_SASS(1,1), ndet,nvec,nconv, f_output, "Vm avant David"," psi_SASS")
        call david_iter_mem(H_dav, ener, ener_np1, deltaE, vect, psi_SASS, &
             Vm, Wm, Dimheff, ncol, fock, hdiag, rspin, &
             det, o_info, int_info, v_info, prog_info, nelact, ndet, nvec, nconv, &
             hcase_info, pcase_info, bdav_info, ener_info, Hmat, intkindlist, mpilist, iter)
-
+       
+       !call recouvrement(Vm(1,1,1), psi_SASS(1,1), ndet,nvec,nconv, f_output, "Vm apres David"," psi_SASS")
+       !call recouvrement(Wm(1,1,1), psi_SASS(1,1), ndet,nvec,nconv, f_output, "Wm apres David"," psi_SASS")
+       
        if (prog_info%iprintHmat .gt. 0) goto 902      
           
        ! On teste la convergence sur les energies
@@ -311,7 +341,7 @@ contains
        ! On met à jour Ener
        Ener(1+nconv:nvectot) = Ener_np1(1:nvec)
 
-       ! Sauvegarde des vecteurs pour restart 
+       ! Sauvegarde des vecteurs pour restart
        rewind(f_restart)
        write(f_restart) nvectot,ndet,nconv, iter+iter0
        write(f_restart) Ecoeur,PotNuc
@@ -320,6 +350,14 @@ contains
           write(f_restart) (psi_SASS(idet,ivec),idet=1,ndet)
        end do
        flush(f_restart)
+            
+       !---  debug
+       if (debugdav) then
+          write(fdav,*) "psi_SASS - Iter =", iter
+          call wrt_vect(f_output,20,ndet,nvectot,psi_SASS)
+          write(fdav,*)
+       end if
+       !--- fin debug
 
        ! Si a on atteint Sizeheff alors on contracte tout sur les premiers vecteurs
        !  (une colonne)
@@ -358,12 +396,12 @@ contains
        end do
        write(f_Wm) Wm(:,:,1)
        flush(f_Wm)
-       write(f_output,*) "debug: je sauve les vecteurs de la colonne:",ncol
+       !write(f_output,*) "debug: je sauve les vecteurs de la colonne:",ncol
        flush(f_output)
        
        ! Puis  Calcul du nouveau vecteur de correction : (HV-EV)/(Ej-Hii) -> Vm(ncol+1)
        ncol = ncol + 1
-       if (debugdav) then 
+       if (debugdav .and. prog_info%id_cpu .eq. 0) then
           write(fdav,*) " vecteurs HV de l'iteration :", iter+iter0
           call wrtmat(Wm, ndet, min(10,ndet), nvec*ncol, fdav)
           write(fdav,*)
@@ -579,7 +617,7 @@ contains
     Ener_np1(:)  = 0.d0
     !if (.not.debugmat) then
 #ifdef VAR_MPI
-    if ((prog_info%id_cpu.eq.0).and.(prog_info%iprint.gt.0)) then
+    if ((prog_info%id_cpu.eq.0).and.(prog_info%iprint.gt.1)) then
        write(*,*) ''
        write(*,*) '>>>> Compute error'
        write(*,*) ''
@@ -753,7 +791,7 @@ contains
     close(f_mat)
     close(f_mat2)
     close(f_bmat)
-    close(666)
+    if (debugdav) close(fdav)
     close(667)
     close(f_Vm,status="delete")
     close(f_Wm,status="delete")
@@ -873,7 +911,7 @@ contains
 
 
 #ifdef VAR_MPI
-    if ((prog_info%id_cpu.eq.0).and.(prog_info%iprint.gt.0)) then
+    if ((prog_info%id_cpu.eq.0).and.(prog_info%iprint.gt.1)) then
        write(*,*) ''
        write(*,*) '>>>> Iteration', iter
        write(*,*) ''
@@ -930,6 +968,10 @@ contains
     endif
 #endif
 
+    ! Correction du defaut d'orthogonalisation de W avec vecteurs converges pour cause d'erreurs numeriques 
+    call orthog_Schmidt(Wm(1,1,1),psi_SASS(1,1),ndet,nvec,nconv,bdav_info%tol_orth)
+
+    
     ! Calcul de l'Hamiltonien effectif  H_dav
     ! bloc iteration iter -> <Vm(ivec,ncol)|H|Vm(jvec,ncol)> 
     nsh = (ncol-1)*nvec
@@ -967,7 +1009,7 @@ contains
     endif
 
     ! si iter=1 ecrire l'énergie des vecteurs d'essais
-    if (iter.eq.1) then
+    if (iter+iter0.eq.1) then
        do ivec = 1, nvec
           ener(ivec) = H_dav(ivec,ivec)
        end do
@@ -1018,6 +1060,7 @@ contains
        if (prog_info%id_cpu .eq. 0) &
             write(fdav,*) ">>> vecteurs vect de l'iteration :", iter + iter0
        call wrtmat(Vect,Dimheff,nvec*ncol,nvec*ncol, fdav)
+       write(fdav,*) 
     end if
 
 
@@ -1068,6 +1111,9 @@ contains
        end do
     end do
     deallocate(Vecttmp)
+!    call orthog_Schmidt(psi_SASS(1,1+nconv),psi_SASS(1,1),ndet,nvec,nconv,bdav_info%tol_orth)
+!    call orthonorm(psi_SASS(1,1+nconv),ndet,nvec,bdav_info%tol_orth,bdav_info%tol_norm,.true.,prog_info, f_output)
+
     
     if (debugdav) then 
        write(fdav,*) " vecteurs psi_SASS  de l'iteration :", iter+iter0, &
@@ -1382,7 +1428,7 @@ contains
 1000 continue
     
     if ((prog_info%id_cpu .eq. 0).and. (prog_info%iprint .ge. 0))  then
-       write(*,*) 'Final Sum Wall time on the dfifferent CPUs'
+       write(*,*) 'Final Sum Wall time on the different CPUs'
        write(*,*) sum_wall(:)
     endif