From 45183fdd8d8606c72128d99d1c245782525be797 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 29 Nov 2016 17:54:10 +0100 Subject: [PATCH] Fixed mrcc_selected --- plugins/Psiref_threshold/psi_ref.irp.f | 13 +- plugins/mrcc_selected/dressing.irp.f | 1008 +++++++++++++++++++ plugins/mrcc_selected/dressing_slave.irp.f | 593 +++++++++++ plugins/mrcc_selected/ezfio_interface.irp.f | 61 ++ plugins/mrcc_selected/mrcc_selected.irp.f | 19 + plugins/mrcc_selected/mrcepa0_general.irp.f | 245 +++++ src/Determinants/determinants.irp.f | 8 +- 7 files changed, 1940 insertions(+), 7 deletions(-) create mode 100644 plugins/mrcc_selected/dressing.irp.f create mode 100644 plugins/mrcc_selected/dressing_slave.irp.f create mode 100644 plugins/mrcc_selected/ezfio_interface.irp.f create mode 100644 plugins/mrcc_selected/mrcc_selected.irp.f create mode 100644 plugins/mrcc_selected/mrcepa0_general.irp.f diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f index 5e722822..ee69ef5c 100644 --- a/plugins/Psiref_threshold/psi_ref.irp.f +++ b/plugins/Psiref_threshold/psi_ref.irp.f @@ -6,19 +6,22 @@ use bitmasks &BEGIN_PROVIDER [ integer, N_det_ref ] implicit none BEGIN_DOC - ! Reference wave function, defined as determinants with coefficients > 0.05 + ! Reference wave function, defined as determinants with amplitudes > 0.05 ! idx_ref gives the indice of the ref determinant in psi_det. END_DOC integer :: i, k, l logical :: good - double precision, parameter :: threshold=0.05d0 + double precision, parameter :: threshold=0.05d0 + double precision :: t(N_states) N_det_ref = 0 - t = threshold * abs_psi_coef_max + do l = 1, N_states + t(l) = threshold * abs_psi_coef_max(l) + enddo do i=1,N_det good = .False. - do l = 1, N_states + do l=1, N_states psi_ref_coef(i,l) = 0.d0 - good = good.or.(dabs(psi_coef(i,l)) > t) + good = good.or.(dabs(psi_coef(i,l)) > t(l)) enddo if (good) then N_det_ref = N_det_ref+1 diff --git a/plugins/mrcc_selected/dressing.irp.f b/plugins/mrcc_selected/dressing.irp.f new file mode 100644 index 00000000..3646b0b2 --- /dev/null +++ b/plugins/mrcc_selected/dressing.irp.f @@ -0,0 +1,1008 @@ +use bitmasks + + + + BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] + use bitmasks + implicit none + integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc + integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2) + integer(bit_kind),allocatable :: buf(:,:,:) + logical :: ok + logical, external :: detEq + + delta_ij_mrcc = 0d0 + delta_ii_mrcc = 0d0 + print *, "Dij", dij(1,1,1) + provide hh_shortcut psi_det_size! lambda_mrcc + !$OMP PARALLEL DO default(none) schedule(dynamic) & + !$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) & + !$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc) & + !$OMP private(h, n, mask, omask, buf, ok, iproc) + do gen= 1, N_det_generators + allocate(buf(N_int, 2, N_det_non_ref)) + iproc = omp_get_thread_num() + 1 + if(mod(gen, 1000) == 0) print *, "mrcc ", gen, "/", N_det_generators + do h=1, hh_shortcut(0) + call apply_hole_local(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int) + if(.not. ok) cycle + omask = 0_bit_kind + if(hh_exists(1, h) /= 0) omask = mask + n = 1 + do p=hh_shortcut(h), hh_shortcut(h+1)-1 + call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int) + if(ok) n = n + 1 + if(n > N_det_non_ref) stop "MRCC..." + end do + n = n - 1 + + if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask) + + end do + deallocate(buf) + end do + !$OMP END PARALLEL DO +END_PROVIDER + + +! subroutine blit(b1, b2) +! double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref) +! b1 = b1 + b2 +! end subroutine + + +subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) + double precision, intent(inout) :: delta_ii_(N_states,N_det_ref) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l,m + integer,allocatable :: idx_alpha(:), degree_alpha(:) + logical :: good, fullMatch + + integer(bit_kind),allocatable :: tq(:,:,:) + integer :: N_tq, c_ref ,degree + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hla(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer :: iint, ipos + integer :: i_state, k_sd, l_sd, i_I, i_alpha + + integer(bit_kind),allocatable :: miniList(:,:,:) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer,allocatable :: idx_miniList(:) + 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 + logical, external :: detEq, is_generable + !double precision, external :: get_dij, get_dij_index + + + leng = max(N_det_generators, N_det_non_ref) + allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref)) + allocate(idx_alpha(0:psi_det_size), degree_alpha(psi_det_size)) + !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) + +! if(fullMatch) then +! 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)) + + if(key_mask(1,1) /= 0) then + call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) + call filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + else + call filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + end if + + + + deallocate(microlist, idx_microlist) + + allocate (dIa_hla(N_states,N_det_non_ref)) + + ! |I> + + ! |alpha> + + if(N_tq > 0) then + 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 + 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 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 + + 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 + + + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd)) + enddo + ! |I> + do i_I=1,N_det_ref + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint) + if (degree > 4) then + cycle + endif + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + + ! |alpha> + do k_sd=1,idx_alpha(0) + ! Loop if lambda == 0 + logical :: loop +! loop = .True. +! do i_state=1,N_states +! if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then +! loop = .False. +! exit +! endif +! enddo +! if (loop) then +! cycle +! endif + + call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) + if (degree > 2) then + cycle + endif + + ! + ! + !hIk = hij_mrcc(idx_alpha(k_sd),i_I) + ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) + + do i_state=1,N_states + dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state) + !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state) + enddo + + + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do k=1,N_int + tmp_det(k,1) = psi_ref(k,1,i_I) + tmp_det(k,2) = psi_ref(k,2,i_I) + enddo + logical :: ok + call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) + if(.not. ok) cycle + + ! + do i_state=1,N_states + 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 + +! loop = .True. +! do i_state=1,N_states +! if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then +! loop = .False. +! exit +! endif +! enddo + loop = .false. + 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) + do i_state=1,N_states + dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2 + !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 + enddo + endif + + exit + endif + enddo + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,N_states + ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) + enddo + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + hla = hij_cache(k_sd) +! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla) + do i_state=1,N_states + dIa_hla(i_state,k_sd) = dIa(i_state) * hla + enddo + enddo + call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + enddo + else + delta_ii_(i_state,i_I) = 0.d0 + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd) + enddo + endif + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) + enddo + enddo + deallocate (dIa_hla,hij_cache) + deallocate(miniList, idx_miniList) +end + + + + + BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ] + use bitmasks + implicit none + integer :: i, j, i_state + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + + do i_state = 1, N_states + if(mrmode == 3) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) + end do + end do +! +! do i = 1, N_det_ref +! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) +! do j = 1, N_det_non_ref +! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) +! end do +! end do + else if(mrmode == 2) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_old(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) + end do + end do + else if(mrmode == 1) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) + end do + end do + else + stop "invalid mrmode" + end if + end do +END_PROVIDER + + +BEGIN_PROVIDER [ integer, HP, (2,N_det_non_ref) ] + integer :: i + do i=1,N_det_non_ref + call getHP(psi_non_ref(1,1,i), HP(1,i), HP(2,i), N_int) + end do +END_PROVIDER + + BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ] +&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_int,2,N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_int,2,N_det_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0, (N_int,2,N_det_non_ref) ] +&BEGIN_PROVIDER [ integer, nlink, (N_det_ref) ] +&BEGIN_PROVIDER [ integer, linked, (N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ integer, blokMwen, (N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, searchance, (N_det_ref) ] +&BEGIN_PROVIDER [ integer, child_num, (N_det_non_ref,N_det_ref) ] + + use bitmasks + implicit none + + integer(bit_kind),allocatable :: det_noactive(:,:,:) + integer, allocatable :: shortcut(:), idx(:) + integer(bit_kind) :: nonactive_sorb(N_int,2), det(N_int, 2) + integer i, II, j, k, n, ni, blok, degree + logical, external :: detEq + + allocate(det_noactive(N_int, 2, N_det_non_ref)) + allocate(idx(N_det_non_ref), shortcut(0:N_det_non_ref+1)) + print *, "pre start" + active_sorb(:,:) = 0_8 + nonactive_sorb(:,:) = not(0_8) + + if(N_det_ref > 1) then + do i=1, N_det_ref + do k=1, N_int + active_sorb(k,1) = ior(psi_ref(k,1,i), active_sorb(k,1)) + active_sorb(k,2) = ior(psi_ref(k,2,i), active_sorb(k,2)) + nonactive_sorb(k,1) = iand(psi_ref(k,1,i), nonactive_sorb(k,1)) + nonactive_sorb(k,2) = iand(psi_ref(k,2,i), nonactive_sorb(k,2)) + end do + end do + do k=1, N_int + active_sorb(k,1) = iand(active_sorb(k,1), not(nonactive_sorb(k,1))) + active_sorb(k,2) = iand(active_sorb(k,2), not(nonactive_sorb(k,2))) + end do + end if + + + do i=1, N_det_non_ref + do k=1, N_int + det_noactive(k,1,i) = iand(psi_non_ref(k,1,i), not(active_sorb(k,1))) + det_noactive(k,2,i) = iand(psi_non_ref(k,2,i), not(active_sorb(k,2))) + end do + end do + + call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int) + + do i=1,N_det_non_ref + det_cepa0(:,:,i) = psi_non_ref(:,:,det_cepa0_idx(i)) + end do + + cepa0_shortcut(0) = 1 + cepa0_shortcut(1) = 1 + do i=2,N_det_non_ref + if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then + cepa0_shortcut(0) += 1 + cepa0_shortcut(cepa0_shortcut(0)) = i + end if + end do + cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1 + + if(.true.) then + do i=1,cepa0_shortcut(0) + n = cepa0_shortcut(i+1) - cepa0_shortcut(i) + call sort_dets_ab(det_cepa0(1,1,cepa0_shortcut(i)), idx, shortcut, n, N_int) + do k=1,n + idx(k) = det_cepa0_idx(cepa0_shortcut(i)-1+idx(k)) + end do + det_cepa0_idx(cepa0_shortcut(i):cepa0_shortcut(i)+n-1) = idx(:n) + end do + end if + + + do i=1,N_det_ref + do k=1, N_int + det_ref_active(k,1,i) = iand(psi_ref(k,1,i), active_sorb(k,1)) + det_ref_active(k,2,i) = iand(psi_ref(k,2,i), active_sorb(k,2)) + end do + end do + + do i=1,N_det_non_ref + do k=1, N_int + det_cepa0_active(k,1,i) = iand(psi_non_ref(k,1,det_cepa0_idx(i)), active_sorb(k,1)) + det_cepa0_active(k,2,i) = iand(psi_non_ref(k,2,det_cepa0_idx(i)), active_sorb(k,2)) + end do + end do + + do i=1,N_det_non_ref + if(.not. detEq(psi_non_ref(1,1,det_cepa0_idx(i)), det_cepa0(1,1,i),N_int)) stop "STOOOP" + end do + + searchance = 0d0 + child_num = 0 + do J = 1, N_det_ref + nlink(J) = 0 + do blok=1,cepa0_shortcut(0) + do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int) + if(degree <= 2) then + nlink(J) += 1 + linked(nlink(J),J) = k + child_num(k, J) = nlink(J) + blokMwen(nlink(J),J) = blok + searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok))) + end if + end do + end do + end do + print *, "pre done" +END_PROVIDER + + +! BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] +! use bitmasks +! implicit none +! integer :: i,j,k +! double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall +! integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref) +! +! ! provide lambda_mrcc +! npres = 0 +! delta_cas = 0d0 +! call wall_time(wall) +! print *, "dcas ", wall +! do i_state = 1, N_states +! !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref) +! do k=1,N_det_non_ref +! if(lambda_mrcc(i_state, k) == 0d0) cycle +! npre = 0 +! do i=1,N_det_ref +! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) +! if(Hki /= 0d0) then +! !!$OMP ATOMIC +! npres(i) += 1 +! npre += 1 +! ipre(npre) = i +! pre(npre) = Hki +! end if +! end do +! +! +! do i=1,npre +! do j=1,i +! !!$OMP ATOMIC +! delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k) +! end do +! end do +! end do +! !!$OMP END PARALLEL DO +! npre=0 +! do i=1,N_det_ref +! npre += npres(i) +! end do +! !stop +! do i=1,N_det_ref +! do j=1,i +! delta_cas(j,i,i_state) = delta_cas(i,j,i_state) +! end do +! end do +! end do +! +! call wall_time(wall) +! print *, "dcas", wall +! ! stop +! END_PROVIDER + + + BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] + use bitmasks + implicit none + integer :: i,j,k + double precision :: Hjk, Hki, Hij + !double precision, external :: get_dij + integer i_state, degree + + provide lambda_mrcc dIj + do i_state = 1, N_states + !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref,dij) + do i=1,N_det_ref + do j=1,i + call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int) + delta_cas(i,j,i_state) = 0d0 + do k=1,N_det_non_ref + + call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) + + delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) + !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int) + end do + delta_cas(j,i,i_state) = delta_cas(i,j,i_state) + end do + end do + !$OMP END PARALLEL DO + end do + END_PROVIDER + + + + +logical function isInCassd(a,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2) + integer(bit_kind) :: inac, virt + integer :: ni, i, deg + + + isInCassd = .false. + + deg = 0 + do i=1,2 + do ni=1,Nint + virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i))) + deg += popcnt(iand(virt, a(ni,i))) + if(deg > 2) return + end do + end do + + deg = 0 + do i=1,2 + do ni=1,Nint + inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i))) + deg += popcnt(xor(iand(inac,a(ni,i)), inac)) + if(deg > 2) return + end do + end do + isInCassd = .true. +end function + + +subroutine getHP(a,h,p,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2) + integer, intent(out) :: h, p + integer(bit_kind) :: inac, virt + integer :: ni, i, deg + + + !isInCassd = .false. + h = 0 + p = 0 + + deg = 0 + lp : do i=1,2 + do ni=1,Nint + virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i))) + deg += popcnt(iand(virt, a(ni,i))) + if(deg > 2) exit lp + end do + end do lp + p = deg + + deg = 0 + lh : do i=1,2 + do ni=1,Nint + inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i))) + deg += popcnt(xor(iand(inac,a(ni,i)), inac)) + if(deg > 2) exit lh + end do + end do lh + h = deg + !isInCassd = .true. +end function + + + BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref) + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) + double precision :: contrib, contrib2, HIIi, HJk, wall + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2) + integer(bit_kind),allocatable :: sortRef(:,:,:) + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit, searchDet + logical, external :: is_in_wavefunction, detEq + !double precision, external :: get_dij + integer :: II, blok + integer*8, save :: notf = 0 + + call wall_time(wall) + allocate(idx_sorted_bit(N_det), sortRef(N_int,2,N_det_ref)) + + sortRef(:,:,:) = det_ref_active(:,:,:) + call sort_det(sortRef, sortRefIdx, N_det_ref, N_int) + + idx_sorted_bit(:) = -1 + do i=1,N_det_non_ref + idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i + enddo + + ! To provide everything + contrib = dij(1, 1, 1) + + do i_state = 1, N_states + delta_mrcepa0_ii(:,:) = 0d0 + delta_mrcepa0_ij(:,:,:) = 0d0 + + !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii) & + !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2) & + !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & + !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) & + !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) + do blok=1,cepa0_shortcut(0) + do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do II=1,N_det_ref + call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int) + if (degree > 2 ) cycle + + do ni=1,N_int + made_hole(ni,1) = iand(det_ref_active(ni,1,II), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II))) + made_hole(ni,2) = iand(det_ref_active(ni,2,II), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II))) + + made_particle(ni,1) = iand(det_cepa0_active(ni,1,i), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II))) + made_particle(ni,2) = iand(det_cepa0_active(ni,2,i), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II))) + end do + + + kloop: do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 !i + !if(lambda_mrcc(i_state, det_cepa0_idx(k)) == 0d0) cycle + + do ni=1,N_int + if(iand(made_hole(ni,1), det_cepa0_active(ni,1,k)) /= 0) cycle kloop + if(iand(made_particle(ni,1), det_cepa0_active(ni,1,k)) /= made_particle(ni,1)) cycle kloop + if(iand(made_hole(ni,2), det_cepa0_active(ni,2,k)) /= 0) cycle kloop + if(iand(made_particle(ni,2), det_cepa0_active(ni,2,k)) /= made_particle(ni,2)) cycle kloop + end do + do ni=1,N_int + myActive(ni,1) = xor(det_cepa0_active(ni,1,k), made_hole(ni,1)) + myActive(ni,1) = xor(myActive(ni,1), made_particle(ni,1)) + myActive(ni,2) = xor(det_cepa0_active(ni,2,k), made_hole(ni,2)) + myActive(ni,2) = xor(myActive(ni,2), made_particle(ni,2)) + end do + + j = searchDet(sortRef, myActive, N_det_ref, N_int) + if(j == -1) then + cycle + end if + j = sortRefIdx(j) + !$OMP ATOMIC + notf = notf+1 + +! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) + contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) + !$OMP ATOMIC + delta_mrcepa0_ii(J,i_state) -= contrib2 + else + contrib = contrib * 0.5d0 + end if + !$OMP ATOMIC + delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib + + end do kloop + end do + end do + end do + !$OMP END PARALLEL DO + end do + deallocate(idx_sorted_bit) + call wall_time(wall) + print *, "cepa0", wall, notf + !stop +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ + logical :: ok + double precision :: phase_Ji, phase_Ik, phase_Ii + double precision :: contrib, contrib2, delta_IJk, HJk, HIk, HIl + integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2) + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + + integer :: II, blok + + provide delta_cas lambda_mrcc + allocate(idx_sorted_bit(N_det)) + idx_sorted_bit(:) = -1 + do i=1,N_det_non_ref + idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i + enddo + + do i_state = 1, N_states + delta_sub_ij(:,:,:) = 0d0 + delta_sub_ii(:,:) = 0d0 + + provide mo_bielec_integrals_in_map + + + !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) & + !$OMP private(i, J, k, degree, degree2, l, deg, ni) & + !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) & + !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & + !$OMP private(det_tmp, det_tmp2, II, blok) & + !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) & + !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb) + do i=1,N_det_non_ref + if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref + do J=1,N_det_ref + call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int) + if(degree == -1) cycle + + + do II=1,N_det_ref + call apply_excitation(psi_ref(1,1,II),exc_Ji,det_tmp,ok,N_int) + + if(.not. ok) cycle + l = get_index_in_psi_det_sorted_bit(det_tmp, N_int) + if(l == 0) cycle + l = idx_sorted_bit(l) + + call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl) + + do k=1,N_det_non_ref + if(lambda_mrcc(i_state, k) == 0d0) cycle + call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int) + + det_tmp(:,:) = 0_bit_kind + det_tmp2(:,:) = 0_bit_kind + + ok = .true. + do ni=1,N_int + det_tmp(ni,1) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,k)), not(active_sorb(ni,1))) + det_tmp(ni,2) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,i)), not(active_sorb(ni,1))) + ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2)))) + + det_tmp(ni,1) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,k)), not(active_sorb(ni,2))) + det_tmp(ni,2) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,i)), not(active_sorb(ni,2))) + ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2)))) + end do + + if(ok) cycle + + + call i_h_j(psi_ref(1,1,J), psi_non_ref(1,1,k), N_int, HJk) + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,k), N_int, HIk) + if(HJk == 0) cycle + !assert HIk == 0 + delta_IJk = HJk * HIk * lambda_mrcc(i_state, k) + call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) + if(ok) cycle + contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) + if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then + contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) + !$OMP ATOMIC + delta_sub_ii(II,i_state) -= contrib2 + else + contrib = contrib * 0.5d0 + endif + !$OMP ATOMIC + delta_sub_ij(II, i, i_state) += contrib + end do + end do + end do + end do + !$OMP END PARALLEL DO + end do + deallocate(idx_sorted_bit) +END_PROVIDER + + +subroutine set_det_bit(det, p, s) + implicit none + integer(bit_kind),intent(inout) :: det(N_int, 2) + integer, intent(in) :: p, s + integer :: ni, pos + + ni = (p-1)/bit_kind_size + 1 + pos = mod(p-1, bit_kind_size) + det(ni,s) = ibset(det(ni,s), pos) +end subroutine + + +BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ] + implicit none + integer :: i,j + do i=1,N_det_ref + do j=1,N_det_non_ref + call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,j), N_int, h_(i,j)) + end do + end do +END_PROVIDER + + + +subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) + + 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,allocatable :: degree(:) + integer,allocatable :: idx(:) + logical :: good + + integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators) + integer,intent(in) :: N_miniList + + allocate(degree(psi_det_size)) + allocate(idx(0:psi_det_size)) + N_tq = 0 + + i_loop : do i=1,N_selected + do k=1, N_minilist + if(is_generable(miniList(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + + ! 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)) 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 + + +subroutine filter_tq_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,allocatable :: degree(:) + integer,allocatable :: idx(:) + logical :: good + + integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + 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 + + + allocate(degree(psi_det_size)) + allocate(idx(0:psi_det_size)) + 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 + do k=ptr_microlist(smallerlist), ptr_microlist(smallerlist)+N_microlist(smallerlist)-1 + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + end if + + if(N_microlist(0) > 0) then + do k=1, N_microlist(0) + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + 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)) 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_selected/dressing_slave.irp.f b/plugins/mrcc_selected/dressing_slave.irp.f new file mode 100644 index 00000000..f1d6f029 --- /dev/null +++ b/plugins/mrcc_selected/dressing_slave.irp.f @@ -0,0 +1,593 @@ +subroutine mrsc2_dressing_slave_tcp(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + call mrsc2_dressing_slave(0,i) +end + + +subroutine mrsc2_dressing_slave_inproc(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + call mrsc2_dressing_slave(1,i) +end + +subroutine mrsc2_dressing_slave(thread,iproc) + use f77_zmq + + implicit none + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + integer, intent(in) :: thread, iproc +! integer :: j,l + integer :: rc + + integer :: worker_id, task_id + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + double precision, allocatable :: delta(:,:,:) + + + + integer :: i_state, i, i_I, J, k, k2, k1, kk, ll, degree, degree2, m, l, deg, ni, m2 + integer :: n(2) + integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, kn + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al + double precision :: diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv(N_states), cj_inv(N_states) + double precision :: contrib, wall, iwall + double precision, allocatable :: dleat(:,:,:) + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt + integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp + logical, external :: is_in_wavefunction, isInCassd, detEq + integer,allocatable :: komon(:) + logical :: komoned + !double precision, external :: get_dij + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_push = new_zmq_push_socket(thread) + + call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) + + allocate (dleat(N_states, N_det_non_ref, 2), delta(N_states,0:N_det_non_ref, 2)) + allocate(komon(0:N_det_non_ref)) + + do + call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) + if (task_id == 0) exit + read (task,*) i_I, J, k1, k2 + do i_state=1, N_states + ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state) + cj_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state) + end do + !delta = 0.d0 + n = 0 + delta(:,0,:) = 0d0 + delta(:,:nlink(J),1) = 0d0 + delta(:,:nlink(i_I),2) = 0d0 + komon(0) = 0 + komoned = .false. + + + + + do kk = k1, k2 + k = det_cepa0_idx(linked(kk, i_I)) + blok = blokMwen(kk, i_I) + + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int) + + if(J /= i_I) then + call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int) + if(.not. ok) cycle + + l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int) + if(l == -1) cycle + ll = cepa0_shortcut(blok)-1+l + l = det_cepa0_idx(ll) + ll = child_num(ll, J) + else + l = k + ll = kk + end if + + + if(.not. komoned) then + m = 0 + m2 = 0 + + do while(m < nlink(i_I) .and. m2 < nlink(J)) + m += 1 + m2 += 1 + if(linked(m, i_I) < linked(m2, J)) then + m2 -= 1 + cycle + else if(linked(m, i_I) > linked(m2, J)) then + m -= 1 + cycle + end if + i = det_cepa0_idx(linked(m, i_I)) + + if(h_(J,i) == 0.d0) cycle + if(h_(i_I,i) == 0.d0) cycle + + !ok = .false. + !do i_state=1, N_states + ! if(lambda_mrcc(i_state, i) /= 0d0) then + ! ok = .true. + ! exit + ! end if + !end do + !if(.not. ok) cycle +! + + komon(0) += 1 + kn = komon(0) + komon(kn) = i + + +! call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int) +! if(I_i /= J) call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ii,N_int) +! if(I_i == J) phase_Ii = phase_Ji + + do i_state = 1,N_states + dkI = h_(J,i) * dij(i_I, i, i_state)!get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int) + !dkI = h_(J,i) * h_(i_I,i) * lambda_mrcc(i_state, i) + dleat(i_state, kn, 1) = dkI + dleat(i_state, kn, 2) = dkI + end do + + end do + + komoned = .true. + end if + + + do m = 1, komon(0) + + i = komon(m) + + call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) + if(.not. ok) cycle + if(HP(1,i) + HP(1,k) <= 2 .and. HP(2,i) + HP(2,k) <= 2) then +! if(is_in_wavefunction(det_tmp, N_int)) cycle + cycle + end if + + !if(isInCassd(det_tmp, N_int)) cycle + + do i_state = 1, N_states + !if(lambda_mrcc(i_state, i) == 0d0) cycle + + + !contrib = h_(i_I,k) * lambda_mrcc(i_state, k) * dleat(i_state, m, 2)! * phase_al + contrib = dij(i_I, k, i_state) * dleat(i_state, m, 2) + delta(i_state,ll,1) += contrib + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then + delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state) + endif + + if(I_i == J) cycle + !contrib = h_(J,l) * lambda_mrcc(i_state, l) * dleat(i_state, m, 1)! * phase_al + contrib = dij(J, l, i_state) * dleat(i_state, m, 1) + delta(i_state,kk,2) += contrib + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state) + end if + enddo !i_state + end do ! while + end do ! kk + + + call push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) + +! end if + + enddo + + deallocate(delta) + + call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + +end + + +subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) + use f77_zmq + implicit none + BEGIN_DOC +! Push integrals in the push socket + END_DOC + + integer, intent(in) :: i_I, J + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + double precision,intent(inout) :: delta(N_states, 0:N_det_non_ref, 2) + integer, intent(in) :: task_id + integer :: rc , i_state, i, kk, li + integer,allocatable :: idx(:,:) + integer :: n(2) + logical :: ok + + allocate(idx(N_det_non_ref,2)) + rc = f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + + do kk=1,2 + n(kk)=0 + if(kk == 1) li = nlink(j) + if(kk == 2) li = nlink(i_I) + do i=1, li + ok = .false. + do i_state=1,N_states + if(delta(i_state, i, kk) /= 0d0) then + ok = .true. + exit + end if + end do + + if(ok) then + n(kk) += 1 +! idx(n,kk) = i + if(kk == 1) then + idx(n(1),1) = det_cepa0_idx(linked(i, J)) + else + idx(n(2),2) = det_cepa0_idx(linked(i, i_I)) + end if + + do i_state=1, N_states + delta(i_state, n(kk), kk) = delta(i_state, i, kk) + end do + end if + end do + + rc = f77_zmq_send( zmq_socket_push, n(kk), 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, n, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + if(n(kk) /= 0) then + rc = f77_zmq_send( zmq_socket_push, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) ! delta(1,0,1) = delta_I delta(1,0,2) = delta_J + if (rc /= (n(kk)+1)*8*N_states) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) + if (rc /= n(kk)*4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, delta, 8*n(kk), ZMQ_SNDMORE)' + stop 'error' + endif + end if + end do + + + rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, task_id, 4, 0)' + stop 'error' + endif + +! ! Activate is zmq_socket_push is a REQ +! integer :: idummy +! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) +! if (rc /= 4) then +! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)' +! stop 'error' +! endif +end + + + +subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) + use f77_zmq + implicit none + BEGIN_DOC +! Push integrals in the push socket + END_DOC + + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + integer, intent(out) :: i_I, J, n(2) + double precision, intent(inout) :: delta(N_states, 0:N_det_non_ref, 2) + integer, intent(out) :: task_id + integer :: rc , i, kk + integer,intent(inout) :: idx(N_det_non_ref,2) + logical :: ok + + rc = f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + do kk = 1, 2 + rc = f77_zmq_recv( zmq_socket_pull, n(kk), 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, n, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + if(n(kk) /= 0) then + rc = f77_zmq_recv( zmq_socket_pull, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) + if (rc /= (n(kk)+1)*8*N_states) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) + if (rc /= n(kk)*4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, delta, n(kk)*4, ZMQ_SNDMORE)' + stop 'error' + endif + end if + end do + + rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)' + stop 'error' + endif + + +! ! Activate is zmq_socket_pull is a REP +! integer :: idummy +! rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0) +! if (rc /= 4) then +! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)' +! stop 'error' +! endif +end + + + +subroutine mrsc2_dressing_collector(delta_ii_,delta_ij_) + use f77_zmq + implicit none + BEGIN_DOC +! Collects results from the AO integral calculation + END_DOC + + double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) + double precision,intent(inout) :: delta_ii_(N_states,N_det_ref) + +! integer :: j,l + integer :: rc + + double precision, allocatable :: delta(:,:,:) + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer*8 :: control, accu + integer :: task_id, more + + integer :: I_i, J, l, i_state, n(2), kk + integer,allocatable :: idx(:,:) + + delta_ii_(:,:) = 0d0 + delta_ij_(:,:,:) = 0d0 + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + + allocate ( delta(N_states,0:N_det_non_ref,2) ) + + allocate(idx(N_det_non_ref,2)) + more = 1 + do while (more == 1) + + call pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) + + + do l=1, n(1) + do i_state=1,N_states + delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1) + end do + end do + + do l=1, n(2) + do i_state=1,N_states + delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2) + end do + end do + + +! +! do l=1,nlink(J) +! do i_state=1,N_states +! delta_ij_(i_state,det_cepa0_idx(linked(l,J)),i_I) += delta(i_state,l,1) +! delta_ij_(i_state,det_cepa0_idx(linked(l,i_I)),j) += delta(i_state,l,2) +! end do +! end do +! + if(n(1) /= 0) then + do i_state=1,N_states + delta_ii_(i_state,i_I) += delta(i_state,0,1) + end do + end if + + if(n(2) /= 0) then + do i_state=1,N_states + delta_ii_(i_state,J) += delta(i_state,0,2) + end do + end if + + + if (task_id /= 0) then + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) + endif + + + enddo + deallocate( delta ) + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + +end + + + + + BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ] + implicit none + + integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2 + integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, nex, nzer, ntot +! integer, allocatable :: linked(:,:), blokMwen(:, :), nlink(:) + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states) + double precision :: contrib, wall, iwall ! , searchance(N_det_ref) + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt + integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp + logical, external :: is_in_wavefunction, isInCassd, detEq + character*(512) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer :: KKsize = 1000000 + + + call new_parallel_job(zmq_to_qp_run_socket,'mrsc2') + + + call wall_time(iwall) +! allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref)) + + +! searchance = 0d0 +! do J = 1, N_det_ref +! nlink(J) = 0 +! do blok=1,cepa0_shortcut(0) +! do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 +! call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int) +! if(degree <= 2) then +! nlink(J) += 1 +! linked(nlink(J),J) = k +! blokMwen(nlink(J),J) = blok +! searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok))) +! end if +! end do +! end do +! end do + + + +! stop + nzer = 0 + ntot = 0 + do nex = 3, 0, -1 + print *, "los ",nex + do I_s = N_det_ref, 1, -1 +! if(mod(I_s,1) == 0) then +! call wall_time(wall) +! wall = wall-iwall +! print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall +! end if + + + do J_s = 1, I_s + + call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int) + if(degree /= nex) cycle + if(nex == 3) nzer = nzer + 1 + ntot += 1 +! if(degree > 3) then +! deg += 1 +! cycle +! else if(degree == -10) then +! KKsize = 100000 +! else +! KKsize = 1000000 +! end if + + + + if(searchance(I_s) < searchance(J_s)) then + i_I = I_s + J = J_s + else + i_I = J_s + J = I_s + end if + + KKsize = nlink(1) + if(nex == 0) KKsize = int(float(nlink(1)) / float(nlink(i_I)) * (float(nlink(1)) / 64d0)) + + !if(KKsize == 0) stop "ZZEO" + + do kk = 1 , nlink(i_I), KKsize + write(task,*) I_i, J, kk, int(min(kk+KKsize-1, nlink(i_I))) + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + end do + + ! do kk = 1 , nlink(i_I) + ! k = linked(kk,i_I) + ! blok = blokMwen(kk,i_I) + ! write(task,*) I_i, J, k, blok + ! call add_task_to_taskserver(zmq_to_qp_run_socket,task) + ! + ! enddo !kk + enddo !J + + enddo !I + end do ! nex + print *, "tasked" +! integer(ZMQ_PTR) ∷ collector_thread +! external ∷ ao_bielec_integrals_in_map_collector +! rc = pthread_create(collector_thread, mrsc2_dressing_collector) + print *, nzer, ntot, float(nzer) / float(ntot) + provide nproc + !$OMP PARALLEL DEFAULT(none) SHARED(delta_ii_old,delta_ij_old) PRIVATE(i) NUM_THREADS(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call mrsc2_dressing_collector(delta_ii_old,delta_ij_old) + else + call mrsc2_dressing_slave_inproc(i) + endif + !$OMP END PARALLEL + +! rc = pthread_join(collector_thread) + call end_parallel_job(zmq_to_qp_run_socket, 'mrsc2') + + +END_PROVIDER + + + diff --git a/plugins/mrcc_selected/ezfio_interface.irp.f b/plugins/mrcc_selected/ezfio_interface.irp.f new file mode 100644 index 00000000..062af449 --- /dev/null +++ b/plugins/mrcc_selected/ezfio_interface.irp.f @@ -0,0 +1,61 @@ +! DO NOT MODIFY BY HAND +! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py +! from file /home/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg + + +BEGIN_PROVIDER [ double precision, thresh_dressed_ci ] + implicit none + BEGIN_DOC +! Threshold on the convergence of the dressed CI energy + END_DOC + + logical :: has + PROVIDE ezfio_filename + + call ezfio_has_mrcc_selected_thresh_dressed_ci(has) + if (has) then + call ezfio_get_mrcc_selected_thresh_dressed_ci(thresh_dressed_ci) + else + print *, 'mrcc_selected/thresh_dressed_ci not found in EZFIO file' + stop 1 + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_it_max_dressed_ci ] + implicit none + BEGIN_DOC +! Maximum number of dressed CI iterations + END_DOC + + logical :: has + PROVIDE ezfio_filename + + call ezfio_has_mrcc_selected_n_it_max_dressed_ci(has) + if (has) then + call ezfio_get_mrcc_selected_n_it_max_dressed_ci(n_it_max_dressed_ci) + else + print *, 'mrcc_selected/n_it_max_dressed_ci not found in EZFIO file' + stop 1 + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, lambda_type ] + implicit none + BEGIN_DOC +! lambda type + END_DOC + + logical :: has + PROVIDE ezfio_filename + + call ezfio_has_mrcc_selected_lambda_type(has) + if (has) then + call ezfio_get_mrcc_selected_lambda_type(lambda_type) + else + print *, 'mrcc_selected/lambda_type not found in EZFIO file' + stop 1 + endif + +END_PROVIDER diff --git a/plugins/mrcc_selected/mrcc_selected.irp.f b/plugins/mrcc_selected/mrcc_selected.irp.f new file mode 100644 index 00000000..91592e62 --- /dev/null +++ b/plugins/mrcc_selected/mrcc_selected.irp.f @@ -0,0 +1,19 @@ +program mrsc2sub + implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + mrmode = 3 + + read_wf = .True. + SOFT_TOUCH read_wf + 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 + diff --git a/plugins/mrcc_selected/mrcepa0_general.irp.f b/plugins/mrcc_selected/mrcepa0_general.irp.f new file mode 100644 index 00000000..e3a2d1f5 --- /dev/null +++ b/plugins/mrcc_selected/mrcepa0_general.irp.f @@ -0,0 +1,245 @@ + + +subroutine run(N_st,energy) + implicit none + + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration + double precision :: E_past(4) + + integer :: n_it_mrcc_max + double precision :: thresh_mrcc + double precision, allocatable :: lambda(:) + allocate (lambda(N_states)) + + + thresh_mrcc = thresh_dressed_ci + n_it_mrcc_max = n_it_max_dressed_ci + + if(n_it_mrcc_max == 1) then + 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_mrcepa0_energy(ci_energy_dressed(1)) + call save_wavefunction + energy(:) = ci_energy_dressed(:) + else + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + lambda = 1.d0 + do while (delta_E > thresh_mrcc) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA0 Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA0 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_mrcepa0_energy(ci_energy_dressed(1)) + if (iteration >= n_it_mrcc_max) then + exit + endif + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") + energy(:) = ci_energy_dressed(:) + endif +end + + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, (psi_cas_coef(i,j), j=1,N_states) + call debug_det(psi_cas(1,1,i),N_int) + enddo + call write_double(6,ci_energy(1),"Initial CI energy") + +end + + + + +subroutine run_pt2_old(N_st,energy) + implicit none + integer :: i,j,k + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + double precision :: pt2_redundant(N_st), pt2(N_st) + double precision :: norm_pert(N_st),H_pert_diag(N_st) + + pt2_redundant = 0.d0 + pt2 = 0d0 + !if(lambda_mrcc_pt2(0) == 0) return + + print*,'Last iteration only to compute the PT2' + + print * ,'Computing the redundant PT2 contribution' + + if (mrmode == 1) then + + N_det_generators = lambda_mrcc_kept(0) + N_det_selectors = lambda_mrcc_kept(0) + + do i=1,N_det_generators + j = lambda_mrcc_kept(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) + psi_selectors(k,1,i) = psi_non_ref(k,1,j) + psi_selectors(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) + psi_selectors_coef(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + + else + + N_det_generators = N_det_non_ref + N_det_selectors = N_det_non_ref + + do i=1,N_det_generators + j = 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) + psi_selectors(k,1,i) = psi_non_ref(k,1,j) + psi_selectors(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) + psi_selectors_coef(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + + endif + + SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized + + call H_apply_mrcepa_PT2(pt2_redundant, norm_pert, H_pert_diag, N_st) + + print * ,'Computing the remaining contribution' + + threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + threshold_generators = max(threshold_generators,threshold_generators_pt2) + + N_det_generators = N_det_non_ref + N_det_ref + N_det_selectors = N_det_non_ref + N_det_ref + + psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + + do i=N_det_ref+1,N_det_generators + j = i-N_det_ref + 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) + psi_selectors(k,1,i) = psi_non_ref(k,1,j) + psi_selectors(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) + psi_selectors_coef(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + + SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized + + call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st) + + + print *, "Redundant PT2 :",pt2_redundant + print *, "Full PT2 :",pt2 + print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1) + pt2 = pt2 - pt2_redundant + + 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_mrcepa0_energy_pt2(energy(1)+pt2(1)) + +end + +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + double precision :: pt2(N_st) + double precision :: norm_pert(N_st),H_pert_diag(N_st) + + pt2 = 0d0 + !if(lambda_mrcc_pt2(0) == 0) return + + print*,'Last iteration only to compute the PT2' + + N_det_generators = N_det_cas + N_det_selectors = N_det_non_ref + + do i=1,N_det_generators + do k=1,N_int + psi_det_generators(k,1,i) = psi_ref(k,1,i) + psi_det_generators(k,2,i) = psi_ref(k,2,i) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_ref_coef(i,k) + enddo + enddo + do i=1,N_det + do k=1,N_int + psi_selectors(k,1,i) = psi_det_sorted(k,1,i) + psi_selectors(k,2,i) = psi_det_sorted(k,2,i) + enddo + do k=1,N_st + psi_selectors_coef(i,k) = psi_coef_sorted(i,k) + enddo + enddo + + SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized + + call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st) + +! call ezfio_set_full_ci_energy_pt2(energy+pt2) + + 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_mrcepa0_energy_pt2(energy(1)+pt2(1)) + +end + diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 39b0f58e..f690c790 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -438,8 +438,12 @@ end do i=1,N_states psi_coef_min(i) = minval(psi_coef(:,i)) psi_coef_max(i) = maxval(psi_coef(:,i)) - abs_psi_coef_min(i) = dabs(psi_coef_min(i)) - abs_psi_coef_max(i) = dabs(psi_coef_max(i)) + abs_psi_coef_min(i) = minval( dabs(psi_coef(:,i)) ) + abs_psi_coef_max(i) = maxval( dabs(psi_coef(:,i)) ) + call write_double(6,psi_coef_max(i), 'Max coef') + call write_double(6,psi_coef_min(i), 'Min coef') + call write_double(6,abs_psi_coef_max(i), 'Max abs coef') + call write_double(6,abs_psi_coef_min(i), 'Min abs coef') enddo END_PROVIDER