non_sym dress: comb

src/davidson/EZFIO.cfg

@@ -1,71 +1,18 @@
-[threshold_davidson]
-type: Threshold
-doc: Thresholds of Davidson's algorithm if threshold_davidson_from_pt2 is false.
-interface: ezfio,provider,ocaml
-default: 1.e-10
-
-[threshold_nonsym_davidson]
-type: Threshold
-doc: Thresholds of non-symetric Davidson's algorithm
-interface: ezfio,provider,ocaml
-default: 1.e-10
-
-[threshold_davidson_from_pt2]
-type: logical
-doc: Thresholds of Davidson's algorithm is set to E(rPT2)*threshold_davidson_from_pt2
-interface: ezfio,provider,ocaml
-default: false
-
-[n_states_diag]
-type: States_number
-doc: Controls the number of states to consider during the Davdison diagonalization. The number of states is n_states * n_states_diag
-default: 4
-interface: ezfio,ocaml
-
-[davidson_sze_max]
-type: Strictly_positive_int
-doc: Number of micro-iterations before re-contracting
-default: 15
-interface: ezfio,provider,ocaml
-
-[state_following]
-type: logical
-doc: If |true|, the states are re-ordered to match the input states
-default: False
-interface: ezfio,provider,ocaml
-
-[disk_based_davidson]
-type: logical
-doc: If |true|, a memory-mapped file may be used to store the W and S2 vectors if not enough RAM is available
-default: True
-interface: ezfio,provider,ocaml
-
 [csf_based]
 type: logical
 doc: If |true|, use the CSF-based algorithm
 default: False
 interface: ezfio,provider,ocaml
 
-[distributed_davidson]
-type: logical
-doc: If |true|, use the distributed algorithm
-default: True 
-interface: ezfio,provider,ocaml
-
 [only_expected_s2]
 type: logical
 doc: If |true|, use filter out all vectors with bad |S^2| values
 default: True
 interface: ezfio,provider,ocaml
 
-[n_det_max_full]
-type: Det_number_max
-doc: Maximum number of determinants where |H| is fully diagonalized 
-interface: ezfio,provider,ocaml
-default: 1000
-
 [without_diagonal]
 type: logical
 doc: If |true|, don't use denominator
 default: False
 interface: ezfio,provider,ocaml
+

src/davidson/NEED

@@ -1 +1,2 @@
 csf
+davidson_keywords

src/davidson/davidson_parallel.irp.f

@@ -546,21 +546,6 @@ end
 
 
 
-BEGIN_PROVIDER [ integer, nthreads_davidson ]
- implicit none
- BEGIN_DOC
- ! Number of threads for Davidson
- END_DOC
- nthreads_davidson = nproc
- character*(32) :: env
- call getenv('QP_NTHREADS_DAVIDSON',env)
- if (trim(env) /= '') then
-   read(env,*) nthreads_davidson
-   call write_int(6,nthreads_davidson,'Target number of threads for <Psi|H|Psi>')
- endif
-END_PROVIDER
-
-
 integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -14,15 +14,6 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
   endif
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, threshold_davidson_pt2 ]
- implicit none
- BEGIN_DOC
- ! Threshold of Davidson's algorithm, using PT2 as a guide
- END_DOC
- threshold_davidson_pt2 = threshold_davidson
-
-END_PROVIDER
-
 
 
 BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ]
@@ -66,7 +57,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   double precision, allocatable  :: H_jj(:)
 
   double precision, external     :: diag_H_mat_elem, diag_S_mat_elem
-  integer                        :: i,k
+  integer                        :: i,k,l
   ASSERT (N_st > 0)
   ASSERT (sze > 0)
   ASSERT (Nint > 0)
@@ -87,9 +78,14 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
 
   if (dressing_state > 0) then
     do k=1,N_st
+
       do i=1,sze
-        H_jj(i)  += u_in(i,k) * dressing_column_h(i,k)
+        H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
       enddo
+
+      !l = dressed_column_idx(k)
+      !H_jj(l) += u_in(l,k) * dressing_column_h(l,k)
+
     enddo
   endif
 

src/davidson/diagonalization_nonsym_h_dressed.irp.f

