From e2162d282d51271fa4db8b7237fbd57250b49354 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 15 Mar 2023 10:23:48 +0100 Subject: [PATCH] non_sym dress: comb --- src/davidson/EZFIO.cfg | 55 +- src/davidson/NEED | 1 + src/davidson/davidson_parallel.irp.f | 15 - .../diagonalization_hs2_dressed.irp.f | 18 +- .../diagonalization_nonsym_h_dressed.irp.f | 541 ++++++++++++++++++ src/davidson/overlap_states.irp.f | 40 ++ .../nonsym_diagonalize_ci.irp.f | 188 ++++++ src/davidson_keywords/EZFIO.cfg | 54 ++ src/davidson_keywords/README.rst | 5 + .../input.irp.f | 10 +- src/davidson_keywords/usef.irp.f | 33 ++ .../null_dressing_vector.irp.f | 2 + .../spindeterminants.ezfio_config | 5 +- ...uet => save_bitcpsileft_for_qmcchem.irp.f} | 0 14 files changed, 883 insertions(+), 84 deletions(-) create mode 100644 src/davidson/diagonalization_nonsym_h_dressed.irp.f create mode 100644 src/davidson/overlap_states.irp.f create mode 100644 src/davidson_dressed/nonsym_diagonalize_ci.irp.f create mode 100644 src/davidson_keywords/EZFIO.cfg create mode 100644 src/davidson_keywords/README.rst rename src/{davidson => davidson_keywords}/input.irp.f (79%) create mode 100644 src/davidson_keywords/usef.irp.f rename src/tc_bi_ortho/{save_bitcpsileft_for_qmcchem.irp.pouet => save_bitcpsileft_for_qmcchem.irp.f} (100%) diff --git a/src/davidson/EZFIO.cfg b/src/davidson/EZFIO.cfg index bfa55526..1152560f 100644 --- a/src/davidson/EZFIO.cfg +++ b/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 + diff --git a/src/davidson/NEED b/src/davidson/NEED index bfe31bd0..bd0abe2f 100644 --- a/src/davidson/NEED +++ b/src/davidson/NEED @@ -1 +1,2 @@ csf +davidson_keywords diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index e627dfc9..4f88506a 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/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 ') - endif -END_PROVIDER - - integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id) use f77_zmq implicit none diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index d37b7386..f4c05595 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/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 diff --git a/src/davidson/diagonalization_nonsym_h_dressed.irp.f b/src/davidson/diagonalization_nonsym_h_dressed.irp.f new file mode 100644 index 00000000..3ff060a6 --- /dev/null +++ b/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> + ! ----------------------------------- + + 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 = = + ! ------------------------------------------- + + 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 + +! --- + + + + + + + diff --git a/src/davidson/overlap_states.irp.f b/src/davidson/overlap_states.irp.f new file mode 100644 index 00000000..797d1210 --- /dev/null +++ b/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 + +! --- + diff --git a/src/davidson_dressed/nonsym_diagonalize_ci.irp.f b/src/davidson_dressed/nonsym_diagonalize_ci.irp.f new file mode 100644 index 00000000..fa4b8b33 --- /dev/null +++ b/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 + +! --- + diff --git a/src/davidson_keywords/EZFIO.cfg b/src/davidson_keywords/EZFIO.cfg new file mode 100644 index 00000000..6337b96f --- /dev/null +++ b/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 + diff --git a/src/davidson_keywords/README.rst b/src/davidson_keywords/README.rst new file mode 100644 index 00000000..9567cdb1 --- /dev/null +++ b/src/davidson_keywords/README.rst @@ -0,0 +1,5 @@ +================= +davidson_keywords +================= + +Keywords used for Davidson algorithms. diff --git a/src/davidson/input.irp.f b/src/davidson_keywords/input.irp.f similarity index 79% rename from src/davidson/input.irp.f rename to src/davidson_keywords/input.irp.f index aba88ae9..4bd79036 100644 --- a/src/davidson/input.irp.f +++ b/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 +! --- diff --git a/src/davidson_keywords/usef.irp.f b/src/davidson_keywords/usef.irp.f new file mode 100644 index 00000000..fed2ba9b --- /dev/null +++ b/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 ') + 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 + +! --- + diff --git a/src/davidson_undressed/null_dressing_vector.irp.f b/src/davidson_undressed/null_dressing_vector.irp.f index faffe964..1989bb6d 100644 --- a/src/davidson_undressed/null_dressing_vector.irp.f +++ b/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 diff --git a/src/determinants/spindeterminants.ezfio_config b/src/determinants/spindeterminants.ezfio_config index 39ccb82b..4fe1333a 100644 --- a/src/determinants/spindeterminants.ezfio_config +++ b/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) diff --git a/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.pouet b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f similarity index 100% rename from src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.pouet rename to src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f