From 35509fb65f57489d7afee50cb6a19ade72333e9d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 30 Sep 2016 22:29:22 +0200 Subject: [PATCH] Added S2 davidson to MRCC --- plugins/MRCC_Utils/davidson.irp.f | 618 ++++++++++++++++++++++++++++ plugins/MRCC_Utils/mrcc_utils.irp.f | 73 +--- 2 files changed, 626 insertions(+), 65 deletions(-) diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 085799f6..a67ca676 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -548,3 +548,621 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8) end + +subroutine davidson_diag_mrcc_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit,istate) + use bitmasks + implicit none + BEGIN_DOC + ! 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 + ! + ! iunit : Unit number for the I/O + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit, istate + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st_diag) + double precision, intent(out) :: energies(N_st) + double precision, allocatable :: H_jj(:), S2_jj(:) + + double precision :: diag_h_mat_elem + integer :: i + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + PROVIDE mo_bielec_integrals_in_map + allocate(H_jj(sze), S2_jj(sze)) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(sze,H_jj,S2_jj, dets_in,Nint,N_det_ref,delta_ii, & + !$OMP idx_ref, istate) & + !$OMP PRIVATE(i) + !$OMP DO SCHEDULE(guided) + do i=1,sze + H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + call get_s2(dets_in(1,1,i),dets_in(1,1,i),Nint,S2_jj(i)) + enddo + !$OMP END DO + !$OMP DO SCHEDULE(guided) + do i=1,N_det_ref + H_jj(idx_ref(i)) += delta_ii(istate,i) + enddo + !$OMP END DO + !$OMP END PARALLEL + + call davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) + deallocate (H_jj,S2_jj) +end + + +subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate ) + use bitmasks + implicit none + BEGIN_DOC + ! Davidson diagonalization with specific diagonal elements of the H matrix + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! S2_jj : specific diagonal S^2 matrix elements + ! + ! 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 : Number of states in which H is diagonalized. Assumed > sze + ! + ! iunit : Unit for the I/O + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, istate + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(in) :: H_jj(sze), S2_jj(sze) + integer, intent(in) :: iunit + double precision, intent(inout) :: u_in(dim_in,N_st_diag) + double precision, intent(out) :: energies(N_st_diag) + + integer :: sze_8 + integer :: iter + integer :: i,j,k,l,m + logical :: converged + + double precision, allocatable :: overlap(:,:) + double precision :: u_dot_v, u_dot_u + + integer, allocatable :: kl_pairs(:,:) + integer :: k_pairs, kl + + integer :: iter2 + double precision, allocatable :: W(:,:), U(:,:), R(:,:), S(:,:) + double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:) + double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:) + double precision :: diag_h_mat_elem + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(3,N_st) + double precision :: cpu, wall + integer :: shift, shift2 + include 'constants.include.F' + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, S, y, h, lambda + if (N_st_diag > sze) then + stop 'error in Davidson : N_st_diag > sze' + endif + + PROVIDE nuclear_repulsion + + call write_time(iunit) + call wall_time(wall) + call cpu_time(cpu) + write(iunit,'(A)') '' + write(iunit,'(A)') 'Davidson Diagonalization' + write(iunit,'(A)') '------------------------' + write(iunit,'(A)') '' + call write_int(iunit,N_st,'Number of states') + call write_int(iunit,N_st_diag,'Number of states in diagonalization') + call write_int(iunit,sze,'Number of determinants') + call write_int(iunit,istate,'Using dressing for state ') + + write(iunit,'(A)') '' + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ =========== ===========' + enddo + write(iunit,'(A)') trim(write_buffer) + write_buffer = ' Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy S^2 Residual' + enddo + write(iunit,'(A)') trim(write_buffer) + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ =========== ===========' + enddo + write(iunit,'(A)') trim(write_buffer) + + integer, external :: align_double + sze_8 = align_double(sze) + + double precision :: delta + + if (s2_eig) then + delta = 1.d0 + else + delta = 0.d0 + endif + + allocate( & + kl_pairs(2,N_st_diag*(N_st_diag+1)/2), & + W(sze_8,N_st_diag*davidson_sze_max), & + U(sze_8,N_st_diag*davidson_sze_max), & + R(sze_8,N_st_diag), & + S(sze_8,N_st_diag*davidson_sze_max), & + h(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & + y(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & + s_(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & + s_tmp(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & + residual_norm(N_st_diag), & + overlap(N_st_diag,N_st_diag), & + c(N_st_diag*davidson_sze_max), & + s2(N_st_diag*davidson_sze_max), & + lambda(N_st_diag*davidson_sze_max)) + + 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=1,N_st + call normalize(u_in(1,k),sze) + enddo + + do k=N_st+1,N_st_diag + do i=1,sze + double precision :: r1, r2 + call random_number(r1) + call random_number(r2) + u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2) + enddo + + ! Gram-Schmidt + ! ------------ + 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), & + c,1,1.d0,u_in(1,k),1) + call normalize(u_in(1,k),sze) + enddo + + + do while (.not.converged) + + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + do iter=1,davidson_sze_max-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------------- + + + call H_S2_u_0_mrcc_nstates(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,& + istate,N_st_diag,sze_8) + + + ! Compute h_kl = = + ! ------------------------------------------- + + +! do l=1,N_st_diag +! do k=1,N_st_diag +! do iter2=1,iter-1 +! h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) +! h(k,iter,l,iter2) = h(k,iter2,l,iter) +! enddo +! enddo +! do k=1,l +! h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) +! h(l,iter,k,iter) = h(k,iter,l,iter) +! enddo +! enddo + + call dgemm('T','N', shift2, N_st_diag, sze, & + 1.d0, U, size(U,1), W(1,shift+1), size(W,1), & + 0.d0, h(1,shift+1), size(h,1)) + + call dgemm('T','N', shift2, N_st_diag, sze, & + 1.d0, U, size(U,1), S(1,shift+1), size(S,1), & + 0.d0, s_(1,shift+1), size(s_,1)) + + ! Diagonalize h + ! ------------- + call lapack_diag(lambda,y,h,size(h,1),shift2) + + ! Compute S2 for each eigenvector + ! ------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, s_, size(s_,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, s_, size(s_,1)) + + do k=1,shift2 + s2(k) = s_(k,k) + S_z2_Sz + enddo + + if (s2_eig) then + logical :: state_ok(N_st_diag*davidson_sze_max) + do k=1,shift2 + state_ok(k) = (dabs(s2(k)-expected_s2) < 0.3d0) + enddo + do k=1,shift2 + if (.not. state_ok(k)) then + do l=k+1,shift2 + if (state_ok(l)) then + call dswap(shift2, y(1,k), 1, y(1,l), 1) + call dswap(1, s2(k), 1, s2(l), 1) + call dswap(1, lambda(k), 1, lambda(l), 1) + state_ok(k) = .True. + state_ok(l) = .False. + exit + endif + enddo + endif + enddo + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + +! do k=1,N_st_diag +! do i=1,sze +! U(i,shift2+k) = 0.d0 +! W(i,shift2+k) = 0.d0 +! S(i,shift2+k) = 0.d0 +! enddo +! do l=1,N_st_diag*iter +! do i=1,sze +! U(i,shift2+k) = U(i,shift2+k) + U(i,l)*y(l,k) +! W(i,shift2+k) = W(i,shift2+k) + W(i,l)*y(l,k) +! S(i,shift2+k) = S(i,shift2+k) + S(i,l)*y(l,k) +! enddo +! enddo +! enddo +! +! + 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)) + 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)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, S, size(S,1), y, size(y,1), 0.d0, S(1,shift2+1), size(S,1)) + + ! Compute residual vector + ! ----------------------- + +! do k=1,N_st_diag +! print *, s2(k) +! s2(k) = u_dot_v(U(1,shift2+k), S(1,shift2+k), sze) + S_z2_Sz +! print *, s2(k) +! print *, '' +! pause +! enddo + do k=1,N_st_diag + do i=1,sze + R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz) + enddo + if (k <= N_st) then + residual_norm(k) = u_dot_u(R(1,k),sze) + to_print(1,k) = lambda(k) + nuclear_repulsion + to_print(2,k) = s2(k) + to_print(3,k) = residual_norm(k) + if (residual_norm(k) > 1.e9) then + stop 'Davidson failed' + endif + endif + enddo + + write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(:,1:N_st) + call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) + if (converged) then + exit + endif + + ! Davidson step + ! ------------- + + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = - R(i,k)/max(H_jj(i) - lambda(k),1.d-2) + enddo + enddo + + ! Gram-Schmidt + ! ------------ + + do k=1,N_st_diag + +! do l=1,N_st_diag*iter +! c(1) = u_dot_v(U(1,shift2+k),U(1,l),sze) +! do i=1,sze +! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,l) +! enddo +! enddo +! + call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & + U(1,shift2+k),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,shift2+k),1) +! +! do l=1,k-1 +! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze) +! do i=1,sze +! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,shift2+l) +! enddo +! enddo +! + call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), & + U(1,shift2+k),1,0.d0,c,1) + call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), & + c,1,1.d0,U(1,shift2+k),1) + + call normalize( U(1,shift2+k), sze ) + enddo + + enddo + + if (.not.converged) then + iter = davidson_sze_max-1 + endif + + ! Re-contract to u_in + ! ----------- + + do k=1,N_st_diag + energies(k) = lambda(k) + enddo + +! do k=1,N_st_diag +! do i=1,sze +! do l=1,iter*N_st_diag +! u_in(i,k) += U(i,l)*y(l,k) +! enddo +! 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), 0.d0, u_in, size(u_in,1)) + + enddo + + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ =========== ===========' + enddo + write(iunit,'(A)') trim(write_buffer) + write(iunit,'(A)') '' + call write_time(iunit) + + deallocate ( & + kl_pairs, & + W, residual_norm, & + U, overlap, & + R, c, S, & + h, & + y, s_, s_tmp, & + lambda & + ) +end + + +subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + ! + ! S2_jj : array of + END_DOC + integer, intent(in) :: N_st,n,Nint, sze_8, istate_in + double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st) + double precision, intent(in) :: u_0(sze_8,N_st) + double precision, intent(in) :: H_jj(n), S2_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + double precision :: hij,s2 + double precision, allocatable :: vt(:,:), ut(:,:), st(:,:) + integer :: i,j,k,l, jj,ii + integer :: i0, j0 + + integer, allocatable :: shortcut(:,:), sort_idx(:,:) + integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) + integer(bit_kind) :: sorted_i(Nint) + + integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate + integer :: N_st_8 + + integer, external :: align_double + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut + + N_st_8 = align_double(N_st) + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (n>0) + PROVIDE ref_bitmask_energy + + allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) + allocate(ut(N_st_8,n)) + + v_0 = 0.d0 + s_0 = 0.d0 + + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(i,istate) + enddo + enddo + + call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) + call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& + !$OMP SHARED(n,keys_tmp,ut,Nint,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8, & + !$OMP N_det_ref, idx_ref, N_det_non_ref, idx_non_ref, delta_ij,istate_in) + allocate(vt(N_st_8,n),st(N_st_8,n)) + Vt = 0.d0 + St = 0.d0 + + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,1) + do sh2=sh,shortcut(0,1) + exa = 0 + do ni=1,Nint + exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) + end do + if(exa > 2) then + cycle + end if + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + if(sh==sh2) then + endi = i-1 + else + endi = shortcut(sh2+1,1)-1 + end if + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh2,1),endi + org_j = sort_idx(j,1) + ext = exa + do ni=1,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + end do + if(ext <= 4) then + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) + vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) + st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) + enddo + endif + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,2) + do i=shortcut(sh,2),shortcut(sh+1,2)-1 + org_i = sort_idx(i,2) + do j=shortcut(sh,2),i-1 + org_j = sort_idx(j,2) + ext = 0 + do ni=1,Nint + ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) + end do + if(ext == 4) then + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) + vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) + st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) + enddo + end if + end do + end do + enddo + !$OMP END DO NOWAIT + +! -------------------------- +! Begin Specific to dressing +! -------------------------- + + !$OMP DO + do ii=1,n_det_ref + i = idx_ref(ii) + do jj = 1, n_det_non_ref + j = idx_non_ref(jj) + do istate=1,N_st + vt (istate,i) = vt (istate,i) + delta_ij(istate_in,jj,ii)*ut(istate,j) + vt (istate,j) = vt (istate,j) + delta_ij(istate_in,jj,ii)*ut(istate,i) + enddo + enddo + enddo + !$OMP END DO + +! ------------------------ +! End Specific to dressing +! ------------------------ + + !$OMP CRITICAL + do istate=1,N_st + do i=n,1,-1 + v_0(i,istate) = v_0(i,istate) + vt(istate,i) + s_0(i,istate) = s_0(i,istate) + st(istate,i) + enddo + enddo + !$OMP END CRITICAL + + deallocate(vt,st) + !$OMP END PARALLEL + + do istate=1,N_st + do i=1,n + v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate) + s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate) + enddo + enddo + deallocate (shortcut, sort_idx, sorted, version, ut) +end + diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 12ae089b..fcfb5dd1 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -148,8 +148,14 @@ END_PROVIDER if (diag_algorithm == "Davidson") then - call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,& - size(CI_eigenvectors_dressed,1),N_det,N_states,N_states_diag,N_int,output_determinants,mrcc_state) +! call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,& +! size(CI_eigenvectors_dressed,1),N_det,N_states,N_states_diag,N_int,output_determinants,mrcc_state) + + call davidson_diag_mrcc_HS2(psi_det,CI_eigenvectors_dressed,& + size(CI_eigenvectors_dressed,1), & + CI_electronic_energy_dressed,N_det,N_states,N_states_diag,N_int, & + output_determinants,mrcc_state) + call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& N_states_diag,size(CI_eigenvectors_dressed,1)) @@ -236,69 +242,6 @@ END_PROVIDER deallocate(eigenvectors,eigenvalues) endif - if( s2_eig.and.(n_states_diag > 1).and.(n_det >= n_states_diag) )then - ! Diagonalizing S^2 within the "n_states_diag" states found - allocate(s2_eigvalues(N_states_diag), e_array(N_states_diag)) - call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors_dressed,n_det,size(psi_det,3),size(CI_eigenvectors_dressed,1),min(n_states_diag,n_det),s2_eigvalues) - - double precision, allocatable :: psi_coef_tmp(:,:) - allocate(psi_coef_tmp(psi_det_size,N_states_diag)) - do j = 1, N_states_diag - do i = 1, N_det - psi_coef_tmp(i,j) = CI_eigenvectors_dressed(i,j) - enddo - enddo - call u_0_H_u_0_mrcc_nstates(e_array,psi_coef_tmp,n_det,psi_det,N_int,mrcc_state,N_states,psi_det_size) - - ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value - ! closer to the "expected_s2" set as input - - allocate(index_good_state_array(N_det),good_state_array(N_det)) - good_state_array = .False. - i_state = 0 - do j = 1, N_states_diag - if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then - good_state_array(j) = .True. - i_state +=1 - index_good_state_array(i_state) = j - endif - enddo - ! Sorting the i_state good states by energy - allocate(iorder(i_state)) - do j = 1, i_state - do i = 1, N_det - CI_eigenvectors_dressed(i,j) = psi_coef_tmp(i,index_good_state_array(j)) - enddo - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j)) - CI_electronic_energy_dressed(j) = e_array(j) - iorder(j) = j - enddo - call dsort(e_array,iorder,i_state) - do j = 1, i_state - CI_electronic_energy_dressed(j) = e_array(j) - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(iorder(j))) - do i = 1, N_det - CI_eigenvectors_dressed(i,j) = psi_coef_tmp(i,index_good_state_array(iorder(j))) - enddo - enddo - - ! Then setting the other states without any specific energy order - i_other_state = 0 - do j = 1, N_states_diag - if(good_state_array(j))cycle - i_other_state +=1 - do i = 1, N_det - CI_eigenvectors_dressed(i,i_state + i_other_state) = psi_coef_tmp(i,j) - enddo - CI_eigenvectors_s2_dressed(i_state + i_other_state) = s2_eigvalues(j) - CI_electronic_energy_dressed(i_state + i_other_state) = e_array(i_state + i_other_state) - enddo - deallocate(iorder,e_array,index_good_state_array,good_state_array,psi_coef_tmp) - - deallocate(s2_eigvalues) - - endif - END_PROVIDER BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]