@@ -0,0 +1,541 @@
+
+! ---
+
+subroutine davidson_diag_nonsym_h(dets_in, u_in, dim_in, energies, sze, N_st, N_st_diag, Nint, dressing_state, converged)
+
+  BEGIN_DOC
+  !
+  ! non-sym Davidson diagonalization.
+  !
+  ! dets_in : bitmasks corresponding to determinants
+  !
+  ! u_in : guess coefficients on the various states. Overwritten on exit
+  !
+  ! dim_in : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  ! Initial guess vectors are not necessarily orthonormal
+  !
+  END_DOC
+
+  use bitmasks
+
+  implicit none
+
+  integer,           intent(in)    :: dim_in, sze, N_st, N_st_diag, Nint
+  integer,           intent(in)    :: dressing_state
+  integer(bit_kind), intent(in)    :: dets_in(Nint,2,sze)
+  logical,           intent(out)   :: converged
+  double precision,  intent(out)   :: energies(N_st_diag)
+  double precision,  intent(inout) :: u_in(dim_in,N_st_diag)
+
+  integer                          :: i, k, l
+  double precision                 :: f
+  double precision, allocatable    :: H_jj(:)
+
+  double precision, external       :: diag_H_mat_elem
+
+  ASSERT (N_st > 0)
+  ASSERT (sze  > 0)
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  PROVIDE mo_two_e_integrals_in_map
+
+  allocate(H_jj(sze))
+
+  H_jj(1) = diag_H_mat_elem(dets_in(1,1,1), Nint)
+ !$OMP PARALLEL DEFAULT(NONE)           &
+ !$OMP SHARED(sze, H_jj, dets_in, Nint) &
+ !$OMP PRIVATE(i)
+ !$OMP DO SCHEDULE(static)
+  do i = 2, sze
+    H_jj(i) = diag_H_mat_elem(dets_in(1,1,i), Nint)
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  if(dressing_state > 0) then
+    do k = 1, N_st
+      do l = 1, N_st
+        f = overlap_states_inv(k,l)
+
+        !do i = 1, N_det
+        !  H_jj(i) += f * dressing_delta(i,k) * psi_coef(i,l)
+        do i = 1, dim_in
+          H_jj(i) += f * dressing_delta(i,k) * u_in(i,l)
+        enddo
+
+      enddo
+    enddo
+  endif
+
+  call davidson_diag_nonsym_hjj(dets_in, u_in, H_jj, energies, dim_in, sze, N_st, N_st_diag, Nint, dressing_state, converged)
+
+  deallocate(H_jj)
+
+end subroutine davidson_diag_nonsym_h
+
+! ---
+
+subroutine davidson_diag_nonsym_hjj(dets_in, u_in, H_jj, energies, dim_in, sze, N_st, N_st_diag_in, Nint, dressing_state, converged)
+
+  BEGIN_DOC
+  !
+  ! non-sym Davidson diagonalization with specific diagonal elements of the H matrix
+  !
+  ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
+  !
+  ! dets_in : bitmasks corresponding to determinants
+  !
+  ! u_in : guess coefficients on the various states. Overwritten on exit
+  !
+  ! dim_in : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze
+  !
+  ! Initial guess vectors are not necessarily orthonormal
+  !
+  END_DOC
+
+  include 'constants.include.F'
+
+  use bitmasks
+  use mmap_module
+
+  implicit none
+
+  integer,           intent(in)    :: dim_in, sze, N_st, N_st_diag_in, Nint
+  integer,           intent(in)    :: dressing_state
+  integer(bit_kind), intent(in)    :: dets_in(Nint,2,sze)
+  double precision,  intent(in)    :: H_jj(sze)
+  double precision,  intent(out)   :: energies(N_st_diag_in)
+  logical,           intent(inout) :: converged
+  double precision,  intent(inout) :: u_in(dim_in,N_st_diag_in)
+
+  logical                          :: disk_based
+  character*(16384)                :: write_buffer
+  integer                          :: i, j, k, l, m
+  integer                          :: iter, N_st_diag, itertot, shift, shift2, itermax, istate
+  integer                          :: nproc_target
+  integer                          :: order(N_st_diag_in)
+  integer                          :: maxab
+  double precision                 :: rss
+  double precision                 :: cmax
+  double precision                 :: to_print(2,N_st)
+  double precision                 :: r1, r2
+  double precision                 :: f
+  double precision, allocatable    :: y(:,:), h(:,:), lambda(:)
+  double precision, allocatable    :: s_tmp(:,:), u_tmp(:,:)
+  double precision, allocatable    :: residual_norm(:)
+  double precision, allocatable    :: U(:,:), overlap(:,:)
+  double precision, pointer        :: W(:,:)
+
+  double precision, external       :: u_dot_u
+
+
+  N_st_diag = N_st_diag_in
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
+  if(N_st_diag*3 > sze) then
+    print *,  'error in Davidson :'
+    print *,  'Increase n_det_max_full to ', N_st_diag*3
+    stop -1
+  endif
+
+  itermax = max(2, min(davidson_sze_max, sze/N_st_diag)) + 1
+  itertot = 0
+
+  if(state_following) then
+    allocate(overlap(N_st_diag*itermax, N_st_diag*itermax))
+  else
+    allocate(overlap(1,1))  ! avoid 'if' for deallocate
+  endif
+  overlap = 0.d0
+
+  PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse threshold_davidson_pt2 threshold_davidson_from_pt2
+  PROVIDE threshold_nonsym_davidson
+
+  call write_time(6)
+  write(6,'(A)') ''
+  write(6,'(A)') 'Davidson Diagonalization'
+  write(6,'(A)') '------------------------'
+  write(6,'(A)') ''
+
+  ! Find max number of cores to fit in memory
+  ! -----------------------------------------
+
+  nproc_target = nproc
+  maxab = max(N_det_alpha_unique, N_det_beta_unique) + 1
+
+  m=1
+  disk_based = .False.
+  call resident_memory(rss)
+  do
+    r1 = 8.d0 *                                   &! bytes
+         ( dble(sze)*(N_st_diag*itermax)          &! U
+         + 1.0d0*dble(sze*m)*(N_st_diag*itermax)  &! W
+         + 3.0d0*(N_st_diag*itermax)**2           &! h,y,s_tmp
+         + 1.d0*(N_st_diag*itermax)               &! lambda
+         + 1.d0*(N_st_diag)                       &! residual_norm
+                                                   ! In H_u_0_nstates_zmq
+         + 2.d0*(N_st_diag*N_det)                 &! u_t, v_t, on collector
+         + 2.d0*(N_st_diag*N_det)                 &! u_t, v_t, on slave
+         + 0.5d0*maxab                            &! idx0 in H_u_0_nstates_openmp_work_*
+         + nproc_target *                         &! In OMP section
+           ( 1.d0*(N_int*maxab)                   &! buffer
+           + 3.5d0*(maxab) )                      &! singles_a, singles_b, doubles, idx
+         ) / 1024.d0**3
+
+    if(nproc_target == 0) then
+      call check_mem(r1, irp_here)
+      nproc_target = 1
+      exit
+    endif
+
+    if(r1+rss < qp_max_mem) then
+      exit
+    endif
+
+    if(itermax > 4) then
+      itermax = itermax - 1
+    else if(m==1 .and. disk_based_davidson) then
+      m = 0
+      disk_based = .True.
+      itermax = 6
+    else
+      nproc_target = nproc_target - 1
+    endif
+
+  enddo
+
+  nthreads_davidson = nproc_target
+  TOUCH nthreads_davidson
+
+  call write_int(6, N_st, 'Number of states')
+  call write_int(6, N_st_diag, 'Number of states in diagonalization')
+  call write_int(6, sze, 'Number of determinants')
+  call write_int(6, nproc_target, 'Number of threads for diagonalization')
+  call write_double(6, r1, 'Memory(Gb)')
+  if(disk_based) then
+    print *, 'Using swap space to reduce RAM'
+  endif
+
+  !---------------
+
+  write(6,'(A)') ''
+  write_buffer = '====='
+  do i = 1, N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6, '(A)') write_buffer(1:6+41*N_st)
+  write_buffer = 'Iter'
+  do i = 1, N_st
+    write_buffer = trim(write_buffer)//'       Energy        Residual '
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = '====='
+  do i = 1, N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+
+
+  if(disk_based) then
+    ! Create memory-mapped files for W and S
+    type(c_ptr) :: ptr_w, ptr_s
+    integer :: fd_s, fd_w
+    call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
+        8, fd_w, .False., ptr_w)
+    call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
+  else
+    allocate(W(sze,N_st_diag*itermax))
+  endif
+
+  allocate(                                                          &
+      ! Large
+      U(sze,N_st_diag*itermax),                                      &
+      ! Small
+      h(N_st_diag*itermax,N_st_diag*itermax),                        &
+      y(N_st_diag*itermax,N_st_diag*itermax),                        &
+      s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
+      residual_norm(N_st_diag),                                      &
+      lambda(N_st_diag*itermax),                                     &
+      u_tmp(N_st,N_st_diag))
+
+  h = 0.d0
+  U = 0.d0
+  y = 0.d0
+  s_tmp = 0.d0
+
+
+  ASSERT (N_st > 0)
+  ASSERT (N_st_diag >= N_st)
+  ASSERT (sze > 0)
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+
+  ! Davidson iterations
+  ! ===================
+
+  converged = .False.
+
+  do k = N_st+1, N_st_diag
+    do i = 1, sze
+        call random_number(r1)
+        call random_number(r2)
+        r1 = dsqrt(-2.d0*dlog(r1))
+        r2 = dtwo_pi*r2
+        u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
+    enddo
+    u_in(k,k) = u_in(k,k) + 10.d0
+  enddo
+  do k = 1, N_st_diag
+    call normalize(u_in(1,k), sze)
+  enddo
+
+  do k = 1, N_st_diag
+    do i = 1, sze
+      U(i,k) = u_in(i,k)
+    enddo
+  enddo
+
+
+  do while (.not.converged)
+    itertot = itertot + 1
+    if(itertot == 8) then
+      exit
+    endif
+
+    do iter = 1, itermax-1
+
+      shift  = N_st_diag*(iter-1)
+      shift2 = N_st_diag*iter
+
+!      if( (iter > 1) .or. (itertot == 1) ) then
+
+        ! Gram-Schmidt to orthogonalize all new guess with the previous vectors 
+        call ortho_qr(U, size(U, 1), sze, shift2)
+        call ortho_qr(U, size(U, 1), sze, shift2)
+
+        ! Compute |W_k> = \sum_i |i><i|H|u_k>
+        ! -----------------------------------
+
+        if( (sze > 100000) .and. distributed_davidson ) then
+          call H_u_0_nstates_zmq   (W(1,shift+1), U(1,shift+1), N_st_diag, sze)
+        else
+          call H_u_0_nstates_openmp(W(1,shift+1), U(1,shift+1), N_st_diag, sze)
+        endif
+!      else
+!         ! Already computed in update below
+!         continue
+!      endif
+
+      if(dressing_state > 0) then
+
+        call dgemm( 'T', 'N', N_st, N_st_diag, sze, 1.d0                   &
+                  , psi_coef, size(psi_coef, 1), U(1, shift+1), size(U, 1) &
+                  , 0.d0, u_tmp, size(u_tmp, 1))
+
+        do istate = 1, N_st_diag
+          do k = 1, N_st
+            do l = 1, N_st
+             f = overlap_states_inv(k,l)
+              do i = 1, sze
+                W(i,shift+istate) += f * dressing_delta(i,k) * u_tmp(l,istate)
+              enddo
+            enddo
+          enddo
+        enddo
+
+      endif
+
+      ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
+      ! -------------------------------------------
+
+      call dgemm( 'T', 'N', shift2, shift2, sze, 1.d0 &
+                , U, size(U, 1), W, size(W, 1)        &
+                , 0.d0, h, size(h, 1))
+
+      ! Diagonalize h
+      ! ---------------
+      call diag_nonsym_right(shift2, h(1,1), size(h, 1), y(1,1), size(y, 1), lambda(1), size(lambda, 1))
+
+
+      if (state_following) then
+
+        overlap = -1.d0
+        do k = 1, shift2
+          do i = 1, shift2
+            overlap(k,i) = dabs(y(k,i))
+          enddo
+        enddo
+        do k = 1, N_st
+          cmax = -1.d0
+          do i = 1, N_st
+            if(overlap(i,k) > cmax) then
+              cmax = overlap(i,k)
+              order(k) = i
+            endif
+          enddo
+          do i = 1, N_st_diag
+            overlap(order(k),i) = -1.d0
+          enddo
+        enddo
+        overlap = y
+        do k = 1, N_st
+          l = order(k)
+          if (k /= l) then
+            y(1:shift2,k) = overlap(1:shift2,l)
+          endif
+        enddo
+        do k = 1, N_st
+          overlap(k,1) = lambda(k)
+        enddo
+
+      endif
+
+      ! Express eigenvectors of h in the determinant basis
+      ! --------------------------------------------------
+
+      call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
+                , U, size(U, 1), y, size(y, 1)           &
+                , 0.d0, U(1,shift2+1), size(U, 1))
+
+      do k = 1, N_st_diag
+        call normalize(U(1,shift2+k), sze)
+      enddo
+
+      call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
+                , W, size(W, 1), y, size(y, 1)           &
+                , 0.d0, W(1,shift2+1), size(W,1))
+
+      ! Compute residual vector and davidson step
+      ! -----------------------------------------
+
+      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
+      do k = 1, N_st_diag
+        do i = 1, sze
+          U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k)) / max(H_jj(i)-lambda(k), 1.d-2)
+        enddo
+
+        if(k <= N_st) then
+          residual_norm(k) = u_dot_u(U(1,shift2+k), sze)
+          to_print(1,k) = lambda(k) + nuclear_repulsion
+          to_print(2,k) = residual_norm(k)
+        endif
+      enddo
+      !$OMP END PARALLEL DO
+
+      if((itertot>1).and.(iter == 1)) then
+        !don't print
+        continue
+      else
+        write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, E11.3))') iter-1, to_print(1:2,1:N_st)
+      endif
+
+      ! Check convergence
+      if(iter > 1) then
+        if(threshold_davidson_from_pt2) then
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson_pt2
+        else
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_nonsym_davidson
+        endif
+      endif
+
+      do k = 1, N_st
+        if(residual_norm(k) > 1.d8) then
+          print *, 'Davidson failed'
+          stop -1
+        endif
+      enddo
+      if(converged) then
+        exit
+      endif
+
+      logical, external :: qp_stop
+      if(qp_stop()) then
+        converged = .True.
+        exit
+      endif
+
+
+    enddo
+
+    ! Re-contract U and update W
+    ! --------------------------------
+
+    call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
+              , W, size(W, 1), y, size(y, 1)           & 
+              , 0.d0, u_in, size(u_in, 1))
+    do k = 1, N_st_diag
+      do i = 1, sze
+        W(i,k) = u_in(i,k)
+      enddo
+    enddo
+
+    call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
+              , U, size(U, 1), y, size(y, 1), 0.d0     &
+              , u_in, size(u_in, 1))
+
+    do k = 1, N_st_diag
+      do i = 1, sze
+        U(i,k) = u_in(i,k)
+      enddo
+    enddo
+
+  enddo
+
+
+  call nullify_small_elements(sze, N_st_diag, U, size(U, 1), threshold_davidson_pt2)
+  do k = 1, N_st_diag
+    do i = 1, sze
+      u_in(i,k) = U(i,k)
+    enddo
+  enddo
+
+  do k = 1, N_st_diag
+    energies(k) = lambda(k)
+  enddo
+  write_buffer = '======'
+  do i = 1, N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') trim(write_buffer)
+  write(6,'(A)') ''
+  call write_time(6)
+
+  if(disk_based) then
+    ! Remove temp files
+    integer, external :: getUnitAndOpen
+    call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 8, fd_w, ptr_w )
+    fd_w = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_w','r')
+    close(fd_w,status='delete')
+  else
+    deallocate(W)
+  endif
+
+  deallocate (                                                       &
+      residual_norm,                                                 &
+      U, overlap,                                                    &
+      h, y, s_tmp,                                                   &
+      lambda,                                                        &
+      u_tmp                                                          &
+      )
+  FREE nthreads_davidson
+
+end subroutine davidson_diag_nonsym_hjj
+
+! ---
+
+
+
+
+
+
+

