From 7f583946c27ccb9b4689790ce29b2f45f56d4553 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Fri, 11 Mar 2016 19:35:57 +0100 Subject: [PATCH] experimental mrcepa0/mrsc2 --- ocaml/.gitignore | 1 - ocaml/qp_edit.ml | 66 ++-- plugins/MRCC_Utils/davidson.irp.f | 1 - plugins/MRCC_Utils/mrcc_utils.irp.f | 21 +- plugins/mrcepa0/dressing.irp.f | 439 +++++++++++++++++++++++--- plugins/mrcepa0/mrcepa0_general.irp.f | 2 +- 6 files changed, 432 insertions(+), 98 deletions(-) diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 45d71ee3..0f0c1ef9 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -4,7 +4,6 @@ ezfio.ml Git.ml Input_auto_generated.ml Input_determinants.ml -Input_foboci.ml Input_hartree_fock.ml Input_integrals_bielec.ml Input_perturbation.ml diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml index 67dc9501..409387b2 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -17,13 +17,12 @@ type keyword = | Electrons | Mo_basis | Nuclei -| Hartree_fock -| Pseudo -| Integrals_bielec -| Perturbation -| Properties -| Foboci | Determinants +| Integrals_bielec +| Pseudo +| Perturbation +| Hartree_fock +| Properties ;; @@ -33,13 +32,12 @@ let keyword_to_string = function | Electrons -> "Electrons" | Mo_basis -> "MO basis" | Nuclei -> "Molecule" -| Hartree_fock -> "Hartree_fock" -| Pseudo -> "Pseudo" -| Integrals_bielec -> "Integrals_bielec" -| Perturbation -> "Perturbation" -| Properties -> "Properties" -| Foboci -> "Foboci" | Determinants -> "Determinants" +| Integrals_bielec -> "Integrals_bielec" +| Pseudo -> "Pseudo" +| Perturbation -> "Perturbation" +| Hartree_fock -> "Hartree_fock" +| Properties -> "Properties" ;; @@ -88,20 +86,18 @@ let get s = f Ao_basis.(read, to_rst) | Determinants_by_hand -> f Determinants_by_hand.(read, to_rst) - | Hartree_fock -> - f Hartree_fock.(read, to_rst) - | Pseudo -> - f Pseudo.(read, to_rst) - | Integrals_bielec -> - f Integrals_bielec.(read, to_rst) - | Perturbation -> - f Perturbation.(read, to_rst) - | Properties -> - f Properties.(read, to_rst) - | Foboci -> - f Foboci.(read, to_rst) | Determinants -> f Determinants.(read, to_rst) + | Integrals_bielec -> + f Integrals_bielec.(read, to_rst) + | Pseudo -> + f Pseudo.(read, to_rst) + | Perturbation -> + f Perturbation.(read, to_rst) + | Hartree_fock -> + f Hartree_fock.(read, to_rst) + | Properties -> + f Properties.(read, to_rst) end with | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "") @@ -139,13 +135,12 @@ let set str s = in let open Input in match s with - | Hartree_fock -> write Hartree_fock.(of_rst, write) s - | Pseudo -> write Pseudo.(of_rst, write) s - | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s - | Perturbation -> write Perturbation.(of_rst, write) s - | Properties -> write Properties.(of_rst, write) s - | Foboci -> write Foboci.(of_rst, write) s | Determinants -> write Determinants.(of_rst, write) s + | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s + | Pseudo -> write Pseudo.(of_rst, write) s + | Perturbation -> write Perturbation.(of_rst, write) s + | Hartree_fock -> write Hartree_fock.(of_rst, write) s + | Properties -> write Properties.(of_rst, write) s | Electrons -> write Electrons.(of_rst, write) s | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s | Nuclei -> write Nuclei.(of_rst, write) s @@ -193,13 +188,12 @@ let run check_only ezfio_filename = Nuclei ; Ao_basis; Electrons ; - Hartree_fock ; - Pseudo ; - Integrals_bielec ; - Perturbation ; - Properties ; - Foboci ; Determinants ; + Integrals_bielec ; + Pseudo ; + Perturbation ; + Hartree_fock ; + Properties ; Mo_basis; Determinants_by_hand ; ] diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 6752afcb..d278ba13 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -359,7 +359,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin y, & lambda & ) - abort_here = abort_all end diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 866a6e3b..155c52ae 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,3 +1,4 @@ + BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ] END_PROVIDER @@ -47,6 +48,7 @@ endif enddo endif + i_pert = 0 if( i_pert == 1)then pert_determinants(k,i) = i_pert endif @@ -110,7 +112,6 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] &BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref,N_det_ref,N_states) ] implicit none BEGIN_DOC ! Dressing matrix in N_det basis @@ -118,7 +119,6 @@ END_PROVIDER integer :: i,j,m delta_ij = 0.d0 delta_ii = 0.d0 - call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref) double precision :: max_delta double precision :: accu @@ -135,6 +135,7 @@ END_PROVIDER endif enddo enddo + !stop "movais delta" print*,'' print*,'' print*,' = ',accu @@ -159,17 +160,6 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) enddo enddo - - !!!!!!!!!! - do ii = 1, N_det_ref - do jj = 1, N_det_ref - i = idx_ref(ii) - j = idx_ref(jj) - h_matrix_dressed(i,j,istate) += delta_cas(ii,jj,istate) - h_matrix_dressed(j,i,istate) += delta_cas(ii,jj,istate) - end do - end do - !!!!!!!!!!!!! do ii = 1, N_det_ref i =idx_ref(ii) h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) @@ -272,11 +262,14 @@ subroutine diagonalize_CI_dressed ! eigenstates of the CI matrix END_DOC integer :: i,j + double precision, parameter :: speed = 1d0 + do j=1,N_states_diag do i=1,N_det - psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) * speed + psi_coef(i,j) * (1d0 - speed) enddo enddo SOFT_TOUCH psi_coef end + diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 53a9417b..30b161d4 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -1,30 +1,63 @@ use bitmasks + + BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] + use bitmasks + implicit none + +! delta_ij(:,:,:) = delta_ij_old(:,:,:) +! delta_ii(:,:) = delta_ii_old(:,:) + delta_ij(:,:,:) = delta_mrcepa0_ij(:,:,:)! - delta_sub_ij(:,:,:) + delta_ii(:,:)= delta_mrcepa0_ii(:,:)! - delta_sub_ii(:,:) +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_det_non_ref) ] &BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_det_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (2) ] use bitmasks implicit none - integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), active_sorb(2) + integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), nonactive_sorb(2) integer i, II, j, k logical, external :: detEq print *, "provide cepa0" - active_sorb = (/b'001100000000', b'001100000000'/) + active_sorb(:) = 0_8 + nonactive_sorb(:) = not(0_8) + + if(N_det_ref > 1) then + do i=1, N_det_ref + active_sorb(1) = ior(psi_ref(1,1,i), active_sorb(1)) + active_sorb(2) = ior(psi_ref(1,2,i), active_sorb(2)) + nonactive_sorb(1) = iand(psi_ref(1,1,i), nonactive_sorb(1)) + nonactive_sorb(2) = iand(psi_ref(1,2,i), nonactive_sorb(2)) + end do + active_sorb(1) = iand(active_sorb(1), not(nonactive_sorb(1))) + active_sorb(2) = iand(active_sorb(2), not(nonactive_sorb(2))) + end if + do i=1, N_det_non_ref det_noactive(1,1,i) = iand(psi_non_ref(1,1,i), not(active_sorb(1))) det_noactive(1,2,i) = iand(psi_non_ref(1,2,i), not(active_sorb(2))) 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 +! print "(B30,B30)", det_noactive(1,1,i), det_noactive(1,2,i) +! end do +! stop do i=1,N_det_ref det_ref_active(i) = iand(psi_ref(1,1,i), active_sorb(1)) det_ref_active(i) = det_ref_active(i) + iand(psi_ref(1,2,i), active_sorb(2)) * 2_8**32_8 end do - + cepa0_shortcut(0) = 1 cepa0_shortcut(1) = 1 det_cepa0_active(1) = iand(psi_non_ref(1,1,det_cepa0_idx(1)), active_sorb(1)) @@ -39,8 +72,8 @@ use bitmasks cepa0_shortcut(cepa0_shortcut(0)) = i end if end do - cepa0_shortcut(0) += 1 - cepa0_shortcut(cepa0_shortcut(0)) = N_det_non_ref+1 + !cepa0_shortcut(0) += 1 + cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1 END_PROVIDER @@ -50,24 +83,22 @@ END_PROVIDER use bitmasks implicit none integer :: i,j,k - double precision :: Hjk, Hki + double precision :: Hjk, Hki, Hij, mat(2,2) integer i_state + provide lambda_mrcc i_state = 1 - do i=1,N_det_ref do j=1,i 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) - call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,k),N_int,Hki) + call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k) end do delta_cas(j,i,i_state) = delta_cas(i,j,i_state) end do end do - - print *, "mrcepa0_cas_dressing", delta_cas(:,:,1) END_PROVIDER logical function detEq(a,b,Nint) @@ -88,8 +119,9 @@ end function - BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] + + 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 @@ -99,14 +131,14 @@ end function 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 integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ - integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle + integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive integer, allocatable :: idx_sorted_bit(:) integer, external :: get_index_in_psi_det_sorted_bit logical, external :: is_in_wavefunction integer :: II, blok - provide det_cepa0_active + provide det_cepa0_active delta_cas lambda_mrcc if(N_int /= 1) then print *, "mrcepa0 experimental N_int==1" @@ -114,11 +146,11 @@ end function end if i_state = 1 - delta_ii(:,:) = 0 - delta_ij(:,:,:) = 0 - + delta_mrcepa0_ii(:,:) = 0d0 + delta_mrcepa0_ij(:,:,:) = 0d0 + ! do i=1,N_det_ref -! delta_ii(i,i_state) = delta_cas(i,i) +! delta_ii(i,i_state) = delta_cas(i,i,i_state) ! end do provide mo_bielec_integrals_in_map @@ -128,53 +160,370 @@ end function 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 - - !sd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii) + !- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii) do blok=1,cepa0_shortcut(0) do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 do II=1,N_det_ref - made_hole = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) - made_particle = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) - - if(popcnt(made_hole) + popcnt(made_particle) > 2) cycle + made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) + made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) + call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int) + if (degree > 2 .or. popcnt(made_hole) * popcnt(made_particle) /= degree*2) cycle do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP" + !do k=1,N_det_non_ref + if(iand(made_hole, det_cepa0_active(k)) /= 0) cycle + if(iand(made_particle, det_cepa0_active(k)) /= made_particle) cycle + myActive = xor(det_cepa0_active(k), made_hole) + myActive = xor(myActive, made_particle) + if(i==k .and. myActive /= det_ref_active(II)) stop "AAAA" + !if(i==k) print *, "i=k" do J=1,N_det_ref - if(iand(made_hole, det_ref_active(J)) /= made_hole) cycle - if(iand(made_particle, det_ref_active(J)) /= 0) cycle - call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,det_cepa0_idx(k)),N_int,Hki) - contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state) - delta_ij(II, det_cepa0_idx(i), i_state) += contrib + if(det_ref_active(J) /= myActive) cycle -! if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then -! !qs$OMP CRITICAL -! delta_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state) -! !dsd$OMP END CRITICAL -! endif - + !!!!! + call get_excitation_degree(psi_ref(1,1,J),psi_non_ref(1,1,det_cepa0_idx(k)),degree,N_int) + if(degree > 2) stop "BBBB" + !!!!!!!!! +! if(i/=k .and. popcnt(made_hole) /= popcnt(made_particle)) then +! print *, "=================", made_hole, made_particle +! call debug_det(psi_ref(1,1,II),N_int) +! call debug_det(psi_non_ref(1,1,det_cepa0_idx(i)),N_int) +! call debug_det(psi_ref(1,1,J),N_int) +! call debug_det(psi_non_ref(1,1,det_cepa0_idx(k)),N_int) +! print *, "=================" +! end if + + call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,II),N_int,Hki) + + contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state) + delta_mrcepa0_ij(II, det_cepa0_idx(i), i_state) += contrib +! + if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then + !-$OMP CRITICAL + delta_mrcepa0_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state) + !-$OMP END CRITICAL + endif end do end do end do end do end do - !qsd $OMP END PARALLEL DO + !- qs $OMP END PARALLEL DO + !print *, "MMMMMMMMMM ", delta_cas(2,2,i_state) , delta_ii(2,i_state) +! do i=1,N_det_ref +! delta_cas(i,i,i_state) += delta_ii(i,i_state) +! end do + + deallocate(idx_sorted_bit) +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, m, 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, 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 det_cepa0_active delta_cas lambda_mrcc + + if(N_int /= 1) then + print *, "mrsc2 experimental N_int==1" + stop + end if + + i_state = 1 + delta_sub_ij(:,:,:) = 0d0 + delta_sub_ii(:,:) = 0d0 + + provide mo_bielec_integrals_in_map + 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 + + !$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_sub_ij, delta_sub_ii) + do i=1,N_det_non_ref + if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪" + 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) + !call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,i),exc_Ii,degree,phase_Ii,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) + + if(psi_non_ref(1,1,l) /= det_tmp(1,1)) stop "sdf" + call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl) + + do k=1,N_det_non_ref + 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 + + det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1))) + det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1))) + ok = (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2)))) + + det_tmp(1,1) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,k)), not(active_sorb(2))) + det_tmp(1,2) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,i)), not(active_sorb(2))) + ok = ok .and. (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2)))) + + if(ok) cycle + + +! call decode_exc(exc_Ii,degree,h1_,p1_,h2_,p2_,s1_,s2_) +! call decode_exc(exc_Ik,degree2,h1,p1,h2,p2,s1,s2) +! +! +! det_tmp(:,:) = 0_bit_kind +! call set_det_bit(det_tmp, p1, s1) +! call set_det_bit(det_tmp, h1, s1) +! call set_det_bit(det_tmp, p1_, s1_) +! call set_det_bit(det_tmp, h1_, s1_) +! if(degree == 2) then +! call set_det_bit(det_tmp, p2_, s2_) +! call set_det_bit(det_tmp, h2_, s2_) +! end if +! if(degree2 == 2) then +! call set_det_bit(det_tmp, p2, s2) +! call set_det_bit(det_tmp, h2, s2) +! end if +! deg = 0 +! do ni = 1, N_int +! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) +! end do +! if(deg == 2*degree2 + 2*degree) cycle + + + + + +! if(degree == -1) 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) + !$OMP CRITICAL + delta_sub_ij(II, i, i_state) += contrib + if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then + delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) + endif + !$OMP END CRITICAL + end do + end do + end do + end do + !$OMP END PARALLEL DO deallocate(idx_sorted_bit) END_PROVIDER subroutine set_det_bit(det, p, s) - use bitmasks - 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) + 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, delta_ij_old, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_det_ref,N_states) ] +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_, x(2), y(2) + 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 + 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) + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + + delta_ii_old(:,:) = 0 + delta_ij_old(:,:,:) = 0 + + i_state = 1 + provide mo_bielec_integrals_in_map + 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 + + !$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij_old, delta_ii_old) + do i = 1 , N_det_non_ref + if(mod(i,1000) == 0) print *, i, N_det_non_ref + do i_I = 1 , N_det_ref + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int) + if(degree2 == -1) cycle + ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state) + + call decode_exc(exc_iI,degree2,h1,p1,h2,p2,s1,s2) + + call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi) + diI = hIi * lambda_mrcc(i_state,i) + do J = 1 , N_det_ref !!! + call get_excitation(psi_ref(1,1,i_I),psi_ref(1,1,J),exc_IJ,degree,phase_IJ,N_int) + call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi) + delta_JI = hJi * diI + do k = 1 , N_det_non_ref + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int) + if(degree == -1) cycle + + call decode_exc(exc_Ik,degree,h1_,p1_,h2_,p2_,s1_,s2_) + + + det_tmp(:,:) = 0_bit_kind + det_tmp2(:,:) = 0_bit_kind + + !!!!!!!!!!!!!!! + + + + + + det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1))) + det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1))) + ok = (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2)))) + + det_tmp(1,1) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,k)), not(active_sorb(2))) + det_tmp(1,2) = iand(xor(HF_bitmask(1,2), psi_non_ref(1,2,i)), not(active_sorb(2))) + ok = ok .and. (popcnt(det_tmp(1,1)) + popcnt(det_tmp(1,2)) == popcnt(xor(det_tmp(1,1), det_tmp(1,2)))) + if(.not. ok) cycle + !if(ok) cycle + + !!!!!!!!!!!!!! + + + +! call set_det_bit(det_tmp, p1, s1) +! +! call set_det_bit(det_tmp, p1_, s1_) +! +! if(degree == 2) then +! call set_det_bit(det_tmp, p2_, s2_) +! +! end if +! if(degree2 == 2) then +! call set_det_bit(det_tmp, p2, s2) +! end if +! +! x(:) = 0 +! do ni=1,N_int +! x(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1))) +! x(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1))) +! end do +! +! +! !det_tmp(:,:) = 0_bit_kind +! +! call set_det_bit(det_tmp, h1, s1) +! call set_det_bit(det_tmp, h1_, s1_) +! if(degree == 2) then +! call set_det_bit(det_tmp, h2_, s2_) +! end if +! if(degree2 == 2) then +! call set_det_bit(det_tmp, h2, s2) +! end if +! +! y(1) = -x(1) +! y(2) = -x(2) +! do ni=1,N_int +! y(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1))) +! y(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1))) +! end do +! +! ! print *, x, y +! +! if(x(1) * y(1) /= 0) cycle +! if(x(2) * y(2) /= 0) cycle +! +! +! +! deg = 0 +! do ni = 1, N_int +! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) +! end do +! if(deg /= 2*degree2 + 2*degree) cycle + + + + call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) + + call get_excitation(psi_non_ref(1,1,i), det_tmp, exc_Ik, degree, phase_al, N_int) + + + if(.not. ok) cycle + if(is_in_wavefunction(det_tmp, N_int)) cycle + + + call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp,ok,N_int) + if(.not. ok) cycle + + call get_excitation(psi_ref(1,1,J), det_tmp, exc_Ik, degree, phase_Jl, N_int) + + l = get_index_in_psi_det_sorted_bit(det_tmp, N_int) + if(l == 0) cycle + l = idx_sorted_bit(get_index_in_psi_det_sorted_bit(det_tmp, N_int)) + + if(l ==-1) cycle + + + call i_h_j(psi_non_ref(1,1,k), psi_ref(1,1,i_I),N_int,HkI) + dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik + + + !$OMP CRITICAL + contrib = dkI(i_state) * delta_JI + !erro += abs(dkI(i_state) - psi_non_ref_coef(k,i_state) / psi_ref_coef(1,i_state)) + delta_ij_old(i_I,l,i_state) += contrib + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then + delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(k,i_state) + endif + !$OMP END CRITICAL + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + deallocate(idx_sorted_bit) +END_PROVIDER + subroutine apply_excitation(det, exc, res, ok, Nint) diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index 3f892887..87ae25f3 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -55,7 +55,7 @@ subroutine mrcepa0_iterations 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) + print*,'delta_ii(i,1) = ',delta_cas(i,i,1) enddo print*,'------------' enddo