Normalize in input of Davidson

plugins/Full_CI/full_ci.irp.f

@@ -40,7 +40,7 @@ program full_ci
 
   integer :: n_det_before
   print*,'Beginning the selection ...'
-  E_CI_before = CI_energy
+  E_CI_before(1:N_states) = CI_energy(1:N_states)
   do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
     n_det_before = N_det
     call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
@@ -49,13 +49,16 @@ program full_ci
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
 
-    if (N_det > N_det_max) then
-       psi_det = psi_det_sorted
-       psi_coef = psi_coef_sorted
-       N_det = N_det_max
-       soft_touch N_det psi_det psi_coef
-    endif
     call diagonalize_CI
+
+    if (N_det > N_det_max) then
+      N_det = N_det_max
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
+    endif
+
+
     call save_wavefunction
     if(n_det_before == N_det)then
      selection_criterion = selection_criterion * 0.5d0
@@ -69,7 +72,6 @@ program full_ci
     print *,  'E(before)+PT2  = ', E_CI_before(k)+pt2(k)
     enddo
     print *,  '-----'
-    E_CI_before = CI_energy
     if(N_states.gt.1)then
      print*,'Variational Energy difference'
      do i = 2, N_states
@@ -82,8 +84,8 @@ program full_ci
       print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
      enddo
     endif
-    E_CI_before = CI_energy
+    E_CI_before(1:N_states) = CI_energy(1:N_states)
-    call ezfio_set_full_ci_energy(CI_energy)
+    call ezfio_set_full_ci_energy(CI_energy(1))
   enddo
    N_det = min(N_det_max,N_det)
    touch N_det psi_det psi_coef
@@ -101,7 +103,7 @@ program full_ci
     print *,  'E        = ', CI_energy
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
-    call ezfio_set_full_ci_energy_pt2(CI_energy+pt2)
+    call ezfio_set_full_ci_energy_pt2(CI_energy(1)+pt2(1))
    endif
    call save_wavefunction
   deallocate(pt2,norm_pert)

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -44,6 +44,8 @@ program fci_zmq
     PROVIDE  psi_coef
     PROVIDE  psi_det
     PROVIDE  psi_det_sorted
+
+    call diagonalize_CI
     
     if (N_det > N_det_max) then
       psi_det = psi_det_sorted
@@ -51,7 +53,6 @@ program fci_zmq
       N_det = N_det_max
       soft_touch N_det psi_det psi_coef
     endif
-    call diagonalize_CI
     call save_wavefunction
     
     print *,  'N_det          = ', N_det

plugins/Full_CI_ZMQ/selection.irp.f

@@ -93,10 +93,19 @@ subroutine select_connected(i_generator,E0,pt2,b)
       particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) )
       particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) )
 
-      hole_mask(k,1) = ior(generators_bitmask(k,1,s_hole,l), generators_bitmask(k,1,s_part,l))
+!      hole_mask(k,1) = ior(generators_bitmask(k,1,s_hole,l), generators_bitmask(k,1,s_part,l))
-      hole_mask(k,2) = ior(generators_bitmask(k,2,s_hole,l), generators_bitmask(k,2,s_part,l))
+!      hole_mask(k,2) = ior(generators_bitmask(k,2,s_hole,l), generators_bitmask(k,2,s_part,l))
-      particle_mask(k,:) = hole_mask(k,:)
+!      particle_mask(k,1) = hole_mask(k,1)
+!      particle_mask(k,2) = hole_mask(k,2)
     enddo
+    print *, 'det'
+    call debug_det(psi_det_generators(1,1,i_generator),N_int)
+    print *, 'hole'
+    call debug_det(hole_mask,N_int)
+    print *, 'particle_mask'
+    call debug_det(particle_mask,N_int)
+    print *,  ''
+    pause
     call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b)
     call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b)
   enddo

plugins/Perturbation/perturbation.template.f

@@ -131,7 +131,6 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
 ! TODO OLD
 !        if(is_connected_to(buffer(1,1,i), microlist_gen(:,:,1:ptr_microlist_gen(1)-1), Nint, N_microlist_gen(0))) then
 ! TODO OLD
-        ASSERT ( N_microlist_gen(0) <= buffer_size)
         if(is_connected_to(buffer(1,1,i), microlist_gen(1,1,1), Nint, N_microlist_gen(0))) then
           cycle
         end if