src/davidson/overlap_states.irp.f

@@ -0,0 +1,40 @@
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, overlap_states,     (N_states,N_states) ]
+&BEGIN_PROVIDER [ double precision, overlap_states_inv, (N_states,N_states) ]
+
+  BEGIN_DOC
+  !
+  ! S_kl = ck.T x cl 
+  !      = psi_coef(:,k).T x psi_coef(:,l)
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i
+  double precision :: o_tmp
+
+  if(N_states == 1) then
+
+    o_tmp = 0.d0
+    do i = 1, N_det
+      o_tmp = o_tmp + psi_coef(i,1) * psi_coef(i,1)
+    enddo
+    overlap_states    (1,1) = o_tmp
+    overlap_states_inv(1,1) = 1.d0 / o_tmp
+
+  else
+
+    call dgemm( 'T', 'N', N_states, N_states, N_det, 1.d0                &
+              , psi_coef, size(psi_coef, 1), psi_coef, size(psi_coef, 1) &
+              , 0.d0, overlap_states, size(overlap_states, 1)            )
+
+    call get_inverse(overlap_states, N_states, N_states, overlap_states_inv, N_states)
+
+  endif
+
+END_PROVIDER
+
+! ---
+

src/davidson_dressed/nonsym_diagonalize_ci.irp.f

