diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 7314674c..1cafc8de 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -24,5 +24,27 @@ s.data["size_max"] = "3072" print s +s = H_apply("mrcepa") +s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" +s.data["declarations"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["decls_main"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["finalization"] = "" +s.data["copy_buffer"] = "" +s.data["generate_psi_guess"] = "" +s.data["size_max"] = "3072" +# print s + + END_SHELL diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 8fe6c411..1e2f974d 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -26,7 +26,11 @@ ! TODO --- Test perturbatif ------ do k=1,N_states lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + ! TODO : i_h_psi peut sortir de la boucle? 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) + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif tmp = psi_non_ref_coef(i,k)/ihpsi_current(k) i_pert = 0 ! Perturbation only if 1st order < 0.5 x second order diff --git a/plugins/MRCC_Utils/mrcepa_dress.irp.f b/plugins/MRCC_Utils/mrcepa_dress.irp.f new file mode 100644 index 00000000..9789b818 --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_dress.irp.f @@ -0,0 +1,260 @@ +use omp_lib +use bitmasks + +subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_ref, Ndet_non_ref + double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(inout) :: delta_ii_(Ndet_ref,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l + integer :: degree_alpha(psi_det_size) + integer :: idx_alpha(0:psi_det_size) + logical :: good, fullMatch + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref ,degree + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hia(:,:) + 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(bit_kind) :: tmp_det_0(Nint,2) + integer :: iint, ipos + integer :: i_state, i_sd, 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 + integer(bit_kind) :: isum + + double precision :: hia + integer, allocatable :: index_sorted(:) + + + leng = max(N_det_generators, N_det_non_ref) + allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det)) + + !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 + + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + + allocate (dIa_hia(N_states,Ndet_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) + 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) + + integer, external :: get_index_in_psi_det_sorted_bit + index_sorted = huge(-1) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j) + end do + + ! |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 + + !TODO: MR + do i_sd=1,idx_alpha(0) + call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint) + if (degree > 2) then + cycle + endif + call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det_0 = 0_bit_kind + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + if (degree == 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_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + endif + + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia) + + ! |alpha> + do k_sd=1,idx_alpha(0) + 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 + + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det = 0_bit_kind + ! 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) = ibset(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 == 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) = ibset(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 + + isum = 0_bit_kind + do iint = 1,N_int + isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) & + + iand(tmp_det(iint,2), tmp_det_0(iint,2)) + enddo + + if (isum /= 0_bit_kind) then + cycle + endif + + ! + ! + 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) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + 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 + ! 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 == 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 + + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + + +! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) ) +! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl) +! do i_state=1,N_states +! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2 +! enddo + + do l_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 + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + 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) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + enddo + 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) = 1.d0/psi_ref_coef(i_I,i_state) + enddo + + k_sd = idx_alpha(i_sd) + do i_state=1,N_states + dIa_hia(i_state,k_sd) = dIa(i_state) * hia + enddo + + call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states + delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd) + + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) + else + delta_ii_(i_I,i_state) = 0.d0 + endif + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) + enddo + enddo + + enddo + deallocate (dIa_hia,index_sorted) + deallocate(miniList, idx_miniList) +end + + + + + + + + diff --git a/plugins/MRCC_Utils/mrcepa_general.irp.f b/plugins/MRCC_Utils/mrcepa_general.irp.f new file mode 100644 index 00000000..3479548b --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_general.irp.f @@ -0,0 +1,97 @@ +subroutine run_mrcepa + implicit none + call set_generators_bitmasks_as_holes_and_particles + call mrcepa_iterations +end + +subroutine mrcepa_iterations + implicit none + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration,i_oscillations + double precision :: E_past(4) + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + j = 1 + i_oscillations = 0 + do while (delta_E > 1.d-7) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + + E_past(j) = E_new + j +=1 + if(j>4)then + j=1 + endif + if(iteration > 4) then + if(delta_E > 1.d-10)then + if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then + print*,'OSCILLATIONS !!!' + oscillations = .True. + i_oscillations +=1 + lambda_mrcc_tmp = lambda_mrcc + endif + endif + endif + call save_wavefunction +! if (i_oscillations > 5) then +! exit +! endif + if (iteration > 200) then + exit + endif + print*,'------------' + print*,'VECTOR' + do i = 1, N_det_ref + print*,'' + print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) + print*,'delta_ii(i,1) = ',delta_ii(i,1) + enddo + print*,'------------' + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + +end + +subroutine set_generators_bitmasks_as_holes_and_particles + implicit none + integer :: i,k + do k = 1, N_generators_bitmask + do i = 1, N_int + ! Pure single part + generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha + generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha + generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta + generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta + + ! Double excitation + generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha + generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha + generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta + generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta + + generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha + generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha + generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta + generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta + + enddo + enddo + touch generators_bitmask + + + +end