diff --git a/plugins/MRCC_CASSD/EZFIO.cfg b/plugins/MRCC_CASSD/EZFIO.cfg index e145c9e0..17ee7f29 100644 --- a/plugins/MRCC_CASSD/EZFIO.cfg +++ b/plugins/MRCC_CASSD/EZFIO.cfg @@ -7,11 +7,11 @@ interface: ezfio type: Threshold doc: Threshold on the convergence of the MRCC energy interface: ezfio,provider,ocaml -default: 1.e-7 +default: 1.e-5 [n_it_mrcc_max] type: Strictly_positive_int doc: Maximum number of MRCC iterations interface: ezfio,provider,ocaml -default: 20 +default: 10 diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f index e784a167..38cd3c55 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -1,13 +1,100 @@ program mrcc implicit none - if (.not.read_wf) then - print *, 'read_wf has to be true.' - stop 1 - endif + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + + read_wf = .True. + SOFT_TOUCH read_wf call print_cas_coefs - call run_mrcc + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end +subroutine run(N_st,energy) + implicit none + + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + + integer :: i + + double precision :: E_new, E_old, delta_e + integer :: iteration + double precision :: E_past(4), lambda + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + lambda = 1.d0 + do while (delta_E > thresh_mrcc) + iteration += 1 + print *, '===========================' + print *, 'MRCC Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCC energy") + call diagonalize_ci_dressed(lambda) + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + call save_wavefunction + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + if (iteration > n_it_mrcc_max) then + exit + endif + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + energy(:) = ci_energy_dressed(:) + +end + + +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + N_det_generators = lambda_mrcc_pt2(0) + do i=1,N_det_generators + j = lambda_mrcc_pt2(i) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + + + call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + + + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) + +end + + subroutine print_cas_coefs implicit none @@ -18,7 +105,7 @@ subroutine print_cas_coefs print *, psi_cas_coef(i,:) call debug_det(psi_cas(1,1,i),N_int) enddo - call write_double(6,ci_energy(1),"Initial CI energy") + end diff --git a/plugins/MRCC_CASSD/mrcc_noiter.irp.f b/plugins/MRCC_CASSD/mrcc_noiter.irp.f new file mode 100644 index 00000000..8d95cea9 --- /dev/null +++ b/plugins/MRCC_CASSD/mrcc_noiter.irp.f @@ -0,0 +1,91 @@ +program mrcc_noiter + implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + read_wf = .True. + threshold_generators = .9999d0 + SOFT_TOUCH read_wf threshold_generators + call print_cas_coefs + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) +end + +subroutine run(N_st,energy) + implicit none + + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef ci_energy_dressed + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + energy(:) = ci_energy_dressed(:) +end + + +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + N_det_generators = lambda_mrcc_pt2(0) + do i=1,N_det_generators + j = lambda_mrcc_pt2(i) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + + + call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) + +end + + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + call write_double(6,ci_energy(1),"Initial CI energy") + +end + diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index decc5a75..df2b67a0 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -24,5 +24,12 @@ s.data["size_max"] = "3072" print s + +s = H_apply_zmq("mrcc_PT2") +s.energy = "ci_electronic_energy_dressed" +s.set_perturbation("epstein_nesbet_2x2") +s.unset_openmp() +print s + END_SHELL diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index d278ba13..354ea389 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -189,10 +189,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin ! Davidson iterations ! =================== - converged = .False. + integer :: iteration + converged = .False. do while (.not.converged) - !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) do k=1,N_st @@ -206,6 +206,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin do iter=1,davidson_sze_max-1 + ! Compute W_k = H |u_k> ! ---------------------- diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index c3f3debf..b2304818 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -24,7 +24,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) - integer :: i,j,k,l + integer :: i,j,k,l,m integer :: degree_alpha(psi_det_size) integer :: idx_alpha(0:psi_det_size) logical :: good, fullMatch @@ -48,9 +48,14 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge integer :: N_miniList, ni, leng double precision, allocatable :: hij_cache(:) + integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:) + integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) + integer :: mobiles(2), smallerlist + + leng = max(N_det_generators, N_det_non_ref) - allocate(miniList(Nint, 2, leng), idx_miniList(leng), hij_cache(N_det_non_ref)) + allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref)) !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) @@ -59,8 +64,21 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge return end if + allocate(ptr_microlist(0:mo_tot_num*2+1), & + N_microlist(0:mo_tot_num*2) ) + allocate( microlist(Nint,2,N_minilist*4), & + idx_microlist(N_minilist*4)) - call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + if(key_mask(1,1) /= 0) then + call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) + call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + else + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + end if + + + + deallocate(microlist, idx_microlist) allocate (dIa_hla(Nstates,Ndet_non_ref)) @@ -69,17 +87,101 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge ! |alpha> if(N_tq > 0) then - call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint) + call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint) + if(N_minilist == 0) return + + + if(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!! + allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist)) + + allocate( microlist(Nint,2,N_minilist*4), & + idx_microlist(N_minilist*4)) + call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) + + + do i=0,mo_tot_num*2 + do k=ptr_microlist(i),ptr_microlist(i+1)-1 + idx_microlist(k) = idx_minilist(idx_microlist(k)) + end do + end do + + do l=1,N_microlist(0) + do k=1,Nint + microlist_zero(k,1,l) = microlist(k,1,l) + microlist_zero(k,2,l) = microlist(k,2,l) + enddo + idx_microlist_zero(l) = idx_microlist(l) + enddo + end if end if + + do i_alpha=1,N_tq - ! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) - call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + if(key_mask(1,1) /= 0) then + call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) + + if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then + smallerlist = mobiles(1) + else + smallerlist = mobiles(2) + end if + - do j=1,idx_alpha(0) - idx_alpha(j) = idx_miniList(idx_alpha(j)) - end do + do l=0,N_microlist(smallerlist)-1 + microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l) + idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l) + end do + + call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_microlist_zero(idx_alpha(j)) + end do + + +! i = 1 +! j = 2 +! do j = 2, idx_alpha_tmp(0) +! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit +! end do +! +! m = j +! +! idx_alpha(0) = idx_alpha_tmp(0) +! +! do l = 1, idx_alpha(0) +! if(j > idx_alpha_tmp(0)) then +! k = i +! i += 1 +! else if(i >= m) then +! k = j +! j += 1 +! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then +! k = i +! i += 1 +! else +! k = j +! j += 1 +! end if +! ! k=l +! idx_alpha(l) = idx_alpha_tmp(k) +! degree_alpha(l) = degree_alpha_tmp(k) +! end do +! + else + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + end do + end if + + +! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) +! do j=1,idx_alpha(0) +! idx_alpha(j) = idx_miniList(idx_alpha(j)) +! end do + !print *, idx_alpha(:idx_alpha(0)) do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) @@ -133,32 +235,17 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge tmp_det(k,1) = psi_ref(k,1,i_I) tmp_det(k,2) = psi_ref(k,2,i_I) enddo - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree_alpha(k_sd) == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - endif + logical :: ok + call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) + if(.not. ok) cycle ! do i_state=1,Nstates dka(i_state) = 0.d0 enddo do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) if (degree == 0) then @@ -172,11 +259,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge if (.not.loop) then call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) hIl = hij_mrcc(idx_alpha(l_sd),i_I) - ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) do i_state=1,Nstates dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 enddo endif + exit endif enddo @@ -211,21 +299,13 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo - deallocate (dIa_hla,hij_cache) - deallocate(miniList, idx_miniList) + !deallocate (dIa_hla,hij_cache) + !deallocate(miniList, idx_miniList) end - BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ] - gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators) - gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators) - call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int) - call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int) -END_PROVIDER + subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) @@ -258,6 +338,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq N_tq = 0 + i_loop : do i=1,N_selected if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then cycle @@ -287,8 +368,84 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq end +subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + + + integer :: nt,ni + logical, external :: is_connected_to + + + integer(bit_kind),intent(in) :: microlist(Nint,2,*) + integer,intent(in) :: ptr_microlist(0:*) + integer,intent(in) :: N_microlist(0:*) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + + integer :: mobiles(2), smallerlist + + N_tq = 0 + + + + i_loop : do i=1,N_selected + call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint) + if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then + smallerlist = mobiles(1) + else + smallerlist = mobiles(2) + end if + + if(N_microlist(smallerlist) > 0) then + if(is_connected_to(det_buffer(1,1,i), microlist(1,1,ptr_microlist(smallerlist)), Nint, N_microlist(smallerlist))) then + cycle + end if + end if + + if(N_microlist(0) > 0) then + if(is_connected_to(det_buffer(1,1,i), microlist, Nint, N_microlist(0))) then + cycle + end if + end if + + ! Select determinants that are triple or quadruple excitations + ! from the ref + good = .True. + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo i_loop +end + + - diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index 5d9acfc1..50343fdb 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -11,12 +11,13 @@ subroutine mrcc_iterations double precision :: E_new, E_old, delta_e integer :: iteration,i_oscillations - double precision :: E_past(4) + double precision :: E_past(4), lambda E_new = 0.d0 delta_E = 1.d0 iteration = 0 j = 1 i_oscillations = 0 + lambda = 1.d0 do while (delta_E > 1.d-7) iteration += 1 print *, '===========================' @@ -25,10 +26,15 @@ subroutine mrcc_iterations print *, '' E_old = sum(ci_energy_dressed) call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed + call diagonalize_ci_dressed(lambda) E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) - +! if (E_new > E_old) then +! lambda = lambda * 0.7d0 +! else +! lambda = min(1.d0, lambda * 1.1d0) +! endif +! print *, 'energy lambda ', lambda E_past(j) = E_new j +=1 call save_wavefunction diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index c41a3431..e6c84a06 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,45 +1,54 @@ -BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] - implicit none - BEGIN_DOC - ! cm/ or perturbative 1/Delta_E(m) - END_DOC - integer :: i,k - double precision :: ihpsi(N_states),ihpsi_current(N_states) - integer :: i_pert_count - - i_pert_count = 0 - lambda_mrcc = 0.d0 - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 - endif - if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then - i_pert_count +=1 - else - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - endif - enddo - enddo + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] +&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] + implicit none + BEGIN_DOC + ! cm/ or perturbative 1/Delta_E(m) + END_DOC + integer :: i,k + double precision :: ihpsi_current(N_states) + integer :: i_pert_count + double precision :: hii, lambda_pert + integer :: N_lambda_mrcc_pt2 + + i_pert_count = 0 + lambda_mrcc = 0.d0 + N_lambda_mrcc_pt2 = 0 + lambda_mrcc_pt2(0) = 0 + + do i=1,N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_normalized, N_int, N_det_ref,& + size(psi_ref_coef,1), N_states,ihpsi_current) + call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif + lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) ) + lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 + if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then + N_lambda_mrcc_pt2 += 1 + lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i + endif + endif + enddo + enddo + lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 + + print*,'N_det_non_ref = ',N_det_non_ref + print*,'Number of ignored determinants = ',i_pert_count + print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) + print*,'lambda max = ',maxval(dabs(lambda_mrcc)) - print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of ignored determinants = ',i_pert_count - print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) END_PROVIDER -!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] -!implicit none -!BEGIN_DOC -!! Dressing matrix in SD basis -!END_DOC -!delta_ij_non_ref = 0.d0 -!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) -!END_PROVIDER + + BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] implicit none @@ -174,17 +183,19 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] END_PROVIDER -subroutine diagonalize_CI_dressed +subroutine diagonalize_CI_dressed(lambda) implicit none BEGIN_DOC ! Replace the coefficients of the CI states by the coefficients of the ! eigenstates of the CI matrix END_DOC + double precision, intent(in) :: lambda integer :: i,j do j=1,N_states_diag do i=1,N_det - psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j) enddo + call normalize(psi_coef(1,j), N_det) enddo SOFT_TOUCH psi_coef diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 94b6b8b0..7bb08c21 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -3,7 +3,7 @@ import perturbation END_SHELL -subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) +subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -14,6 +14,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) integer(bit_kind),intent(in) :: key_mask(Nint,2) double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) + double precision, intent(in) :: electronic_energy(N_st) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) @@ -151,7 +152,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c idx_microlist_zero(ptr_microlist(1)+l) = idx_microlist(ptr_microlist(smallerlist)+l) enddo end if - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_microlist(smallerlist)+N_microlist(0), & n_st,microlist_zero,idx_microlist_zero,N_microlist(smallerlist)+N_microlist(0)) else @@ -160,11 +161,11 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c cycle end if - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) end if -! call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & +! call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & ! c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st @@ -182,7 +183,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c end -subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) +subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -193,6 +194,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) integer(bit_kind),intent(in) :: key_mask(Nint,2) double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) + double precision, intent(in) :: electronic_energy(N_st) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) @@ -241,7 +243,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ cycle endif - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 72d03808..e990a37c 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -29,11 +29,11 @@ subroutine pt2_epstein_nesbet ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st - if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then + if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then c_pert(i) = -1.d0 e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) + else if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h) H_pert_diag(i) = h*c_pert(i)*c_pert(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else @@ -66,7 +66,6 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) double precision :: i_H_psi_array(N_st) ASSERT (Nint == N_int) ASSERT (Nint > 0) - PROVIDE CI_electronic_energy !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) @@ -74,7 +73,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st if (i_H_psi_array(i) /= 0.d0) then - delta_e = h - CI_electronic_energy(i) + delta_e = h - electronic_energy(i) if (delta_e > 0.d0) then e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) else @@ -165,7 +164,7 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments) ! ! that can be repeated by repeating all the double excitations ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! : you repeat all the correlation energy already taken into account in electronic_energy(1) ! ! that could be repeated to this determinant. ! @@ -195,16 +194,16 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments) enddo h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) h = h + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + delta_E = 1.d0/(electronic_energy(1) - h) - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h) e_2_pert(1) = i_H_psi_array(1) * c_pert(1) do i =2,N_st H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h)) e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) else c_pert(i) = i_H_psi_array(i) @@ -248,7 +247,7 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) ! ! that can be repeated by repeating all the double excitations ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! : you repeat all the correlation energy already taken into account in electronic_energy(1) ! ! that could be repeated to this determinant. ! @@ -278,16 +277,16 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) enddo h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) h = h + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + delta_E = 1.d0/(electronic_energy(1) - h) - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h) e_2_pert(1) = i_H_psi_array(1) * c_pert(1) do i =2,N_st H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h)) e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) else c_pert(i) = i_H_psi_array(i) @@ -328,11 +327,11 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st - if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then + if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then c_pert(i) = -1.d0 e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) + else if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h) H_pert_diag(i) = h*c_pert(i)*c_pert(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else @@ -348,7 +347,7 @@ end SUBST [ arguments, declarations ] -det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; +electronic_energy,det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; integer, intent(in) :: Nint integer, intent(in) :: ndet @@ -357,6 +356,7 @@ det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minili integer(bit_kind), intent(in) :: det_ref (Nint,2) integer(bit_kind), intent(in) :: det_pert(Nint,2) double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + double precision , intent(in) :: electronic_energy(N_st) double precision , intent(out) :: c_pert(N_st) double precision , intent(out) :: e_2_pert(N_st) double precision, intent(out) :: H_pert_diag(N_st) diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index f67f0587..c7b9867f 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -24,6 +24,21 @@ use bitmasks enddo enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized, (psi_det_size,n_states) ] + implicit none + BEGIN_DOC +! Normalized coefficients of the reference + END_DOC + integer :: i,j,k + do k=1,N_states + do j=1,N_det_ref + psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k) + enddo + call normalize(psi_ref_coef_normalized(1,k), N_det_ref) + enddo + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ] diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index c6466569..436f092d 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -61,6 +61,7 @@ class H_apply(object): s["params_post"] = "" self.selection_pt2 = None + self.energy = "CI_electronic_energy" self.perturbation = None self.do_double_exc = do_double_exc #s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) & @@ -264,13 +265,13 @@ class H_apply(object): self.data["keys_work"] = """ ! if(check_double_excitation)then call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) - """%(pert) + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s) + """%(pert,self.energy) else: self.data["keys_work"] = """ call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) - """%(pert) + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s) + """%(pert,self.energy) self.data["finalization"] = """ diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index b63ae69f..4ab9deda 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1728,3 +1728,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) deallocate (shortcut, sort_idx, sorted, version) end + +subroutine apply_excitation(det, exc, res, ok, Nint) + use bitmasks + implicit none + + integer, intent(in) :: Nint + integer, intent(in) :: exc(0:2,2,2) + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: h1,p1,h2,p2,s1,s2,degree + integer :: ii, pos + + + ok = .false. + degree = exc(0,1,1) + exc(0,1,2) + + if(.not. (degree > 0 .and. degree <= 2)) then + print *, degree + print *, "apply ex" + STOP + endif + + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + res = det + + ii = (h1-1)/bit_kind_size + 1 + pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s1) = ibclr(res(ii, s1), pos) + + ii = (p1-1)/bit_kind_size + 1 + pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) + + if(degree == 2) then + ii = (h2-1)/bit_kind_size + 1 + pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s2) = ibclr(res(ii, s2), pos) + + ii = (p2-1)/bit_kind_size + 1 + pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s2) = ibset(res(ii, s2), pos) + endif + + ok = .true. +end subroutine + + diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 1be43412..1ced9e1d 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -53,10 +53,10 @@ function test_exe() { } function run_HF() { - thresh=1.e-8 + thresh=1.e-7 test_exe SCF || skip ezfio set_file $1 - ezfio set hartree_fock thresh_scf 1.e-10 + ezfio set hartree_fock thresh_scf 1.e-11 qp_run SCF $1 energy="$(ezfio get hartree_fock energy)" eq $energy $2 $thresh @@ -155,7 +155,7 @@ function run_all_1h_1p() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2289109271715 1.E-3 + eq $energy -76.2284994316618 1.e-4 } @@ -166,7 +166,7 @@ function run_all_1h_1p() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -0.169483703904991E+02 + run_HF h2o_pseudo.ezfio -16.9483708495521 } @test "FCI H2O VDZ pseudo" {