@@ -0,0 +1,188 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, CI_energy_nonsym_dressed, (N_states_diag) ]
+
+  BEGIN_DOC
+  ! N_states lowest eigenvalues of the CI matrix
+  END_DOC
+
+  implicit none
+  integer       :: j
+  character*(8) :: st
+
+  call write_time(6)
+  do j = 1, min(N_det, N_states_diag)
+    CI_energy_nonsym_dressed(j) = CI_electronic_energy_nonsym_dressed(j) + nuclear_repulsion
+  enddo
+
+  do j = 1, min(N_det, N_states)
+    write(st, '(I4)') j
+    call write_double(6, CI_energy_nonsym_dressed(j), 'Energy of state '//trim(st))
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, CI_electronic_energy_nonsym_dressed, (N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_eigenvectors_nonsym_dressed, (N_det,N_states_diag) ]
+
+  BEGIN_DOC
+  ! Eigenvectors/values of the CI matrix
+  END_DOC
+
+  implicit none
+  logical                       :: converged
+  integer                       :: i, j, k
+  integer                       :: i_other_state
+  integer                       :: i_state
+  logical,          allocatable :: good_state_array(:)
+  integer,          allocatable :: index_good_state_array(:)
+  double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
+
+  PROVIDE threshold_nonsym_davidson nthreads_davidson
+
+  ! Guess values for the "N_states" states of the CI_eigenvectors_nonsym_dressed
+  do j = 1, min(N_states, N_det)
+    do i = 1, N_det
+      CI_eigenvectors_nonsym_dressed(i,j) = psi_coef(i,j)
+    enddo
+  enddo
+
+  do j = min(N_states, N_det)+1, N_states_diag
+    do i = 1, N_det
+      CI_eigenvectors_nonsym_dressed(i,j) = 0.d0
+    enddo
+  enddo
+
+  ! ---
+
+  if(diag_algorithm == "Davidson") then
+
+    ASSERT(n_states_diag .lt. n_states)
+
+    do j = 1, min(N_states, N_det)
+      do i = 1, N_det
+        CI_eigenvectors_nonsym_dressed(i,j) = psi_coef(i,j)
+      enddo
+    enddo
+
+    converged = .False.
+    call davidson_diag_nonsym_h( psi_det, CI_eigenvectors_nonsym_dressed &
+                               , size(CI_eigenvectors_nonsym_dressed, 1) &
+                               , CI_electronic_energy_nonsym_dressed     &
+                               , N_det, min(N_det, N_states), min(N_det, N_states_diag), N_int, 1, converged )
+
+  else if(diag_algorithm == "Lapack") then
+
+    allocate(eigenvectors(size(H_matrix_nonsym_dressed, 1),N_det))
+    allocate(eigenvalues(N_det))
+
+    call diag_nonsym_right( N_det, H_matrix_nonsym_dressed, size(H_matrix_nonsym_dressed, 1)       &
+                          , eigenvectors, size(eigenvectors, 1), eigenvalues, size(eigenvalues, 1) )
+
+    CI_electronic_energy_nonsym_dressed(:) = 0.d0
+
+    ! Select the "N_states_diag" states of lowest energy
+    do j = 1, min(N_det, N_states_diag)
+      do i = 1, N_det
+        CI_eigenvectors_nonsym_dressed(i,j) = eigenvectors(i,j)
+      enddo
+      CI_electronic_energy_nonsym_dressed(j) = eigenvalues(j)
+    enddo
+
+    deallocate(eigenvectors, eigenvalues)
+
+    ! --- ---
+
+  endif
+
+  ! ---
+
+END_PROVIDER
+
+! ---
+
+subroutine diagonalize_CI_nonsym_dressed()
+
+  BEGIN_DOC
+  !  Replace the coefficients of the CI states by the coefficients of the
+  !  eigenstates of the CI matrix
+  END_DOC
+
+  implicit none
+  integer :: i, j
+
+  PROVIDE dressing_delta
+
+  do j = 1, N_states
+    do i = 1, N_det
+      psi_coef(i,j) = CI_eigenvectors_nonsym_dressed(i,j)
+    enddo
+  enddo
+
+  SOFT_TOUCH psi_coef
+
+end subroutine diagonalize_CI_nonsym_dressed
+
+! ---
+
+BEGIN_PROVIDER [ double precision, H_matrix_nonsym_dressed, (N_det,N_det) ]
+
+  BEGIN_DOC
+  ! Dressed H with Delta_ij
+  END_DOC
+
+  implicit none
+  integer          :: i, j, l, k
+  double precision :: f
+
+  H_matrix_nonsym_dressed(1:N_det,1:N_det) = h_matrix_all_dets(1:N_det,1:N_det)
+
+  if(N_states == 1) then
+
+!    !symmetric formula
+!    l = dressed_column_idx(1)
+!    f = 1.0d0/psi_coef(l,1)
+!    do i=1,N_det
+!      h_matrix_nonsym_dressed(i,l) +=  dressing_column_h(i,1) *f
+!      h_matrix_nonsym_dressed(l,i) +=  dressing_column_h(i,1) *f
+!    enddo
+
+!    l = dressed_column_idx(1)
+!    f = 1.0d0 / psi_coef(l,1)
+!    do j = 1, N_det
+!      H_matrix_nonsym_dressed(j,l) += f * dressing_delta(j,1) 
+!    enddo
+
+    k = 1
+    l = 1
+    f = overlap_states_inv(k,l)
+    do j = 1, N_det
+      do i = 1, N_det
+        H_matrix_nonsym_dressed(i,j) = H_matrix_nonsym_dressed(i,j) + f * dressing_delta(i,k) * psi_coef(j,l)
+      enddo
+    enddo
+
+  else
+
+    do k = 1, N_states
+      do l = 1, N_states
+        f = overlap_states_inv(k,l)
+
+        do j = 1, N_det
+          do i = 1, N_det
+            H_matrix_nonsym_dressed(i,j) = H_matrix_nonsym_dressed(i,j) + f * dressing_delta(i,k) * psi_coef(j,l)
+          enddo
+        enddo
+
+      enddo
+    enddo
+
+  endif
+
+END_PROVIDER
+
+! ---
+

src/davidson_keywords/EZFIO.cfg

@@ -0,0 +1,54 @@
+[threshold_davidson]
+type: Threshold
+doc: Thresholds of Davidson's algorithm if threshold_davidson_from_pt2 is false.
+interface: ezfio,provider,ocaml
+default: 1.e-10
+
+[threshold_nonsym_davidson]
+type: Threshold
+doc: Thresholds of non-symetric Davidson's algorithm
+interface: ezfio,provider,ocaml
+default: 1.e-10
+
+[davidson_sze_max]
+type: Strictly_positive_int
+doc: Number of micro-iterations before re-contracting
+default: 15
+interface: ezfio,provider,ocaml
+
+[state_following]
+type: logical
+doc: If |true|, the states are re-ordered to match the input states
+default: False
+interface: ezfio,provider,ocaml
+
+[disk_based_davidson]
+type: logical
+doc: If |true|, a memory-mapped file may be used to store the W and S2 vectors if not enough RAM is availabl
+default: True
+interface: ezfio,provider,ocaml
+
+[n_states_diag]
+type: States_number
+doc: Controls the number of states to consider during the Davdison diagonalization. The number of states is n_states * n_states_diag
+default: 4
+interface: ezfio,ocaml
+
+[n_det_max_full]
+type: Det_number_max
+doc: Maximum number of determinants where |H| is fully diagonalized
+interface: ezfio,provider,ocaml
+default: 1000
+
+[threshold_davidson_from_pt2]
+type: logical
+doc: Thresholds of Davidson's algorithm is set to E(rPT2)*threshold_davidson_from_pt2
+interface: ezfio,provider,ocaml
+default: false
+
+[distributed_davidson]
+type: logical
+doc: If |true|, use the distributed algorithm
+default: True
+interface: ezfio,provider,ocaml
+

src/davidson_keywords/README.rst

@@ -0,0 +1,5 @@
+=================
+davidson_keywords
+=================
+
+Keywords used for Davidson algorithms.

src/davidson_keywords/input.irp.f

@@ -1,3 +1,6 @@
+
+! ---
+
 BEGIN_PROVIDER [ integer, n_states_diag  ]
   implicit none
   BEGIN_DOC
@@ -8,11 +11,11 @@ BEGIN_PROVIDER [ integer, n_states_diag  ]
   PROVIDE ezfio_filename
   if (mpi_master) then
 
-    call ezfio_has_davidson_n_states_diag(has)
+    call ezfio_has_davidson_keywords_n_states_diag(has)
     if (has) then
-      call ezfio_get_davidson_n_states_diag(n_states_diag)
+      call ezfio_get_davidson_keywords_n_states_diag(n_states_diag)
     else
-      print *, 'davidson/n_states_diag not found in EZFIO file'
+      print *, 'davidson_keywords/n_states_diag not found in EZFIO file'
       stop 1
     endif
     n_states_diag = max(2,N_states * N_states_diag)
@@ -37,3 +40,4 @@ BEGIN_PROVIDER [ integer, n_states_diag  ]
 
 END_PROVIDER
 
+! ---

src/davidson_keywords/usef.irp.f

@@ -0,0 +1,33 @@
+use bitmasks
+use f77_zmq
+
+
+! ---
+
+BEGIN_PROVIDER [ integer, nthreads_davidson ]
+ implicit none
+ BEGIN_DOC
+ ! Number of threads for Davidson
+ END_DOC
+ nthreads_davidson = nproc
+ character*(32) :: env
+ call getenv('QP_NTHREADS_DAVIDSON',env)
+ if (trim(env) /= '') then
+   read(env,*) nthreads_davidson
+   call write_int(6,nthreads_davidson,'Target number of threads for <Psi|H|Psi>')
+ endif
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, threshold_davidson_pt2 ]
+ implicit none
+ BEGIN_DOC
+ ! Threshold of Davidson's algorithm, using PT2 as a guide
+ END_DOC
+ threshold_davidson_pt2 = threshold_davidson
+
+END_PROVIDER
+
+! ---
+