src/Davidson/diagonalization.irp.f

@@ -381,25 +381,28 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
   ! ===================
   
   converged = .False.
+  
+  do k=1,N_st_diag
 
-  do k=N_st+1,N_st_diag
+    if (k > N_st) then
-    double precision :: r1, r2
+      do i=1,sze
-    do i=1,sze
+        double precision               :: r1, r2
-      call random_number(r1)
+        call random_number(r1)
-      call random_number(r2)
+        call random_number(r2)
-      u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2)
+        u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2)
-    enddo
+      enddo
+    endif
     
     ! Gram-Schmidt
     ! ------------
-    call dgemv('T',sze,k-1,1.d0,u_in,size(u_in,1),   &
+    call dgemv('T',sze,k-1,1.d0,u_in,size(u_in,1),                   &
         u_in(1,k),1,0.d0,c,1)
-    call dgemv('N',sze,k-1,-1.d0,u_in,size(u_in,1),        &
+    call dgemv('N',sze,k-1,-1.d0,u_in,size(u_in,1),                  &
         c,1,1.d0,u_in(1,k),1)
-    call normalize( u_in(1,k), sze )
+    call normalize(u_in(1,k),sze)
-
   enddo
 
+
   
   do while (.not.converged)
     
@@ -461,8 +464,8 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
 !          enddo
 !        enddo
 !      enddo
-
+!
-
+!
       call dgemm('N','N', sze, N_st_diag, N_st_diag*iter,            &
           1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1))
       call dgemm('N','N',sze,N_st_diag,N_st_diag*iter,               &
@@ -511,19 +514,19 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
 !            enddo
 !          enddo
 !        enddo
-
+!
         call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1),  &
               U(1,k,iter+1),1,0.d0,c,1)
         call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
               c,1,1.d0,U(1,k,iter+1),1)
-
+!
 !        do l=1,k-1
 !          c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
 !          do i=1,sze
 !            U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter+1)
 !          enddo
 !        enddo
-
+!
         call dgemv('T',sze,k-1,1.d0,U(1,1,iter+1),size(U,1),   &
             U(1,k,iter+1),1,0.d0,c,1)
         call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1),        &

src/Determinants/determinants.irp.f

@@ -90,6 +90,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
   logical                        :: exists
   character*64                   :: label
   
+  psi_det = 0_bit_kind
   if (read_wf) then
     call ezfio_has_determinants_N_int(exists)
     if (exists) then
@@ -255,7 +256,7 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
   character*(64)                 :: label
 
   psi_coef = 0.d0
-  do i=1,N_states
+  do i=1,min(N_states,psi_det_size)
     psi_coef(i,i) = 1.d0
   enddo
 
@@ -331,7 +332,6 @@ END_PROVIDER
    iorder(i) = i
  enddo
  call dsort(psi_average_norm_contrib_sorted,iorder,N_det)
- !DIR$ IVDEP
  do i=1,N_det
   do j=1,N_int
     psi_det_sorted(j,1,i) = psi_det(j,1,iorder(i))

src/Determinants/s2.irp.f

@@ -381,7 +381,8 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nsta
 
   call lapack_diagd(s2_eigvalues,eigvectors,s2,nstates,nstates)
   print*,'Eigenvalues'
-  double precision :: t(nstates), iorder(nstates)
+  double precision :: t(nstates)
+  integer :: iorder(nstates)
   do i = 1, nstates
     t(i) = dabs(s2_eigvalues(i))
     iorder(i) = i

tests/bats/qp.bats

@@ -68,7 +68,7 @@ function run_FCI() {
   ezfio set_file $1
   ezfio set perturbation do_pt2_end True
   ezfio set determinants n_det_max $2
-  ezfio set determinants threshold_davidson 1.e-10
+  ezfio set davidson threshold_davidson 1.e-10
 
   qp_run full_ci $1 
   energy="$(ezfio get full_ci energy)"
@@ -83,7 +83,7 @@ function run_all_1h_1p() {
   ezfio set_file $1
   ezfio set determinants n_det_max $2
   ezfio set perturbation pt2_max $3
-  ezfio set determinants threshold_davidson 1.e-10
+  ezfio set davidson threshold_davidson 1.e-10
 
   qp_run all_1h_1p $1 | tee $1.F1h1p.out
   energy="$(ezfio get all_singles energy)"