src/davidson_undressed/null_dressing_vector.irp.f

@@ -1,10 +1,12 @@
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
+&BEGIN_PROVIDER [ double precision, dressing_delta   , (N_det,N_states) ]
  implicit none
  BEGIN_DOC
  ! Null dressing vectors
  END_DOC
  dressing_column_h(:,:) = 0.d0
  dressing_column_s(:,:) = 0.d0
+ dressing_delta   (:,:) = 0.d0
 END_PROVIDER
 

src/determinants/spindeterminants.ezfio_config

@@ -9,8 +9,11 @@ spindeterminants
   psi_det_beta    integer*8 (spindeterminants_n_int*spindeterminants_bit_kind/8,spindeterminants_n_det_beta)
   psi_coef_matrix_rows    integer (spindeterminants_n_det)
   psi_coef_matrix_columns integer (spindeterminants_n_det)
-  psi_coef_matrix_values  double precision (spindeterminants_n_det,spindeterminants_n_states)
+  psi_coef_matrix_values      double precision (spindeterminants_n_det,spindeterminants_n_states)
+  psi_left_coef_matrix_values double precision (spindeterminants_n_det,spindeterminants_n_states)
   n_svd_coefs     integer
+  n_svd_alpha     integer
+  n_svd_beta      integer
   psi_svd_alpha   double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs,spindeterminants_n_states)
   psi_svd_beta    double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs,spindeterminants_n_states)
   psi_svd_coefs   double precision (spindeterminants_n_svd_coefs,spindeterminants_n_states)