diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 0f0c1ef9..45d71ee3 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -4,6 +4,7 @@ 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 53e3ea59..67dc9501 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -22,6 +22,7 @@ type keyword = | Integrals_bielec | Perturbation | Properties +| Foboci | Determinants ;; @@ -37,6 +38,7 @@ let keyword_to_string = function | Integrals_bielec -> "Integrals_bielec" | Perturbation -> "Perturbation" | Properties -> "Properties" +| Foboci -> "Foboci" | Determinants -> "Determinants" ;; @@ -96,6 +98,8 @@ let get s = f Perturbation.(read, to_rst) | Properties -> f Properties.(read, to_rst) + | Foboci -> + f Foboci.(read, to_rst) | Determinants -> f Determinants.(read, to_rst) end @@ -140,6 +144,7 @@ let set str 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 | Electrons -> write Electrons.(of_rst, write) s | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s @@ -193,6 +198,7 @@ let run check_only ezfio_filename = Integrals_bielec ; Perturbation ; Properties ; + Foboci ; Determinants ; Mo_basis; Determinants_by_hand ; diff --git a/plugins/CAS_SD/H_apply.irp.f b/plugins/CAS_SD/H_apply.irp.f index 35c45fb6..aa393bc7 100644 --- a/plugins/CAS_SD/H_apply.irp.f +++ b/plugins/CAS_SD/H_apply.irp.f @@ -20,18 +20,22 @@ print s s = H_apply("CAS_S",do_double_exc=False) +s.unset_double_excitations() print s s = H_apply("CAS_S_selected_no_skip",do_double_exc=False) +s.unset_double_excitations() s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() print s s = H_apply("CAS_S_selected",do_double_exc=False) +s.unset_double_excitations() s.set_selection_pt2("epstein_nesbet_2x2") print s s = H_apply("CAS_S_PT2",do_double_exc=False) +s.unset_double_excitations() s.set_perturbation("epstein_nesbet_2x2") print s diff --git a/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f b/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f new file mode 100644 index 00000000..2e25f431 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f @@ -0,0 +1,37 @@ +program Dressed_Ref_Hamiltonian implicit none + BEGIN_DOC +! TODO + END_DOC + print *, ' _/ ' + print *, ' -:\_?, _Jm####La ' + print *, 'J"(:" > _]#AZ#Z#UUZ##, ' + print *, '_,::./ %(|i%12XmX1*1XL _?, ' + print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( ' + print *, ' .:< ]J=mQD?WXn|,)nr" ' + print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" ' + print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 ' + print *, ' "23Z#1S2oo2XXSnnnoSo2>v" ' + print *, ' miX#L -~`""!!1}oSoe|i7 ' + print *, ' 4cn#m, v221=|v[ ' + print *, ' ]hI3Zma,;..__wXSe=+vo ' + print *, ' ]Zov*XSUXXZXZXSe||vo2 ' + print *, ' ]Z#>=|< ' + print *, ' -ziiiii||||||+||==+> ' + print *, ' -%|+++||=|=+|=|==/ ' + print *, ' -a>====+|====-:- ' + print *, ' "~,- -- /- ' + print *, ' -. )> ' + print *, ' .~ +- ' + print *, ' . .... : . ' + print *, ' -------~ ' + print *, '' +end diff --git a/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES b/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..55c429bc --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +MRCC_Utils diff --git a/plugins/Dressed_Ref_Hamiltonian/README.rst b/plugins/Dressed_Ref_Hamiltonian/README.rst new file mode 100644 index 00000000..71f2f099 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/README.rst @@ -0,0 +1,16 @@ +======================= +Dressed_Ref_Hamiltonian +======================= +The following modules proposes to build an effective Hamiltonian +spanned on the reference determinants supposed to be the CAS ones. +The effective matrix Hamiltonian are built using the multi parentage +proposal used in the MR-CCSD formalism of Giner et. al. (JCP, 144, 064101 (2016); doi: 10.1063/1.4940781) + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f b/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f new file mode 100644 index 00000000..2e66ff16 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f @@ -0,0 +1,41 @@ + BEGIN_PROVIDER [double precision, psi_ref_coef_dressed, (n_det_ref,N_states) ] +&BEGIN_PROVIDER [double precision, energies_ref_dressed, (N_states) ] + implicit none + integer :: i,j,k,l,istate,igoodstate + double precision, allocatable :: H_matrix_tmp(:,:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:),psi_coef_ref_tmp(:) + double precision :: accu, accu1 + allocate(H_matrix_tmp(n_det_ref,n_det_ref)) + allocate(eigvalues(n_det_ref)) + allocate(eigvectors(n_det_ref,n_det_ref)) + allocate(psi_coef_ref_tmp(n_det_ref)) + do istate = 1, N_states + accu1 = 0.d0 + do j = 1, n_det_ref + accu1 += psi_ref_coef(j,istate)**2 ! norm of the "istate" eigenvector in the projected in the reference space + do k = 1, n_det_ref + H_matrix_tmp(j,k) = hamiltonian_total_dressed(j,k,istate) + enddo + enddo + accu1 = 1.d0/dsqrt(accu1) + do j = 1, n_det_ref + psi_coef_ref_tmp(j) = psi_ref_coef(j,istate) * accu1 + enddo + call lapack_diagd(eigvalues,eigvectors,H_matrix_tmp,n_det_ref,n_det_ref) + do j = 1, n_det_ref + accu = 0.d0 + do k = 1, n_det_ref + accu += eigvectors(k,j) * psi_coef_ref_tmp(k) + enddo + if(dabs(accu).gt.0.9d0)then + igoodstate = j + exit + endif + enddo + energies_ref_dressed(istate) = eigvalues(igoodstate) + do j = 1,n_det_ref + psi_ref_coef_dressed(j,istate) = eigvectors(j,igoodstate) + enddo + enddo + +END_PROVIDER diff --git a/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f b/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f new file mode 100644 index 00000000..90b27e07 --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f @@ -0,0 +1,46 @@ +BEGIN_PROVIDER [double precision, dressing_ref_hamiltonian, (n_det_ref,n_det_ref,N_states)] + implicit none + integer :: i,j,k,l + integer :: ii,jj,istate + double precision :: hij,sec_order,H_ref(N_det_ref),hik,hkl + integer :: idx(0:N_det_ref) + double precision :: accu_negative,accu_positive,phase + integer :: degree_exc_ionic,degree_exc_neutral,exc(0:2,2,2) + dressing_ref_hamiltonian = 0.d0 + accu_negative = 0.d0 + accu_positive = 0.d0 + integer :: h1,p1,h2,p2,s1,s2 + do istate = 1, N_states + do i = 1, N_det_non_ref + call filter_connected_i_H_psi0(psi_ref,psi_non_ref(1,1,i),N_int,N_det_ref,idx) + H_ref = 0.d0 + do ii=1,idx(0) + k = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(psi_ref(1,1,k),psi_non_ref(1,1,i),N_int,hij) + H_ref(k) = hij + enddo + do ii= 1, idx(0) + k = idx(ii) + hik = H_ref(k) * lambda_mrcc(istate,i) + do jj = 1, idx(0) + l = idx(jj) + dressing_ref_hamiltonian(k,l,istate) += hik * H_ref(l) + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, hamiltonian_total_dressed, (n_det_ref,n_det_ref,N_states)] + implicit none + integer :: i,j,k + do k = 1, N_states + do i = 1, N_det_ref + do j = 1, N_det_ref + hamiltonian_total_dressed(j,i,k) = dressing_ref_hamiltonian(j,i,k) + ref_hamiltonian_matrix(j,i) + enddo + enddo + enddo + +END_PROVIDER diff --git a/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f b/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f new file mode 100644 index 00000000..aab716ac --- /dev/null +++ b/plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f @@ -0,0 +1,108 @@ +program print + read_wf = .True. + touch read_wf + call provide_all_stuffs +end +subroutine provide_all_stuffs + implicit none + provide ref_hamiltonian_matrix dressing_ref_hamiltonian + integer :: i,j,istate + double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + double precision, allocatable :: H_dressed(:,:) + double precision, allocatable :: H_print(:,:) + double precision :: accu_norm + allocate (H_dressed(N_det_ref,N_det_ref)) + allocate (H_print(N_det_ref,N_det_ref)) + allocate (psi_restart_ref_normalized(N_det_ref)) + allocate (psi_ref_zeroth_order(N_det_ref)) + print*,'# nuclear_repulsion = ',nuclear_repulsion + allocate (psi_ref_dressed(N_det_ref)) + allocate (eigvalues(N_det_ref)) + allocate (eigvectors(N_det_ref,N_det_ref)) + + + + do istate= 1, N_states + do i = 1, N_det_ref + do j = 1, N_det_ref + H_print(i,j) = ref_hamiltonian_matrix(j,i) + enddo + enddo + do i = 1, N_det_ref + H_print(i,i) -= ref_hamiltonian_matrix(1,1) + enddo + print*,'Ref Hamiltonian matrix emelent = ',ref_hamiltonian_matrix(1,1) + print*,'ISTATE = ',istate + accu_norm = 0.d0 + do i = 1, N_det_ref + accu_norm += psi_ref_coef(i,1) * psi_ref_coef(i,1) + enddo + print*,'accu_norm = ',accu_norm + accu_norm = 1.d0/dsqrt(accu_norm) + do i = 1, N_det_ref + psi_restart_ref_normalized(i) = psi_ref_coef(i,istate)* accu_norm + enddo + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') H_print(i,:) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX DRESSING' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') dressing_ref_hamiltonian(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + do i = 1, N_det_ref + do j = 1, N_det_ref + H_dressed(j,i) = ref_hamiltonian_matrix(j,i) + dressing_ref_hamiltonian(j,i,istate) + H_print(i,j) += dressing_ref_hamiltonian(j,i,istate) + enddo + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'TOTAL DRESSED H MATRIX ' + print*,'' + do i = 1, N_det_ref + write(*,'(10(F8.5 ,4X))') H_print(i,:) + enddo + print*,'' + print*,'' + print*,'' + + + call lapack_diagd(eigvalues,eigvectors,ref_hamiltonian_matrix,n_det_ref,n_det_ref) + do i = 1, N_det_ref + psi_ref_zeroth_order(i) = eigvectors(i,istate) + enddo + + + call lapack_diagd(eigvalues,eigvectors,H_dressed,n_det_ref,n_det_ref) + do i = 1, N_det_ref + psi_ref_dressed(i) = eigvectors(i,istate) + enddo + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = 1, N_det_ref + write(*,'(10(F10.7 ,4X))') psi_ref_coef(i,istate)/psi_ref_coef(1,istate), psi_ref_dressed(i)/psi_ref_dressed(1),psi_ref_zeroth_order(i)/psi_ref_zeroth_order(1) + enddo + enddo + + deallocate (H_dressed) + deallocate (H_print) + deallocate (psi_restart_ref_normalized) + deallocate (psi_ref_zeroth_order) + deallocate (psi_ref_dressed) + + deallocate (eigvalues) + deallocate (eigvectors) + +end diff --git a/plugins/FOBOCI/EZFIO.cfg b/plugins/FOBOCI/EZFIO.cfg new file mode 100644 index 00000000..d4a10add --- /dev/null +++ b/plugins/FOBOCI/EZFIO.cfg @@ -0,0 +1,30 @@ +[threshold_singles] +type: double precision +doc: threshold to select the pertinent single excitations at second order +interface: ezfio,provider,ocaml +default: 0.01 + +[threshold_fobo_dm] +type: double precision +doc: threshold to eliminate small density matrix elements in the fobo procedure +interface: ezfio,provider,ocaml +default: 0.00001 + +[do_it_perturbative] +type: logical +doc: if true, you do the FOBOCI calculation perturbatively +interface: ezfio,provider,ocaml +default: .False. + +[second_order_h] +type: logical +doc: if true, you do the FOBOCI calculation using second order intermediate Hamiltonian +interface: ezfio,provider,ocaml +default: .False. + +[do_all_2p] +type: logical +doc: if true, you do all 2p type excitation on the LMCT +interface: ezfio,provider,ocaml +default: .True. + diff --git a/plugins/FOBOCI/H_apply.irp.f b/plugins/FOBOCI/H_apply.irp.f new file mode 100644 index 00000000..0a488753 --- /dev/null +++ b/plugins/FOBOCI/H_apply.irp.f @@ -0,0 +1,50 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("just_1h_1p") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_only_1h1p() +print s + + +s = H_apply("all_but_1h_and_1p") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_1h() +s.filter_1p() +print s + + + +s = H_apply("standard") +s.set_selection_pt2("epstein_nesbet") +s.unset_skip() +print s + +s = H_apply("just_mono",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +print s + + + +s = H_apply("just_mono_no_1h_no_1p",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_1h() +s.filter_1p() +print s + +s = H_apply("just_mono_no_1h_no_1p_no_2p",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_1h() +s.filter_1p() +s.filter_2p() +print s + + +END_SHELL + diff --git a/plugins/FOBOCI/H_apply_dressed_autonom.irp.f b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f new file mode 100644 index 00000000..c5b0aa5c --- /dev/null +++ b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f @@ -0,0 +1,605 @@ +subroutine H_apply_dressed_pert_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref ) + use omp_lib + use bitmasks + implicit none + BEGIN_DOC + ! Generate all double excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC + integer,parameter :: size_max = 3072 + + integer, intent(in) :: Ndet_generators + double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + + integer ,intent(in) :: i_generator + integer(bit_kind),intent(in) :: key_in(N_int,2) + integer(bit_kind),allocatable :: keys_out(:,:,:) + integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) + integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) + integer, intent(in) :: iproc_in + integer(bit_kind), allocatable :: hole_save(:,:) + integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) + integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) + integer(bit_kind), allocatable :: key_union_hole_part(:) + integer :: ii,i,jj,j,k,ispin,l + integer, allocatable :: occ_particle(:,:), occ_hole(:,:) + integer, allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:) + integer :: kk,pp,other_spin,key_idx + integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) + integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) + + double precision :: mo_bielec_integral + logical :: is_a_two_holes_two_particles + integer, allocatable :: ia_ja_pairs(:,:,:) + integer, allocatable :: ib_jb_pairs(:,:) + double precision :: diag_H_mat_elem + integer :: iproc + integer :: jtest_vvvv + integer(omp_lock_kind), save :: lck, ifirst=0 + if (ifirst == 0) then +!$ call omp_init_lock(lck) + ifirst=1 + endif + + logical :: check_double_excitation + logical :: is_a_1h1p + logical :: b_cycle + check_double_excitation = .True. + iproc = iproc_in + + + + + PROVIDE elec_num_tab +! !$OMP PARALLEL DEFAULT(SHARED) & +! !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, & +! !$OMP occ_particle,occ_hole,j_a,k_a,other_spin, & +! !$OMP hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs, & +! !$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, & +! !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,& +! !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, & +! !$OMP N_elec_in_key_hole_2,ia_ja_pairs,key_union_hole_part) & +! !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & +! !$OMP hole_1, particl_1, hole_2, particl_2, & +! !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc) +!$ iproc = omp_get_thread_num() + allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & + key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& + particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & + occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& + occ_hole_tmp(N_int*bit_kind_size,2),key_union_hole_part(N_int)) + + + + + !!!! First couple hole particle + do j = 1, N_int + hole(j,1) = iand(hole_1(j,1),key_in(j,1)) + hole(j,2) = iand(hole_1(j,2),key_in(j,2)) + particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1)) + particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2)) + enddo + call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int) + call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int) + call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int) + call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int) + allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2), & + ib_jb_pairs(2,0:(elec_alpha_num)*mo_tot_num)) + + do ispin=1,2 + i=0 + do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole + i_a = occ_hole(ii,ispin) + ASSERT (i_a > 0) + ASSERT (i_a <= mo_tot_num) + + do jj=1,N_elec_in_key_part_1(ispin) !particle + j_a = occ_particle(jj,ispin) + ASSERT (j_a > 0) + ASSERT (j_a <= mo_tot_num) + i += 1 + ia_ja_pairs(1,i,ispin) = i_a + ia_ja_pairs(2,i,ispin) = j_a + enddo + enddo + ia_ja_pairs(1,0,ispin) = i + enddo + + key_idx = 0 + + integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b + integer(bit_kind) :: test(N_int,2) + double precision :: accu + logical, allocatable :: array_pairs(:,:) + allocate(array_pairs(mo_tot_num,mo_tot_num)) + accu = 0.d0 + do ispin=1,2 + other_spin = iand(ispin,1)+1 + if (abort_here) then + exit + endif +! !$OMP DO SCHEDULE (static) + do ii=1,ia_ja_pairs(1,0,ispin) + if (abort_here) then + cycle + endif + i_a = ia_ja_pairs(1,ii,ispin) + ASSERT (i_a > 0) + ASSERT (i_a <= mo_tot_num) + j_a = ia_ja_pairs(2,ii,ispin) + ASSERT (j_a > 0) + ASSERT (j_a <= mo_tot_num) + hole = key_in + k = ishft(i_a-1,-bit_kind_shift)+1 + j = i_a-ishft(k-1,bit_kind_shift)-1 + hole(k,ispin) = ibclr(hole(k,ispin),j) + k_a = ishft(j_a-1,-bit_kind_shift)+1 + l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 + hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + + !!!! Second couple hole particle + do j = 1, N_int + hole_tmp(j,1) = iand(hole_2(j,1),hole(j,1)) + hole_tmp(j,2) = iand(hole_2(j,2),hole(j,2)) + particle_tmp(j,1) = iand(xor(particl_2(j,1),hole(j,1)),particl_2(j,1)) + particle_tmp(j,2) = iand(xor(particl_2(j,2),hole(j,2)),particl_2(j,2)) + enddo + + call bitstring_to_list(particle_tmp(1,1),occ_particle_tmp(1,1),N_elec_in_key_part_2(1),N_int) + call bitstring_to_list(particle_tmp(1,2),occ_particle_tmp(1,2),N_elec_in_key_part_2(2),N_int) + call bitstring_to_list(hole_tmp (1,1),occ_hole_tmp (1,1),N_elec_in_key_hole_2(1),N_int) + call bitstring_to_list(hole_tmp (1,2),occ_hole_tmp (1,2),N_elec_in_key_hole_2(2),N_int) + + ! hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin) + hole_save = hole + + ! Build array of the non-zero integrals of second excitation + array_pairs = .True. + if (ispin == 1) then + integer :: jjj + + i=0 + do kk = 1,N_elec_in_key_hole_2(other_spin) + i_b = occ_hole_tmp(kk,other_spin) + ASSERT (i_b > 0) + ASSERT (i_b <= mo_tot_num) + do jjj=1,N_elec_in_key_part_2(other_spin) ! particule + j_b = occ_particle_tmp(jjj,other_spin) + ASSERT (j_b > 0) + ASSERT (j_b <= mo_tot_num) + if (array_pairs(i_b,j_b)) then + + i+= 1 + ib_jb_pairs(1,i) = i_b + ib_jb_pairs(2,i) = j_b + endif + enddo + enddo + ib_jb_pairs(1,0) = i + + do kk = 1,ib_jb_pairs(1,0) + hole = hole_save + i_b = ib_jb_pairs(1,kk) + j_b = ib_jb_pairs(2,kk) + k = ishft(i_b-1,-bit_kind_shift)+1 + j = i_b-ishft(k-1,bit_kind_shift)-1 + hole(k,other_spin) = ibclr(hole(k,other_spin),j) + key = hole + k = ishft(j_b-1,-bit_kind_shift)+1 + l = j_b-ishft(k-1,bit_kind_shift)-1 + key(k,other_spin) = ibset(key(k,other_spin),l) + + + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref) + key_idx = 0 + endif + if (abort_here) then + exit + endif + enddo + endif + + ! does all the mono excitations of the same spin + i=0 + do kk = 1,N_elec_in_key_hole_2(ispin) + i_b = occ_hole_tmp(kk,ispin) + if (i_b <= i_a.or.i_b == j_a) cycle + ASSERT (i_b > 0) + ASSERT (i_b <= mo_tot_num) + do jjj=1,N_elec_in_key_part_2(ispin) ! particule + j_b = occ_particle_tmp(jjj,ispin) + ASSERT (j_b > 0) + ASSERT (j_b <= mo_tot_num) + if (j_b <= j_a) cycle + if (array_pairs(i_b,j_b)) then + + i+= 1 + ib_jb_pairs(1,i) = i_b + ib_jb_pairs(2,i) = j_b + endif + enddo + enddo + ib_jb_pairs(1,0) = i + + do kk = 1,ib_jb_pairs(1,0) + hole = hole_save + i_b = ib_jb_pairs(1,kk) + j_b = ib_jb_pairs(2,kk) + k = ishft(i_b-1,-bit_kind_shift)+1 + j = i_b-ishft(k-1,bit_kind_shift)-1 + hole(k,ispin) = ibclr(hole(k,ispin),j) + key = hole + k = ishft(j_b-1,-bit_kind_shift)+1 + l = j_b-ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibset(key(k,ispin),l) + + + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref) + key_idx = 0 + endif + if (abort_here) then + exit + endif + enddo ! kk + + enddo ! ii +! !$OMP ENDDO NOWAIT + enddo ! ispin + call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref) + + deallocate (ia_ja_pairs, ib_jb_pairs, & + keys_out, hole_save, & + key,hole, particle, hole_tmp,& + particle_tmp, occ_particle, & + occ_hole, occ_particle_tmp,& + occ_hole_tmp,array_pairs,key_union_hole_part) +! !$OMP END PARALLEL + +end + +subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref ) + use omp_lib + use bitmasks + implicit none + BEGIN_DOC + ! Generate all single excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC + integer,parameter :: size_max = 3072 + + integer, intent(in) :: Ndet_generators + double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + + integer ,intent(in) :: i_generator + integer(bit_kind),intent(in) :: key_in(N_int,2) + integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2) + integer, intent(in) :: iproc_in + integer(bit_kind),allocatable :: keys_out(:,:,:) + integer(bit_kind),allocatable :: hole_save(:,:) + integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:) + integer(bit_kind),allocatable :: hole_tmp(:,:), particle_tmp(:,:) + integer(bit_kind),allocatable :: hole_2(:,:), particl_2(:,:) + integer :: ii,i,jj,j,k,ispin,l + integer,allocatable :: occ_particle(:,:), occ_hole(:,:) + integer,allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:) + integer,allocatable :: ib_jb_pairs(:,:) + integer :: kk,pp,other_spin,key_idx + integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) + integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) + logical :: is_a_two_holes_two_particles + integer(bit_kind), allocatable :: key_union_hole_part(:) + + integer, allocatable :: ia_ja_pairs(:,:,:) + logical, allocatable :: array_pairs(:,:) + double precision :: diag_H_mat_elem + integer(omp_lock_kind), save :: lck, ifirst=0 + integer :: iproc + + logical :: check_double_excitation + logical :: is_a_1h1p + logical :: is_a_1h + logical :: is_a_1p + iproc = iproc_in + + check_double_excitation = .True. + + check_double_excitation = .False. + + + + + if (ifirst == 0) then + ifirst=1 +!!$ call omp_init_lock(lck) + endif + + + + PROVIDE elec_num_tab +! !$OMP PARALLEL DEFAULT(SHARED) & +! !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, & +! !$OMP occ_particle,occ_hole,j_a,k_a,other_spin, & +! !$OMP hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs, & +! !$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, & +! !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,& +! !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, & +! !$OMP N_elec_in_key_hole_2,ia_ja_pairs,key_union_hole_part) & +! !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & +! !$OMP hole_1, particl_1, hole_2, particl_2, & +! !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc) +!!$ iproc = omp_get_thread_num() + allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & + key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& + particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & + occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& + occ_hole_tmp(N_int*bit_kind_size,2),key_union_hole_part(N_int)) + + !!!! First couple hole particle + do j = 1, N_int + hole(j,1) = iand(hole_1(j,1),key_in(j,1)) + hole(j,2) = iand(hole_1(j,2),key_in(j,2)) + particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1)) + particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2)) + enddo + + call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int) + call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int) + call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int) + call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int) + allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2)) + + do ispin=1,2 + i=0 + do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole + i_a = occ_hole(ii,ispin) + do jj=1,N_elec_in_key_part_1(ispin) !particule + j_a = occ_particle(jj,ispin) + i += 1 + ia_ja_pairs(1,i,ispin) = i_a + ia_ja_pairs(2,i,ispin) = j_a + enddo + enddo + ia_ja_pairs(1,0,ispin) = i + enddo + + key_idx = 0 + + integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b + integer(bit_kind) :: test(N_int,2) + double precision :: accu + accu = 0.d0 + integer :: jjtest,na,nb + do ispin=1,2 + other_spin = iand(ispin,1)+1 +! !$OMP DO SCHEDULE (static) + do ii=1,ia_ja_pairs(1,0,ispin) + i_a = ia_ja_pairs(1,ii,ispin) + j_a = ia_ja_pairs(2,ii,ispin) + hole = key_in + k = ishft(i_a-1,-bit_kind_shift)+1 + j = i_a-ishft(k-1,bit_kind_shift)-1 + + hole(k,ispin) = ibclr(hole(k,ispin),j) + k_a = ishft(j_a-1,-bit_kind_shift)+1 + l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 + + hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + na = 0 + nb = 0 +! if (is_a_1h(hole)) then +! cycle +! endif +! if (is_a_1p(hole)) then +! cycle +! endif + + + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = hole(k,1) + keys_out(k,2,key_idx) = hole(k,2) + enddo + if (key_idx == size_max) then + call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref) + key_idx = 0 + endif + enddo ! ii +! !$OMP ENDDO NOWAIT + enddo ! ispin + call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref) + + deallocate (ia_ja_pairs, & + keys_out, hole_save, & + key,hole, particle, hole_tmp,& + particle_tmp, occ_particle, & + occ_hole, occ_particle_tmp,& + occ_hole_tmp,key_union_hole_part) +! !$OMP END PARALLEL + + +end + + +subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref) + implicit none + use omp_lib + use bitmasks + BEGIN_DOC + ! Calls H_apply on the HF determinant and selects all connected single and double + ! excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + END_DOC + + + integer, intent(in) :: Ndet_generators + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators),E_ref + double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators) + + + integer :: i_generator, nmax + double precision :: wall_0, wall_1 + integer(omp_lock_kind) :: lck + integer(bit_kind), allocatable :: mask(:,:,:) + integer :: ispin, k + integer :: iproc + + + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map + + nmax = mod( Ndet_generators,nproc ) + + +! !$ call omp_init_lock(lck) + call start_progress(Ndet_generators,'Selection (norm)',0.d0) + + call wall_time(wall_0) + + iproc = 0 + allocate( mask(N_int,2,6) ) + do i_generator=1,nmax + + progress_bar(1) = i_generator + + if (abort_here) then + exit + endif + + + +! ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_det_generators_input(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_det_generators_input(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not(psi_det_generators_input(k,ispin,i_generator)) ) + enddo + enddo + if(.False.)then + call H_apply_dressed_pert_diexc(psi_det_generators_input(1,1,i_generator), & + mask(1,1,d_hole1), mask(1,1,d_part1), & + mask(1,1,d_hole2), mask(1,1,d_part2), & + i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref) + endif + if(.True.)then + call H_apply_dressed_pert_monoexc(psi_det_generators_input(1,1,i_generator), & + mask(1,1,s_hole ), mask(1,1,s_part ), & + i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref) + endif + call wall_time(wall_1) + + if (wall_1 - wall_0 > 2.d0) then + write(output_determinants,*) & + 100.*float(i_generator)/float(Ndet_generators), '% in ', wall_1-wall_0, 's' + wall_0 = wall_1 + endif + enddo + + deallocate( mask ) + +! !$OMP PARALLEL DEFAULT(SHARED) & +! !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) + call wall_time(wall_0) +! !$ iproc = omp_get_thread_num() + allocate( mask(N_int,2,6) ) +! !$OMP DO SCHEDULE(dynamic,1) + do i_generator=nmax+1,Ndet_generators + if (iproc == 0) then + progress_bar(1) = i_generator + endif + if (abort_here) then + cycle + endif + + + + ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_det_generators_input(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_det_generators_input(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_det_generators_input(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not (psi_det_generators_input(k,ispin,i_generator)) ) + enddo + enddo + + if(.False.)then + call H_apply_dressed_pert_diexc(psi_det_generators_input(1,1,i_generator), & + mask(1,1,d_hole1), mask(1,1,d_part1), & + mask(1,1,d_hole2), mask(1,1,d_part2), & + i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref) + endif + if(.True.)then + call H_apply_dressed_pert_monoexc(psi_det_generators_input(1,1,i_generator), & + mask(1,1,s_hole ), mask(1,1,s_part ), & + i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref) + endif +! !$ call omp_set_lock(lck) + call wall_time(wall_1) + + if (wall_1 - wall_0 > 2.d0) then + write(output_determinants,*) & + 100.*float(i_generator)/float(Ndet_generators), '% in ', wall_1-wall_0, 's' + wall_0 = wall_1 + endif +! !$ call omp_unset_lock(lck) + enddo +! !$OMP END DO + deallocate( mask ) +! !$OMP END PARALLEL +! !$ call omp_destroy_lock(lck) + + abort_here = abort_all + call stop_progress + + + + +end + + diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..adeefe99 --- /dev/null +++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Perturbation Generators_restart Selectors_no_sorted diff --git a/plugins/FOBOCI/README.rst b/plugins/FOBOCI/README.rst new file mode 100644 index 00000000..95a09211 --- /dev/null +++ b/plugins/FOBOCI/README.rst @@ -0,0 +1,12 @@ +====== +FOBOCI +====== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f new file mode 100644 index 00000000..e2c4c01e --- /dev/null +++ b/plugins/FOBOCI/all_singles.irp.f @@ -0,0 +1,362 @@ +subroutine all_single + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 1.d-8 + soft_touch selection_criterion + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'Doing all the mono excitations !' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + enddo + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo +! call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_single_no_1h_or_1p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'Doing all the mono excitations !' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_mono_no_1h_no_1p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + enddo + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo +! call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_single_no_1h_or_1p_or_2p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'Doing all the mono excitations !' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_mono_no_1h_no_1p_no_2p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + enddo + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo +! call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end + + +subroutine all_2p + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion + threshold_davidson = 1.d-5 + soft_touch threshold_davidson davidson_criterion + i = 0 + print*,'' + print*,'' + print*,'' + print*,'' + print*,'' + print*,'*****************************' + print*,'Doing all the 2P excitations' + print*,'*****************************' + print*,'' + print*,'' + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + print*,'N_det_generators = ',N_det_generators + pt2=-1.d0 + E_before = ref_bitmask_energy + + print*,'Initial Step ' + print*,'Inital determinants ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo + n_det_max = 100000 + i = 0 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_standard(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + + enddo + print*,'Final Step ' + print*,'N_det = ',N_det + do i = 1, N_states_diag + print*,'' + print*,'i = ',i + print*,'E = ',CI_energy(i) + print*,'S^2 = ',CI_eigenvectors_s2(i) + enddo +! call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_1h_1p_routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision :: E_before + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_1h_1p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + E_before = CI_energy(1) + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + enddo + deallocate(pt2,norm_pert) +end +subroutine all_but_1h_1p_routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision :: E_before + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_all_but_1h_and_1p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + E_before = CI_energy(1) + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + enddo + deallocate(pt2,norm_pert) +end diff --git a/plugins/FOBOCI/all_singles_split.irp.f b/plugins/FOBOCI/all_singles_split.irp.f new file mode 100644 index 00000000..e7b0943f --- /dev/null +++ b/plugins/FOBOCI/all_singles_split.irp.f @@ -0,0 +1,243 @@ +subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,Ndet_generators_input,dressing_matrix) + implicit none + use bitmasks + integer, intent(in) :: Ndet_generators_input + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators_input) + double precision, intent(inout) :: dressing_matrix(Ndet_generators_input,Ndet_generators_input) + double precision, intent(in) :: psi_coef_generators_input(ndet_generators_input,n_states) + integer :: i,i_hole + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + do i = 1, n_inact_orb + i_hole = list_inact(i) + print*,'' + print*,'Doing all the single excitations from the orbital ' + print*,i_hole + print*,'' + print*,'' + threshold_davidson = 1.d-4 + soft_touch threshold_davidson davidson_criterion + call modify_bitmasks_for_hole(i_hole) + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input) + call all_single + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) + enddo + n_det_max_jacobi = 1000 + soft_touch n_det_max_jacobi +end + + +subroutine all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + implicit none + use bitmasks + double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) + integer :: i,i_hole + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + integer :: n_det_1h1p,n_det_2h1p + integer(bit_kind), allocatable :: psi_ref_out(:,:,:) + integer(bit_kind), allocatable :: psi_1h1p(:,:,:) + integer(bit_kind), allocatable :: psi_2h1p(:,:,:) + double precision, allocatable :: psi_ref_coef_out(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_2h1p(:,:) + call all_single_no_1h_or_1p + + threshold_davidson = 1.d-12 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) + allocate(psi_ref_out(N_int,2,N_det_generators)) + allocate(psi_1h1p(N_int,2,n_det_1h1p)) + allocate(psi_2h1p(N_int,2,n_det_2h1p)) + allocate(psi_ref_coef_out(N_det_generators,N_states)) + allocate(psi_coef_1h1p(n_det_1h1p,N_states)) + allocate(psi_coef_2h1p(n_det_2h1p,N_states)) + call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) + call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h1p,psi_coef_1h1p,n_det_1h1p) + call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_2h1p,psi_coef_2h1p,n_det_2h1p) + + deallocate(psi_ref_out) + deallocate(psi_1h1p) + deallocate(psi_2h1p) + deallocate(psi_ref_coef_out) + deallocate(psi_coef_1h1p) + deallocate(psi_coef_2h1p) + +end + + + + +subroutine all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + implicit none + use bitmasks + double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) + integer :: i,i_hole + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + integer :: n_det_1h1p,n_det_2h1p + integer(bit_kind), allocatable :: psi_ref_out(:,:,:) + integer(bit_kind), allocatable :: psi_1h1p(:,:,:) + integer(bit_kind), allocatable :: psi_2h1p(:,:,:) + double precision, allocatable :: psi_ref_coef_out(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_2h1p(:,:) + do i = 1, n_inact_orb + i_hole = list_inact(i) + print*,'' + print*,'Doing all the single excitations from the orbital ' + print*,i_hole + print*,'' + print*,'' + threshold_davidson = 1.d-4 + soft_touch threshold_davidson davidson_criterion + selection_criterion_factor = 1.d-4 + soft_touch selection_criterion_factor selection_criterion selection_criterion_min + call modify_bitmasks_for_hole(i_hole) + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_generators_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators) + call set_psi_det_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators) + call all_single_no_1h_or_1p + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) + allocate(psi_ref_out(N_int,2,N_det_generators)) + allocate(psi_1h1p(N_int,2,n_det_1h1p)) + allocate(psi_2h1p(N_int,2,n_det_2h1p)) + allocate(psi_ref_coef_out(N_det_generators,N_states)) + allocate(psi_coef_1h1p(n_det_1h1p,N_states)) + allocate(psi_coef_2h1p(n_det_2h1p,N_states)) + call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) + call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h1p,psi_coef_1h1p,n_det_1h1p) + call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_2h1p,psi_coef_2h1p,n_det_2h1p) + + deallocate(psi_ref_out) + deallocate(psi_1h1p) + deallocate(psi_2h1p) + deallocate(psi_ref_coef_out) + deallocate(psi_coef_1h1p) + deallocate(psi_coef_2h1p) + enddo + n_det_max_jacobi = 1000 + soft_touch n_det_max_jacobi +end + + +subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + implicit none + use bitmasks + double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) + integer :: i,i_hole + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + integer :: n_det_1h1p,n_det_1h2p + integer(bit_kind), allocatable :: psi_ref_out(:,:,:) + integer(bit_kind), allocatable :: psi_1h1p(:,:,:) + integer(bit_kind), allocatable :: psi_1h2p(:,:,:) + double precision, allocatable :: psi_ref_coef_out(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_1h2p(:,:) + do i = 1, n_inact_orb + i_hole = list_inact(i) + print*,'' + print*,'Doing all the single excitations from the orbital ' + print*,i_hole + print*,'' + print*,'' + threshold_davidson = 1.d-4 + soft_touch threshold_davidson davidson_criterion + selection_criterion_factor = 1.d-4 + soft_touch selection_criterion_factor selection_criterion selection_criterion_min + call modify_bitmasks_for_hole(i_hole) + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_generators_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators) + call set_psi_det_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators) + call all_single_no_1h_or_1p + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) + allocate(psi_ref_out(N_int,2,N_det_generators)) + allocate(psi_1h1p(N_int,2,n_det_1h1p)) + allocate(psi_1h2p(N_int,2,n_det_1h2p)) + allocate(psi_ref_coef_out(N_det_generators,N_states)) + allocate(psi_coef_1h1p(n_det_1h1p,N_states)) + allocate(psi_coef_1h2p(n_det_1h2p,N_states)) + call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) + call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h1p,psi_coef_1h1p,n_det_1h1p) + call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h2p,psi_coef_1h2p,n_det_1h2p) + + deallocate(psi_ref_out) + deallocate(psi_1h1p) + deallocate(psi_1h2p) + deallocate(psi_ref_coef_out) + deallocate(psi_coef_1h1p) + deallocate(psi_coef_1h2p) + enddo + n_det_max_jacobi = 1000 + soft_touch n_det_max_jacobi +end + +subroutine all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + implicit none + use bitmasks + double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) + integer :: i,i_hole + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + integer :: n_det_1h1p,n_det_1h2p + integer(bit_kind), allocatable :: psi_ref_out(:,:,:) + integer(bit_kind), allocatable :: psi_1h1p(:,:,:) + integer(bit_kind), allocatable :: psi_1h2p(:,:,:) + double precision, allocatable :: psi_ref_coef_out(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_1h2p(:,:) + call all_single_no_1h_or_1p_or_2p + + threshold_davidson = 1.d-12 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) + allocate(psi_ref_out(N_int,2,N_det_generators)) + allocate(psi_1h1p(N_int,2,n_det_1h1p)) + allocate(psi_1h2p(N_int,2,n_det_1h2p)) + allocate(psi_ref_coef_out(N_det_generators,N_states)) + allocate(psi_coef_1h1p(n_det_1h1p,N_states)) + allocate(psi_coef_1h2p(n_det_1h2p,N_states)) + call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) + call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h1p,psi_coef_1h1p,n_det_1h1p) + call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_1h2p,psi_coef_1h2p,n_det_1h2p) + + deallocate(psi_ref_out) + deallocate(psi_1h1p) + deallocate(psi_1h2p) + deallocate(psi_ref_coef_out) + deallocate(psi_coef_1h1p) + deallocate(psi_coef_1h2p) + +end + + diff --git a/plugins/FOBOCI/create_1h_or_1p.irp.f b/plugins/FOBOCI/create_1h_or_1p.irp.f new file mode 100644 index 00000000..140ed504 --- /dev/null +++ b/plugins/FOBOCI/create_1h_or_1p.irp.f @@ -0,0 +1,218 @@ +subroutine create_restart_and_1h(i_hole) + implicit none + use bitmasks + integer, intent(in) :: i_hole + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,i_part_act,ispin,k,l,i_ok + integer :: n_new_det + integer(bit_kind), allocatable :: new_det(:,:,:) + integer(bit_kind), allocatable :: old_psi_det(:,:,:) + allocate (old_psi_det(N_int,2,n_det)) + do i = 1, N_det + do j = 1, N_int + old_psi_det(j,1,i) = psi_det(j,1,i) + old_psi_det(j,2,i) = psi_det(j,2,i) + enddo + enddo + n_new_det = 0 + do j = 1, n_act_orb + i_part_act = list_act(j) ! index of the particle in the active space + do i = 1, N_det + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + enddo + enddo + enddo + + integer :: N_det_old + N_det_old = N_det + N_det += n_new_det + allocate (new_det(N_int,2,n_new_det)) + if (psi_det_size < N_det) then + psi_det_size = N_det + TOUCH psi_det_size + endif + do i = 1, N_det_old + do k = 1, N_int + psi_det(k,1,i) = old_psi_det(k,1,i) + psi_det(k,2,i) = old_psi_det(k,2,i) + enddo + enddo + + n_new_det = 0 + do j = 1, n_act_orb + i_part_act = list_act(j) ! index of the particle in the active space + do i = 1, N_det_old + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + do k = 1, N_int + psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1) + psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2) + enddo + psi_coef(n_det_old+n_new_det,:) = 0.d0 + enddo + enddo + enddo + + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) +end + +subroutine create_restart_and_1p(i_particle) + implicit none + integer, intent(in) :: i_particle + use bitmasks + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,i_hole_act,ispin,k,l,i_ok + integer :: n_new_det + integer(bit_kind), allocatable :: new_det(:,:,:) + integer(bit_kind), allocatable :: old_psi_det(:,:,:) + allocate (old_psi_det(N_int,2,n_det)) + do i = 1, N_det + do j = 1, N_int + old_psi_det(j,1,i) = psi_det(j,1,i) + old_psi_det(j,2,i) = psi_det(j,2,i) + enddo + enddo + n_new_det = 0 + do j = 1, n_act_orb + i_hole_act = list_act(j) ! index of the particle in the active space + do i = 1, N_det + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole_act,i_particle,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + enddo + enddo + enddo + + integer :: N_det_old + N_det_old = N_det + N_det += n_new_det + allocate (new_det(N_int,2,n_new_det)) + if (psi_det_size < N_det) then + psi_det_size = N_det + TOUCH psi_det_size + endif + do i = 1, N_det_old + do k = 1, N_int + psi_det(k,1,i) = old_psi_det(k,1,i) + psi_det(k,2,i) = old_psi_det(k,2,i) + enddo + enddo + + n_new_det = 0 + do j = 1, n_act_orb + i_hole_act = list_act(j) ! index of the particle in the active space + do i = 1, N_det_old + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole_act,i_particle,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + do k = 1, N_int + psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1) + psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2) + enddo + psi_coef(n_det_old+n_new_det,:) = 0.d0 + enddo + enddo + enddo + + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) +end + +subroutine create_restart_1h_1p(i_hole,i_part) + implicit none + use bitmasks + integer, intent(in) :: i_hole + integer, intent(in) :: i_part + + integer :: i,j,i_part_act,ispin,k,l,i_ok + integer(bit_kind) :: key_tmp(N_int,2) + integer :: n_new_det + integer(bit_kind), allocatable :: new_det(:,:,:) + integer(bit_kind), allocatable :: old_psi_det(:,:,:) + + allocate (old_psi_det(N_int,2,n_det)) + do i = 1, N_det + do j = 1, N_int + old_psi_det(j,1,i) = psi_det(j,1,i) + old_psi_det(j,2,i) = psi_det(j,2,i) + enddo + enddo + n_new_det = 0 + i_part_act = i_part ! index of the particle in the active space + do i = 1, N_det + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + enddo + enddo + + integer :: N_det_old + N_det_old = N_det + N_det += n_new_det + allocate (new_det(N_int,2,n_new_det)) + if (psi_det_size < N_det) then + psi_det_size = N_det + TOUCH psi_det_size + endif + do i = 1, N_det_old + do k = 1, N_int + psi_det(k,1,i) = old_psi_det(k,1,i) + psi_det(k,2,i) = old_psi_det(k,2,i) + enddo + enddo + + n_new_det = 0 + i_part_act = i_part ! index of the particle in the active space + do i = 1, N_det_old + do ispin = 1,2 + do k = 1, N_int + key_tmp(k,1) = psi_det(k,1,i) + key_tmp(k,2) = psi_det(k,2,i) + enddo + call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok) + if(i_ok .ne. 1)cycle + n_new_det +=1 + do k = 1, N_int + psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1) + psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2) + enddo + psi_coef(n_det_old+n_new_det,:) = 0.d0 + enddo + enddo + + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) + +end diff --git a/plugins/FOBOCI/density_matrix.irp.f b/plugins/FOBOCI/density_matrix.irp.f new file mode 100644 index 00000000..aaf80c4f --- /dev/null +++ b/plugins/FOBOCI/density_matrix.irp.f @@ -0,0 +1,133 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_generators_restart, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_generators_restart, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, norm_generators_restart] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for the generators restart + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + integer :: degree_respect_to_HF_k + integer :: degree_respect_to_HF_l,index_ref_generators_restart + double precision :: inv_coef_ref_generators_restart + integer :: i + + do i = 1, N_det_generators_restart + ! Find the reference determinant for intermediate normalization + call get_excitation_degree(ref_generators_restart,psi_det_generators_restart(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart = i + inv_coef_ref_generators_restart = 1.d0/psi_coef_generators_restart(i,1) + exit + endif + enddo + norm_generators_restart = 0.d0 + do i = 1, N_det_generators_restart + psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart + norm_generators_restart += psi_coef_generators_restart(i,1)**2 + enddo + + + one_body_dm_mo_alpha_generators_restart = 0.d0 + one_body_dm_mo_beta_generators_restart = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP tmp_a, tmp_b, n_occ_alpha)& + !$OMP SHARED(psi_det_generators_restart,psi_coef_generators_restart,N_int,elec_alpha_num,& + !$OMP elec_beta_num,one_body_dm_mo_alpha_generators_restart,one_body_dm_mo_beta_generators_restart,N_det_generators_restart,mo_tot_num_align,& + !$OMP mo_tot_num,N_states, state_average_weight) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do k=1,N_det_generators_restart + call bitstring_to_list(psi_det_generators_restart(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(psi_det_generators_restart(1,2,k), occ(1,2), n_occ_alpha, N_int) + do m=1,N_states + ck = psi_coef_generators_restart(k,m)*psi_coef_generators_restart(k,m) * state_average_weight(m) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + enddo + do l=1,k-1 + call get_excitation_degree(psi_det_generators_restart(1,1,k),psi_det_generators_restart(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det_generators_restart(1,1,k),psi_det_generators_restart(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef_generators_restart(k,m) * psi_coef_generators_restart(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_dm_mo_alpha_generators_restart = one_body_dm_mo_alpha_generators_restart + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_dm_mo_beta_generators_restart = one_body_dm_mo_beta_generators_restart + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL + + do i = 1, mo_tot_num + print*,'DM restat',i,one_body_dm_mo_beta_generators_restart(i,i) + one_body_dm_mo_alpha_generators_restart(i,i) + enddo + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, one_body_dm_mo_generators_restart, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix for the generators_restart + END_DOC + one_body_dm_mo_generators_restart = one_body_dm_mo_alpha_generators_restart + one_body_dm_mo_beta_generators_restart +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_spin_density_mo_generators_restart, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! rho(alpha) - rho(beta) + END_DOC + one_body_spin_density_mo_generators_restart = one_body_dm_mo_alpha_generators_restart - one_body_dm_mo_beta_generators_restart +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_osoci, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_osoci, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix that will be used for the OSOCI approach + END_DOC +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_1h1p, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_1h1p, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix that will be used for the 1h1p approach + END_DOC +END_PROVIDER + diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f new file mode 100644 index 00000000..a4c6b652 --- /dev/null +++ b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f @@ -0,0 +1,35 @@ +subroutine diag_inactive_virt_and_update_mos + implicit none + integer :: i,j,i_inact,j_inact,i_virt,j_virt + double precision :: tmp(mo_tot_num_align,mo_tot_num) + character*(64) :: label + tmp = 0.d0 + do i = 1, mo_tot_num + tmp(i,i) = Fock_matrix_mo(i,i) + enddo + + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = i+1, n_inact_orb + j_inact = list_inact(j) + tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact) + tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact) + enddo + enddo + + do i = 1, n_virt_orb + i_virt = list_virt(i) + do j = i+1, n_virt_orb + j_virt = list_virt(j) + tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt) + tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt) + enddo + enddo + + + label = "Canonical" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + soft_touch mo_coef + + +end diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f new file mode 100644 index 00000000..2f662f4d --- /dev/null +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -0,0 +1,358 @@ + +subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_generator,n_selected,det_buffer,Nint,iproc,psi_det_generators_input,E_ref) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_generators,size_buffer + double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,size_buffer) + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + integer :: i,j,k,m + integer :: new_size + integer :: degree(Ndet_generators) + integer :: idx(0:Ndet_generators) + logical :: good + + integer :: c_ref + integer :: connected_to_ref + + + double precision :: hka, haa + double precision :: haj + double precision :: f + integer :: connected_to_ref_by_mono + logical :: is_in_wavefunction + double precision :: H_array(Ndet_generators) + double precision :: H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) + double precision :: eigenvectors(Ndet_generators+1,Ndet_generators+1), eigenvalues(Ndet_generators+1) + double precision :: contrib,lambda_i,accu + + do k = 1, Ndet_generators + call i_h_j(psi_det_generators_input(1,1,k),psi_det_generators_input(1,1,k),Nint,hka) + H_matrix_tmp(k,k) = hka + do j = k+1, Ndet_generators + call i_h_j(psi_det_generators_input(1,1,k),psi_det_generators_input(1,1,j),Nint,hka) + H_matrix_tmp(k,j) = hka + H_matrix_tmp(j,k) = hka + enddo + H_matrix_tmp(k,Ndet_generators+1) = 0.d0 + enddo + + do i=1,n_selected + c_ref = connected_to_ref_by_mono(det_buffer(1,1,i),psi_det_generators_input,N_int,i_generator,Ndet_generators) + if (c_ref /= 0) then + cycle + endif + if (is_in_wavefunction(det_buffer(1,1,i),Nint)) then + cycle + endif + call get_excitation_degree_vector(psi_det_generators_input,det_buffer(1,1,i),degree,N_int,Ndet_generators,idx) + H_array = 0.d0 + do k=1,idx(0) + call i_h_j(det_buffer(1,1,i),psi_det_generators_input(1,1,idx(k)),Nint,hka) + H_array(idx(k)) = hka + enddo + + call i_h_j(det_buffer(1,1,i),det_buffer(1,1,i),Nint,haa) + f = 1.d0/(E_ref-haa) + + if(second_order_h)then + lambda_i = f + else + ! You write the new Hamiltonian matrix + do k = 1, Ndet_generators + H_matrix_tmp(k,Ndet_generators+1) = H_array(k) + H_matrix_tmp(Ndet_generators+1,k) = H_array(k) + enddo + H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa + ! Then diagonalize it + call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1) + ! Then you extract the effective denominator + accu = 0.d0 + do k = 1, Ndet_generators + accu += eigenvectors(k,1) * H_array(k) + enddo + lambda_i = eigenvectors(Ndet_generators+1,1)/accu + endif + do k=1,idx(0) + contrib = H_array(idx(k)) * H_array(idx(k)) * lambda_i + delta_ij_generators_(idx(k), idx(k)) += contrib + do j=k+1,idx(0) + contrib = H_array(idx(k)) * H_array(idx(j)) * lambda_i + delta_ij_generators_(idx(k), idx(j)) += contrib + delta_ij_generators_(idx(j), idx(k)) += contrib + enddo + enddo +! H_matrix_tmp_bis(idx(k),idx(k)) += contrib +! H_matrix_tmp_bis(idx(k),idx(j)) += contrib +! H_matrix_tmp_bis(idx(j),idx(k)) += contrib +! do k = 1, Ndet_generators +! do j = 1, Ndet_generators +! H_matrix_tmp_bis(k,j) = H_matrix_tmp(k,j) +! enddo +! enddo +! double precision :: H_matrix_tmp_bis(Ndet_generators,Ndet_generators) +! double precision :: eigenvectors_bis(Ndet_generators,Ndet_generators), eigenvalues_bis(Ndet_generators) +! call lapack_diag(eigenvalues_bis,eigenvectors_bis,H_matrix_tmp_bis,Ndet_generators,Ndet_generators) +! print*,'f,lambda_i = ',f,lambda_i +! print*,'eigenvalues_bi(1)',eigenvalues_bis(1) +! print*,'eigenvalues ',eigenvalues(1) +! do k = 1, Ndet_generators +! print*,'coef,coef_dres = ', eigenvectors(k,1), eigenvectors_bis(k,1) +! enddo +! pause +! accu = 0.d0 +! do k = 1, Ndet_generators +! do j = 1, Ndet_generators +! accu += eigenvectors(k,1) * eigenvectors(j,1) * (H_matrix_tmp(k,j) + delta_ij_generators_(k,j)) +! enddo +! enddo +! print*,'accu,eigv = ',accu,eigenvalues(1) +! pause + + enddo +end + + +subroutine is_a_good_candidate(threshold,is_ok,verbose) + use bitmasks + implicit none + double precision, intent(in) :: threshold + logical, intent(out) :: is_ok + logical, intent(in) :: verbose + + integer :: l,k,m + double precision,allocatable :: dressed_H_matrix(:,:) + double precision,allocatable :: psi_coef_diagonalized_tmp(:,:) + integer(bit_kind), allocatable :: psi_det_generators_input(:,:,:) + + allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators)) + allocate(psi_coef_diagonalized_tmp(N_det_generators,N_states)) + dressed_H_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_int + psi_det_generators_input(l,1,k) = psi_det_generators(l,1,k) + psi_det_generators_input(l,2,k) = psi_det_generators(l,2,k) + enddo + enddo +!call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input) + call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose) + if(do_it_perturbative)then + if(is_ok)then + N_det = N_det_generators + do m = 1, N_states + do k = 1, N_det_generators + do l = 1, N_int + psi_det(l,1,k) = psi_det_generators_input(l,1,k) + psi_det(l,2,k) = psi_det_generators_input(l,2,k) + enddo + psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m) + enddo + enddo + touch psi_coef psi_det N_det + endif + endif + + deallocate(psi_det_generators_input,dressed_H_matrix,psi_coef_diagonalized_tmp) + + + + +end + +subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose) + use bitmasks + implicit none + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + integer, intent(in) :: Ndet_generators + double precision, intent(in) :: threshold + logical, intent(in) :: verbose + logical, intent(out) :: is_ok + double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states) + double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators) + + + integer :: i,j,degree,index_ref_generators_restart,i_count,k,i_det_no_ref + double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij + double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average + logical :: is_a_ref_det(Ndet_generators) + + is_a_ref_det = .False. + do i = 1, N_det_generators + do j = 1, N_det_generators_restart + call get_excitation_degree(psi_det_generators_input(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int) + if(degree == 0)then + is_a_ref_det(i) = .True. + exit + endif + enddo + enddo + + + do i = 1, Ndet_generators + call get_excitation_degree(ref_generators_restart,psi_det_generators_input(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart = i + endif + do j = 1, Ndet_generators + call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij) ! Fill the zeroth order H matrix + dressed_H_matrix(i,j) = hij + enddo + enddo + i_det_no_ref = 0 + diag_h_mat_average = 0.d0 + do i = 1, Ndet_generators + if(is_a_ref_det(i))cycle + i_det_no_ref +=1 + diag_h_mat_average+=dressed_H_matrix(i,i) + enddo + diag_h_mat_average = diag_h_mat_average/dble(i_det_no_ref) + print*,'diag_h_mat_average = ',diag_h_mat_average + print*,'ref h_mat = ',dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) + integer :: number_of_particles, number_of_holes + ! Filter the the MLCT that are higher than 27.2 eV in energy with respect to the reference determinant + do i = 1, Ndet_generators + if(is_a_ref_det(i))cycle + if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then + if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then + is_ok = .False. + return + endif + endif + + ! Filter the the LMCT that are higher than 54.4 eV in energy with respect to the reference determinant + if(number_of_holes(psi_det_generators_input(1,1,i)).eq.1 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.0)then + if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then + is_ok = .False. + return + endif + endif + exit + enddo + + call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix + + double precision :: s2,E_ref(N_states) + integer :: i_state(N_states) + integer :: n_state_good + n_state_good = 0 + if(s2_eig)then + do i = 1, Ndet_generators + call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2) + print*,'s2 = ',s2 + print*,dabs(s2-expected_s2) + if(dabs(s2-expected_s2).le.0.3d0)then + n_state_good +=1 + i_state(n_state_good) = i + E_ref(n_state_good) = eigvalues(i) + endif + if(n_state_good==N_states)then + exit + endif + enddo + else + do i = 1, N_states + i_state(i) = i + E_ref(i) = eigvalues(i) + enddo + endif + do i = 1,N_states + print*,'i_state = ',i_state(i) + enddo + do k = 1, N_states + print*,'state ',k + do i = 1, Ndet_generators + psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k)) + psi_coef_ref(i,k) = eigvectors(i,i_state(k)) + print*,'psi_coef_ref(i) = ',psi_coef_ref(i,k) + enddo + enddo + if(verbose)then + print*,'Zeroth order space :' + do i = 1, Ndet_generators + write(*,'(10(F16.8),X)')dressed_H_matrix(i,:) + enddo + print*,'' + print*,'Zeroth order space Diagonalized :' + do k = 1, N_states + print*,'state ',k + do i = 1, Ndet_generators + print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i) + enddo + enddo + endif + double precision :: E_ref_average + E_ref_average = 0.d0 + do i = 1, N_states + E_ref_average += E_ref(i) + enddo + E_ref_average = E_ref_average / dble(N_states) + + call H_apply_dressed_pert(dressed_H_matrix,Ndet_generators,psi_det_generators_input,E_ref_average) ! Calculate the dressing of the H matrix + if(verbose)then + print*,'Zeroth order space Dressed by outer space:' + do i = 1, Ndet_generators + write(*,'(10(F16.8),X)')dressed_H_matrix(i,:) + enddo + endif + call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix + integer :: i_good_state(0:N_states) + i_good_state(0) = 0 + do i = 1, Ndet_generators + call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2) + ! State following + do k = 1, N_states + accu = 0.d0 + do j =1, Ndet_generators + accu += eigvectors(j,i) * psi_coef_ref(j,k) + enddo + if(dabs(accu).ge.0.8d0)then + i_good_state(0) +=1 + i_good_state(i_good_state(0)) = i + endif + enddo + if(i_good_state(0)==N_states)then + exit + endif + enddo + do i = 1, N_states + i_state(i) = i_good_state(i) + E_ref(i) = eigvalues(i_good_state(i)) + enddo + double precision :: accu + accu = 0.d0 + do k = 1, N_states + do i = 1, Ndet_generators + psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k)) + enddo + enddo + if(verbose)then + do k = 1, N_states + print*,'state ',k + do i = 1, Ndet_generators + print*,'coef, = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i) + enddo + enddo + endif + is_ok = .False. + do i = 1, Ndet_generators + if(is_a_ref_det(i))cycle + do k = 1, N_states + if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold)then + is_ok = .True. + exit + endif + enddo + if(is_ok)then + exit + endif + enddo + if(verbose)then + print*,'is_ok = ',is_ok + endif + + +end + diff --git a/plugins/FOBOCI/fobo_coupled_ci.irp.f b/plugins/FOBOCI/fobo_coupled_ci.irp.f new file mode 100644 index 00000000..29513f25 --- /dev/null +++ b/plugins/FOBOCI/fobo_coupled_ci.irp.f @@ -0,0 +1,5 @@ +program osoci_program +implicit none + call new_approach +! call save_natural_mos +end diff --git a/plugins/FOBOCI/fobo_diff_dm.irp.f b/plugins/FOBOCI/fobo_diff_dm.irp.f new file mode 100644 index 00000000..b0368007 --- /dev/null +++ b/plugins/FOBOCI/fobo_diff_dm.irp.f @@ -0,0 +1,18 @@ +program osoci_program +call debug_det(ref_bitmask,N_int) + +implicit none + call FOBOCI_lmct_mlct_old_thr + call provide_all_the_rest +end +subroutine provide_all_the_rest +implicit none +integer :: i + call update_one_body_dm_mo + call provide_properties + call save_osoci_natural_mos + call save_mos + + + +end diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f new file mode 100644 index 00000000..087f791b --- /dev/null +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -0,0 +1,315 @@ + +subroutine FOBOCI_lmct_mlct_old_thr + use bitmasks + implicit none + integer :: i,j,k,l + integer(bit_kind),allocatable :: unpaired_bitmask(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + double precision :: norm_tmp(N_states),norm_total(N_states) + logical :: test_sym + double precision :: thr,hij + double precision :: threshold + double precision, allocatable :: dressing_matrix(:,:) + logical :: verbose,is_ok + verbose = .True. + threshold = threshold_singles + print*,'threshold = ',threshold + thr = 1.d-12 + allocate(unpaired_bitmask(N_int,2)) + allocate (occ(N_int*bit_kind_size,2)) + do i = 1, N_int + unpaired_bitmask(i,1) = unpaired_alpha_electrons(i) + unpaired_bitmask(i,2) = unpaired_alpha_electrons(i) + enddo + norm_total = 0.d0 + call initialize_density_matrix_osoci + call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int) + print*,'' + print*,'' + print*,'mulliken spin population analysis' + accu =0.d0 + do i = 1, nucl_num + accu += mulliken_spin_densities(i) + print*,i,nucl_charge(i),mulliken_spin_densities(i) + enddo + print*,'' + print*,'' + print*,'DOING FIRST LMCT !!' + do i = 1, n_inact_orb + integer :: i_hole_osoci + i_hole_osoci = list_inact(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_hole_osoci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_hole_osoci = ',i_hole_osoci + call create_restart_and_1h(i_hole_osoci) + call set_generators_to_psi_det + print*,'Passed set generators' + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + if(.not.is_ok)cycle + ! so all the mono excitation on the new generators + allocate(dressing_matrix(N_det_generators,N_det_generators)) + if(.not.do_it_perturbative)then +! call all_single + dressing_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) + dressing_matrix(k,l) = hkl + enddo + enddo + double precision :: hkl +! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) +! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) + call debug_det(reunion_of_bitmask,N_int) + call all_single + endif + call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) + do k = 1, N_states + print*,'norm_tmp = ',norm_tmp(k) + norm_total(k) += norm_tmp(k) + enddo + call update_density_matrix_osoci + deallocate(dressing_matrix) + enddo + + if(.True.)then + print*,'' + print*,'DOING THEN THE MLCT !!' + do i = 1, n_virt_orb + integer :: i_particl_osoci + i_particl_osoci = list_virt(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_particl_osoci= ',i_particl_osoci + ! Initialize the bitmask to the restart ones + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call modify_bitmasks_for_particl(i_particl_osoci) + call print_generators_bitmasks_holes + ! Impose that only the active part can be reached + call set_bitmask_hole_as_input(unpaired_bitmask) +!! call all_single_h_core + call create_restart_and_1p(i_particl_osoci) +!! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) +!! ! so all the mono excitation on the new generators + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + if(.not.is_ok)cycle + allocate(dressing_matrix(N_det_generators,N_det_generators)) + if(.not.do_it_perturbative)then + dressing_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) + dressing_matrix(k,l) = hkl + enddo + enddo + ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) + ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) + call all_single + endif + call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) + do k = 1, N_states + print*,'norm_tmp = ',norm_tmp(k) + norm_total(k) += norm_tmp(k) + enddo + call update_density_matrix_osoci + deallocate(dressing_matrix) + enddo + endif + if(.False.)then + print*,'LAST loop for all the 1h-1p' + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call set_bitmask_particl_as_input(inact_virt_bitmask) + call set_bitmask_hole_as_input(inact_virt_bitmask) +! call set_bitmask_particl_as_input(reunion_of_bitmask) +! call set_bitmask_hole_as_input(reunion_of_bitmask) + call all_single + call set_intermediate_normalization_1h1p(norm_tmp) + norm_total += norm_tmp + call update_density_matrix_osoci + endif + + + print*,'norm_total = ',norm_total + norm_total = norm_generators_restart + norm_total = 1.d0/norm_total +! call rescale_density_matrix_osoci(norm_total) + double precision :: accu + accu = 0.d0 + do i = 1, mo_tot_num + accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i) + enddo + print*,'accu = ',accu +end + + +subroutine FOBOCI_mlct_old + use bitmasks + implicit none + integer :: i,j,k,l + integer(bit_kind),allocatable :: unpaired_bitmask(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + double precision :: norm_tmp,norm_total + logical :: test_sym + double precision :: thr + double precision :: threshold + logical :: verbose,is_ok + verbose = .False. + threshold = 1.d-2 + thr = 1.d-12 + allocate(unpaired_bitmask(N_int,2)) + allocate (occ(N_int*bit_kind_size,2)) + do i = 1, N_int + unpaired_bitmask(i,1) = unpaired_alpha_electrons(i) + unpaired_bitmask(i,2) = unpaired_alpha_electrons(i) + enddo + norm_total = 0.d0 + call initialize_density_matrix_osoci + call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int) + print*,'' + print*,'' + print*,'' + print*,'DOING FIRST MLCT !!' + do i = 1, n_virt_orb + integer :: i_particl_osoci + i_particl_osoci = list_virt(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_particl_osoci= ',i_particl_osoci + ! Initialize the bitmask to the restart ones + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call modify_bitmasks_for_particl(i_particl_osoci) + call print_generators_bitmasks_holes + ! Impose that only the active part can be reached + call set_bitmask_hole_as_input(unpaired_bitmask) +! call all_single_h_core + call create_restart_and_1p(i_particl_osoci) +! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) +! ! so all the mono excitation on the new generators + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + is_ok =.True. + if(.not.is_ok)cycle + call all_single + call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) + print*,'norm_tmp = ',norm_tmp + norm_total += norm_tmp + call update_density_matrix_osoci + enddo + + print*,'norm_total = ',norm_total + norm_total += 1.d0 + norm_total = 1.d0/norm_total + call rescale_density_matrix_osoci(norm_total) + double precision :: accu + accu = 0.d0 + do i = 1, mo_tot_num + accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i) + enddo + print*,'accu = ',accu +end + + +subroutine FOBOCI_lmct_old + use bitmasks + implicit none + integer :: i,j,k,l + integer(bit_kind),allocatable :: unpaired_bitmask(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + double precision :: norm_tmp,norm_total + logical :: test_sym + double precision :: thr + double precision :: threshold + logical :: verbose,is_ok + verbose = .False. + threshold = 1.d-2 + thr = 1.d-12 + allocate(unpaired_bitmask(N_int,2)) + allocate (occ(N_int*bit_kind_size,2)) + do i = 1, N_int + unpaired_bitmask(i,1) = unpaired_alpha_electrons(i) + unpaired_bitmask(i,2) = unpaired_alpha_electrons(i) + enddo + norm_total = 0.d0 + call initialize_density_matrix_osoci + call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int) + print*,'' + print*,'' + print*,'DOING FIRST LMCT !!' + do i = 1, n_inact_orb + integer :: i_hole_osoci + i_hole_osoci = list_inact(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_hole_osoci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_hole_osoci = ',i_hole_osoci + ! Initialize the bitmask to the restart ones + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call modify_bitmasks_for_hole(i_hole_osoci) + call print_generators_bitmasks_holes + ! Impose that only the active part can be reached + call set_bitmask_particl_as_input(unpaired_bitmask) +! call all_single_h_core + call create_restart_and_1h(i_hole_osoci) +! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + if(.not.is_ok)cycle +! ! so all the mono excitation on the new generators + call all_single +! call set_intermediate_normalization_lmct_bis(norm_tmp,i_hole_osoci) + call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) + print*,'norm_tmp = ',norm_tmp + norm_total += norm_tmp + call update_density_matrix_osoci + enddo + + print*,'norm_total = ',norm_total + norm_total += 1.d0 + norm_total = 1.d0/norm_total + call rescale_density_matrix_osoci(norm_total) + double precision :: accu + accu = 0.d0 + do i = 1, mo_tot_num + accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i) + enddo + print*,'accu = ',accu +end diff --git a/plugins/FOBOCI/generators_restart_save.irp.f b/plugins/FOBOCI/generators_restart_save.irp.f new file mode 100644 index 00000000..dca4c901 --- /dev/null +++ b/plugins/FOBOCI/generators_restart_save.irp.f @@ -0,0 +1,126 @@ +use bitmasks + +BEGIN_PROVIDER [ integer, N_det_generators_restart ] + implicit none + BEGIN_DOC + ! Number of determinants in the wave function + END_DOC + logical :: exists + character*64 :: label + integer, save :: ifirst = 0 +!if(ifirst == 0)then + PROVIDE ezfio_filename + call ezfio_has_determinants_n_det(exists) + print*,'exists = ',exists + if(.not.exists)then + print*,'The OSOCI needs a restart WF' + print*,'There are none in the EZFIO file ...' + print*,'Stopping ...' + stop + endif + print*,'passed N_det_generators_restart' + call ezfio_get_determinants_n_det(N_det_generators_restart) + ASSERT (N_det_generators_restart > 0) + ifirst = 1 +!endif +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ] + implicit none + BEGIN_DOC + ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file + ! is empty + END_DOC + integer :: i + logical :: exists + character*64 :: label + + integer, save :: ifirst = 0 +!if(ifirst == 0)then + provide N_det_generators_restart + if(.True.)then + call ezfio_has_determinants_N_int(exists) + if (exists) then + call ezfio_has_determinants_bit_kind(exists) + if (exists) then + call ezfio_has_determinants_N_det(exists) + if (exists) then + call ezfio_has_determinants_N_states(exists) + if (exists) then + call ezfio_has_determinants_psi_det(exists) + endif + endif + endif + endif + + if(.not.exists)then + print*,'The OSOCI needs a restart WF' + print*,'There are none in the EZFIO file ...' + print*,'Stopping ...' + stop + endif + print*,'passed psi_det_generators_restart' + + call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) + do i = 1, N_int + ref_generators_restart(i,1) = psi_det_generators_restart(i,1,1) + ref_generators_restart(i,2) = psi_det_generators_restart(i,2,1) + enddo + endif + ifirst = 1 +!endif + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (psi_det_size,N_states_diag) ] + implicit none + BEGIN_DOC + ! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file + ! is empty + END_DOC + + integer :: i,k, N_int2 + logical :: exists + double precision, allocatable :: psi_coef_read(:,:) + character*(64) :: label + + integer, save :: ifirst = 0 +!if(ifirst == 0)then + psi_coef_generators_restart = 0.d0 + do i=1,N_states_diag + psi_coef_generators_restart(i,i) = 1.d0 + enddo + + call ezfio_has_determinants_psi_coef(exists) + + if(.not.exists)then + print*,'The OSOCI needs a restart WF' + print*,'There are none in the EZFIO file ...' + print*,'Stopping ...' + stop + endif + print*,'passed psi_coef_generators_restart' + + if (exists) then + + allocate (psi_coef_read(N_det_generators_restart,N_states)) + call ezfio_get_determinants_psi_coef(psi_coef_read) + do k=1,N_states + do i=1,N_det_generators_restart + psi_coef_generators_restart(i,k) = psi_coef_read(i,k) + enddo + enddo + deallocate(psi_coef_read) + + endif + ifirst = 1 +!endif + + + +END_PROVIDER + diff --git a/plugins/FOBOCI/modify_generators.irp.f b/plugins/FOBOCI/modify_generators.irp.f new file mode 100644 index 00000000..c756f0c2 --- /dev/null +++ b/plugins/FOBOCI/modify_generators.irp.f @@ -0,0 +1,157 @@ +subroutine set_generators_to_psi_det + implicit none + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det_generators = N_det + integer :: i,k + do i=1,N_det_generators + do k=1,N_int + psi_det_generators(k,1,i) = psi_det(k,1,i) + psi_det_generators(k,2,i) = psi_det(k,2,i) + enddo + do k = 1, N_states + psi_coef_generators(i,k) = psi_coef(i,k) + enddo + enddo + + touch N_det_generators psi_coef_generators psi_det_generators + +end + +subroutine set_generators_as_input_psi(ndet_input,psi_det_input,psi_coef_input) + implicit none + integer, intent(in) :: ndet_input + integer(bit_kind), intent(in) :: psi_det_input(N_int,2,ndet_input) + double precision, intent(in) :: psi_coef_input(ndet_input,N_states) + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det_generators = ndet_input + integer :: i,k + do i=1,N_det_generators + do k=1,N_int + psi_det_generators(k,1,i) = psi_det_input(k,1,i) + psi_det_generators(k,2,i) = psi_det_input(k,2,i) + enddo + do k = 1, N_states + psi_coef_generators(i,k) = psi_coef_input(i,k) + enddo + enddo + + touch N_det_generators psi_coef_generators psi_det_generators + +end + +subroutine set_psi_det_as_input_psi(ndet_input,psi_det_input,psi_coef_input) + implicit none + integer, intent(in) :: ndet_input + integer(bit_kind), intent(in) :: psi_det_input(N_int,2,ndet_input) + double precision, intent(in) :: psi_coef_input(ndet_input,N_states) + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det= ndet_input + if (psi_det_size < N_det) then + psi_det_size = N_det + TOUCH psi_det_size + endif + + integer :: i,k + do i=1,N_det + do k=1,N_int + psi_det(k,1,i) = psi_det_input(k,1,i) + psi_det(k,2,i) = psi_det_input(k,2,i) + enddo + do k = 1, N_states + psi_coef(i,k) = psi_coef_input(i,k) + enddo + enddo + + soft_touch N_det psi_coef psi_det + +end + + +subroutine set_psi_det_to_generators + implicit none + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det= N_det_generators + integer :: i,k + do i = 1, psi_det_size + do k=1,N_int + psi_det(k,1,i) = 0_bit_kind + psi_det(k,2,i) = 0_bit_kind + enddo + do k = 1, N_states + psi_coef(i,k) = 0.d0 + enddo + enddo + do i=1,N_det_generators + do k=1,N_int + psi_det(k,1,i) = psi_det_generators(k,1,i) + psi_det(k,2,i) = psi_det_generators(k,2,i) + enddo + do k = 1, N_states + psi_coef(i,k) = psi_coef_generators(i,k) + enddo + enddo + + touch N_det psi_coef psi_det + +end + + + +subroutine set_generators_to_generators_restart + implicit none + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det_generators = N_det_generators_restart + integer :: i,k + do i=1,N_det_generators + do k=1,N_int + psi_det_generators(k,1,i) = psi_det_generators_restart(k,1,i) + psi_det_generators(k,2,i) = psi_det_generators_restart(k,2,i) + enddo + do k = 1, N_states + psi_coef_generators(i,k) = psi_coef_generators_restart(i,k) + enddo + enddo + + touch N_det_generators psi_coef_generators psi_det_generators + +end + + +subroutine set_psi_det_to_generators_restart + implicit none + BEGIN_DOC +! subroutines that sets psi_det_generators to +! the current psi_det + END_DOC + N_det = N_det_generators_restart + integer :: i,k + do i=1,N_det_generators + do k=1,N_int + psi_det(k,1,i) = psi_det_generators_restart(k,1,i) + psi_det(k,2,i) = psi_det_generators_restart(k,2,i) + enddo + do k = 1, N_states + psi_coef(i,k) = psi_coef_generators_restart(i,k) + enddo + enddo + + touch N_det psi_coef psi_det + +end + + diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f new file mode 100644 index 00000000..49dcafc3 --- /dev/null +++ b/plugins/FOBOCI/new_approach.irp.f @@ -0,0 +1,413 @@ + +subroutine new_approach + use bitmasks + implicit none + integer :: n_max_good_det + n_max_good_det = n_inact_orb * n_act_orb *n_det_generators_restart + n_virt_orb * n_act_orb * n_det_generators_restart + integer :: n_good_det,n_good_hole, n_good_particl + n_good_det = 0 + n_good_hole = 0 + n_good_particl = 0 + integer(bit_kind), allocatable :: psi_good_det(:,:,:) + double precision, allocatable :: dressing_restart_good_det(:,:) + double precision, allocatable :: dressing_matrix_restart_1h1p(:,:) + double precision, allocatable :: dressing_matrix_restart_2h1p(:,:) + double precision, allocatable :: dressing_matrix_restart_1h2p(:,:) + double precision, allocatable :: dressing_diag_good_det(:) + + double precision :: hjk + + integer :: i,j,k,l,i_hole_foboci + logical :: test_sym + double precision :: thr,hij + double precision :: threshold,accu + double precision, allocatable :: dressing_matrix_1h1p(:,:) + double precision, allocatable :: dressing_matrix_2h1p(:,:) + double precision, allocatable :: dressing_matrix_1h2p(:,:) + double precision, allocatable :: H_matrix_tmp(:,:) + logical :: verbose,is_ok + + double precision,allocatable :: eigenvectors(:,:), eigenvalues(:) + + + allocate(psi_good_det(N_int,2,n_max_good_det)) + allocate(dressing_restart_good_det(n_max_good_det,n_det_generators_restart)) + allocate(dressing_matrix_restart_1h1p(N_det_generators_restart, N_det_generators_restart)) + allocate(dressing_matrix_restart_2h1p(N_det_generators_restart, N_det_generators_restart)) + allocate(dressing_matrix_restart_1h2p(N_det_generators_restart, N_det_generators_restart)) + allocate(dressing_diag_good_det(n_max_good_det)) + + dressing_restart_good_det = 0.d0 + dressing_matrix_restart_1h1p = 0.d0 + dressing_matrix_restart_2h1p = 0.d0 + dressing_matrix_restart_1h2p = 0.d0 + dressing_diag_good_det = 0.d0 + + + verbose = .True. + threshold = threshold_singles + print*,'threshold = ',threshold + thr = 1.d-12 + print*,'' + print*,'' + print*,'mulliken spin population analysis' + accu =0.d0 + do i = 1, nucl_num + accu += mulliken_spin_densities(i) + print*,i,nucl_charge(i),mulliken_spin_densities(i) + enddo + print*,'' + print*,'' + print*,'DOING FIRST LMCT !!' + integer :: i_particl_osoci + + do i = 1, n_inact_orb + i_hole_foboci = list_inact(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_hole_foboci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_hole_foboci = ',i_hole_foboci + call create_restart_and_1h(i_hole_foboci) +! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + if(.not.is_ok)cycle + ! so all the mono excitation on the new generators + allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_2h1p(N_det_generators,N_det_generators)) + dressing_matrix_1h1p = 0.d0 + dressing_matrix_2h1p = 0.d0 + if(.not.do_it_perturbative)then + n_good_hole +=1 +! call all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + call all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + allocate(H_matrix_tmp(N_det_generators,N_det_generators)) + do j = 1,N_det_generators + do k = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,j),psi_det_generators(1,1,k),N_int,hjk) + H_matrix_tmp(j,k) = hjk + enddo + enddo + do j = 1, N_det_generators + do k = 1, N_det_generators + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + enddo + enddo + hjk = H_matrix_tmp(1,1) + do j = 1, N_det_generators + H_matrix_tmp(j,j) -= hjk + enddo + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'Dressed matrix :' + do j = 1, N_det_generators + write(*,'(100(X,F8.5))') H_matrix_tmp(j,:) + enddo + allocate(eigenvectors(N_det_generators,N_det_generators), eigenvalues(N_det_generators)) + call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,N_det_generators,N_det_generators) + print*,'Eigenvector of the dressed matrix :' + do j = 1, N_det_generators + print*,'coef = ',eigenvectors(j,1) + enddo + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + deallocate(eigenvectors, eigenvalues) + deallocate(H_matrix_tmp) + call update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_restart_good_det,dressing_matrix_restart_1h1p, & + dressing_matrix_restart_2h1p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det) + endif + deallocate(dressing_matrix_1h1p) + deallocate(dressing_matrix_2h1p) + enddo + + print*,'' + print*,'' + print*,'DOING THEN THE MLCT !!' + do i = 1, n_virt_orb + i_particl_osoci = list_virt(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call set_generators_to_generators_restart + call set_psi_det_to_generators + call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle + print*,'i_part_foboci = ',i_particl_osoci + call create_restart_and_1p(i_particl_osoci) + ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) + call is_a_good_candidate(threshold,is_ok,verbose) + print*,'is_ok = ',is_ok + if(.not.is_ok)cycle + ! so all the mono excitation on the new generators + allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_1h2p(N_det_generators,N_det_generators)) + dressing_matrix_1h1p = 0.d0 + dressing_matrix_1h2p = 0.d0 + if(.not.do_it_perturbative)then + n_good_hole +=1 +! call all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + call all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + allocate(H_matrix_tmp(N_det_generators,N_det_generators)) + do j = 1,N_det_generators + do k = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,j),psi_det_generators(1,1,k),N_int,hjk) + H_matrix_tmp(j,k) = hjk + enddo + enddo + do j = 1, N_det_generators + do k = 1, N_det_generators + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + enddo + enddo + hjk = H_matrix_tmp(1,1) + do j = 1, N_det_generators + H_matrix_tmp(j,j) -= hjk + enddo + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'Dressed matrix :' + do j = 1, N_det_generators + write(*,'(100(F8.5))') H_matrix_tmp(j,:) + enddo + allocate(eigenvectors(N_det_generators,N_det_generators), eigenvalues(N_det_generators)) + call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,N_det_generators,N_det_generators) + print*,'Eigenvector of the dressed matrix :' + do j = 1, N_det_generators + print*,'coef = ',eigenvectors(j,1) + enddo + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + print*,'-----------------------' + deallocate(eigenvectors, eigenvalues) + deallocate(H_matrix_tmp) + call update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_restart_good_det,dressing_matrix_restart_1h1p, & + dressing_matrix_restart_1h2p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det) + + endif + deallocate(dressing_matrix_1h1p) + deallocate(dressing_matrix_1h2p) + enddo + double precision, allocatable :: H_matrix_total(:,:) + integer :: n_det_total + n_det_total = N_det_generators_restart + n_good_det + allocate(H_matrix_total(n_det_total, n_det_total)) + ! Building of the effective Hamiltonian + ! We assume that the first determinants are the n_det_generators_restart ones + ! and then come the n_good_det determinants in psi_good_det + H_matrix_total = 0.d0 + do i = 1, N_det_generators_restart + do j = 1, N_det_generators_restart + call i_H_j(psi_det_generators_restart(1,1,i),psi_det_generators_restart(1,1,j),N_int,hij) + H_matrix_total(i,j) = hij + !!! Adding the averaged dressing coming from the 1h1p that are redundant for each of the "n_good_hole" 1h + H_matrix_total(i,j) += dressing_matrix_restart_1h1p(i,j)/dble(n_good_hole+n_good_particl) + !!! Adding the dressing coming from the 2h1p that are not redundant for the any of CI calculations + H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + enddo + enddo + do i = 1, n_good_det + call i_H_j(psi_good_det(1,1,i),psi_good_det(1,1,i),N_int,hij) + !!! Adding the diagonal dressing coming from the singles + H_matrix_total(n_det_generators_restart+i,n_det_generators_restart+i) = hij + dressing_diag_good_det(i) + do j = 1, N_det_generators_restart + !!! Adding the extra diagonal dressing between the references and the singles + print*,' dressing_restart_good_det = ',dressing_restart_good_det(i,j) + call i_H_j(psi_good_det(1,1,i),psi_det_generators_restart(1,1,j),N_int,hij) + H_matrix_total(n_det_generators_restart+i,j) += hij + H_matrix_total(j,n_det_generators_restart+i) += hij + H_matrix_total(j,n_det_generators_restart+i) += dressing_restart_good_det(i,j) + H_matrix_total(n_det_generators_restart+i,j) += dressing_restart_good_det(i,j) + enddo + do j = i+1, n_good_det + !!! Building the naked Hamiltonian matrix between the singles + call i_H_j(psi_good_det(1,1,i),psi_good_det(1,1,j),N_int,hij) + H_matrix_total(n_det_generators_restart+i,n_det_generators_restart+j) = hij + H_matrix_total(n_det_generators_restart+j,n_det_generators_restart+i) = hij + enddo + enddo + print*,'H matrix to diagonalize' + double precision :: href + href = H_matrix_total(1,1) + do i = 1, n_det_total + H_matrix_total(i,i) -= href + enddo + do i = 1, n_det_total + write(*,'(100(X,F16.8))')H_matrix_total(i,:) + enddo + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total)) + call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total) + print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion + href + do i = 1, n_det_total + print*,'coef = ',eigvectors(i,1) + enddo + integer(bit_kind), allocatable :: psi_det_final(:,:,:) + double precision, allocatable :: psi_coef_final(:,:) + double precision :: norm + allocate(psi_coef_final(n_det_total, N_states)) + allocate(psi_det_final(N_int,2,n_det_total)) + do i = 1, N_det_generators_restart + do j = 1,N_int + psi_det_final(j,1,i) = psi_det_generators_restart(j,1,i) + psi_det_final(j,2,i) = psi_det_generators_restart(j,2,i) + enddo + enddo + do i = 1, n_good_det + do j = 1,N_int + psi_det_final(j,1,n_det_generators_restart+i) = psi_good_det(j,1,i) + psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i) + enddo + enddo + norm = 0.d0 + do i = 1, n_det_total + do j = 1, N_states + psi_coef_final(i,j) = eigvectors(i,j) + enddo + norm += psi_coef_final(i,1)**2 +! call debug_det(psi_det_final(1, 1, i), N_int) + enddo + print*,'norm = ',norm + + call set_psi_det_as_input_psi(n_det_total,psi_det_final,psi_coef_final) + print*,'' +!do i = 1, N_det +! call debug_det(psi_det(1,1,i),N_int) +! print*,'coef = ',psi_coef(i,1) +!enddo + provide one_body_dm_mo + + integer :: i_core,iorb,jorb,i_inact,j_inact,i_virt,j_virt,j_core + do i = 1, n_core_orb + i_core = list_core(i) + one_body_dm_mo(i_core,i_core) = 10.d0 + do j = i+1, n_core_orb + j_core = list_core(j) + one_body_dm_mo(i_core,j_core) = 0.d0 + one_body_dm_mo(j_core,i_core) = 0.d0 + enddo + do j = 1, n_inact_orb + iorb = list_inact(j) + one_body_dm_mo(i_core,iorb) = 0.d0 + one_body_dm_mo(iorb,i_core) = 0.d0 + enddo + do j = 1, n_act_orb + iorb = list_act(j) + one_body_dm_mo(i_core,iorb) = 0.d0 + one_body_dm_mo(iorb,i_core) = 0.d0 + enddo + do j = 1, n_virt_orb + iorb = list_virt(j) + one_body_dm_mo(i_core,iorb) = 0.d0 + one_body_dm_mo(iorb,i_core) = 0.d0 + enddo + enddo + ! Set to Zero the inact-inact part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = i+1, n_inact_orb + j_inact = list_inact(j) + one_body_dm_mo(i_inact,j_inact) = 0.d0 + one_body_dm_mo(j_inact,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the inact-virt part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = 1, n_virt_orb + j_virt = list_virt(j) + one_body_dm_mo(i_inact,j_virt) = 0.d0 + one_body_dm_mo(j_virt,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the virt-virt part to avoid arbitrary rotations + do i = 1, n_virt_orb + i_virt = list_virt(i) + do j = i+1, n_virt_orb + j_virt = list_virt(j) + one_body_dm_mo(i_virt,j_virt) = 0.d0 + one_body_dm_mo(j_virt,i_virt) = 0.d0 + enddo + enddo + + + print*,'' + print*,'Inactive-active Part of the One body DM' + print*,'' + do i = 1,n_act_orb + iorb = list_act(i) + print*,'' + print*,'ACTIVE ORBITAL ',iorb + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + print*,'INACTIVE ' + print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) + endif + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + print*,'VIRT ' + print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) + endif + enddo + enddo + do i = 1, mo_tot_num + do j = i+1, mo_tot_num + if(dabs(one_body_dm_mo(i,j)).le.threshold_fobo_dm)then + one_body_dm_mo(i,j) = 0.d0 + one_body_dm_mo(j,i) = 0.d0 + endif + enddo + enddo + + + + + + + label = "Natural" + character*(64) :: label + integer :: sign + sign = -1 + + call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign) + soft_touch mo_coef + call save_mos + + deallocate(eigvalues,eigvectors,psi_det_final,psi_coef_final) + + + + + deallocate(H_matrix_total) + deallocate(psi_good_det) + deallocate(dressing_restart_good_det) + deallocate(dressing_matrix_restart_1h1p) + deallocate(dressing_matrix_restart_2h1p) + deallocate(dressing_diag_good_det) + +end + + diff --git a/plugins/FOBOCI/routine_new_approach.irp.f b/plugins/FOBOCI/routine_new_approach.irp.f new file mode 100644 index 00000000..34274b76 --- /dev/null +++ b/plugins/FOBOCI/routine_new_approach.irp.f @@ -0,0 +1,56 @@ +subroutine update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_restart_good_det,dressing_matrix_restart_1h1p, & + dressing_matrix_restart_2h1p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det) + implicit none + integer, intent(in) :: n_max_good_det + integer, intent(inout) :: n_good_det + integer(bit_kind), intent(inout) :: psi_good_det(N_int,2,n_max_good_det) + double precision, intent(in) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) + double precision, intent(in) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) + double precision, intent(inout) :: dressing_matrix_restart_1h1p(N_det_generators_restart, N_det_generators_restart) + double precision, intent(inout) :: dressing_matrix_restart_2h1p(N_det_generators_restart, N_det_generators_restart) + double precision, intent(inout) :: dressing_restart_good_det(n_max_good_det,N_det_generators_restart) + double precision, intent(inout) :: dressing_diag_good_det(n_max_good_det) + use bitmasks + integer :: k,l,degree + logical, allocatable :: is_a_ref_det(:) + integer, allocatable :: index_restart_generators(:) + allocate(is_a_ref_det(N_det_generators),index_restart_generators(N_det_generators)) + is_a_ref_det = .False. + do k = 1, N_det_generators + do l = 1, N_det_generators_restart + call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators_restart(1,1,l), degree, N_int) + if(degree==0)then + is_a_ref_det(k) = .True. + index_restart_generators(k) = l + endif + enddo + enddo + do k = 1, N_det_generators + if(is_a_ref_det(k))then + do l = 1, N_det_generators + if(.not.is_a_ref_det(l))cycle + !!!! Dressing of the reference space in the order of the restart determinants + dressing_matrix_restart_1h1p(index_restart_generators(l),index_restart_generators(k)) += dressing_matrix_1h1p(k,l) + print*,' dressing_matrix_1h1p(k,l) = ',dressing_matrix_1h1p(k,l) + dressing_matrix_restart_2h1p(index_restart_generators(l),index_restart_generators(k)) += dressing_matrix_2h1p(k,l) + enddo + else + !!!! Incrementing the counting of the good determinants + n_good_det +=1 + !!!! Adding the good determinant to the global_list (psi_good_det) + do l = 1, N_int + psi_good_det(l,1,n_good_det) = psi_det_generators(l,1,k) + psi_good_det(l,2,n_good_det) = psi_det_generators(l,2,k) + enddo + !!! Storing the diagonal dressing of the good det + dressing_diag_good_det(n_good_det) = dressing_matrix_1h1p(k,k) + dressing_matrix_2h1p(k,k) + do l = 1, N_det_generators + if(.not.is_a_ref_det(l))cycle + !!! Storing the extra diagonal dressing of the good det with the restart determinants + dressing_restart_good_det(n_good_det,index_restart_generators(l)) = dressing_matrix_1h1p(k,l) + dressing_matrix_2h1p(k,l) + enddo + endif + enddo + deallocate(is_a_ref_det,index_restart_generators) + +end diff --git a/plugins/FOBOCI/routines_dressing.irp.f b/plugins/FOBOCI/routines_dressing.irp.f new file mode 100644 index 00000000..910f1109 --- /dev/null +++ b/plugins/FOBOCI/routines_dressing.irp.f @@ -0,0 +1,457 @@ +subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) + use bitmasks + implicit none + integer, intent(in) :: ndet_generators_input + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input) + double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input) + double precision :: H_array(N_det),hka + logical :: is_a_ref_det(N_det) + integer :: i,j,n_det_ref_tmp + integer :: connected_to_ref_by_mono,degree + double precision :: coef_ref(Ndet_generators_input) + double precision :: accu,lambda_i + integer :: k + integer :: index_ref_tmp(N_det) + is_a_ref_det = .False. + n_det_ref_tmp = 0 + do i = 1, N_det + do j = 1, Ndet_generators_input + call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int) + if(degree == 0)then + is_a_ref_det(i) = .True. + n_det_ref_tmp +=1 + index_ref_tmp(n_det_ref_tmp) = i + coef_ref(n_det_ref_tmp) = psi_coef(i,1) + exit + endif + enddo + enddo + if( ndet_generators_input .ne. n_det_ref_tmp)then + print*,'Problem !!!! ' + print*,' ndet_generators .ne. n_det_ref_tmp !!!' + print*,'ndet_generators,n_det_ref_tmp' + print*,ndet_generators_input,n_det_ref_tmp + stop + endif + + call i_h_j(psi_det_generators_input(1,1,1),psi_det_generators_input(1,1,1),N_int,href) + integer :: i_pert, i_pert_count + i_pert_count = 0 + do i = 1, N_det + if(is_a_ref_det(i))cycle + call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hka) + double precision :: f,href + f = 1.d0/(href - hka) + H_array = 0.d0 + accu = 0.d0 + do j=1,ndet_generators_input + call i_h_j(psi_det(1,1,i),psi_det_generators_input(1,1,j),N_int,hka) + H_array(j) = hka + accu += coef_ref(j) * hka + enddo + lambda_i = psi_coef(i,1)/accu + i_pert = 1 + if(accu * f / psi_coef(i,1) .gt. 0.5d0 .and. accu * f/psi_coef(i,1).gt.0.d0)then + i_pert = 0 + endif + do j = 1, ndet_generators_input + if(dabs(H_array(j)*lambda_i).gt.0.5d0)then + i_pert = 1 + exit + endif + enddo +! print*,'' +! print*,'lambda_i,f = ',lambda_i,f +! print*,'i_pert = ',i_pert +! print*,'' + if(i_pert==1)then + lambda_i = f + i_pert_count +=1 + endif + do k=1,ndet_generators_input + double precision :: contrib + contrib = H_array(k) * H_array(k) * lambda_i + dressing_matrix(k, k) += contrib + do j=k+1,ndet_generators_input + contrib = H_array(k) * H_array(j) * lambda_i + dressing_matrix(k, j) += contrib + dressing_matrix(j, k) += contrib + enddo + enddo + enddo +!print*,'i_pert_count = ',i_pert_count +end + +subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, & + psi_det_outer_input,psi_coef_outer_input,n_det_outer_input) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_input + integer(bit_kind), intent(in) :: psi_det_ref_input(N_int,2,n_det_ref_input) + double precision, intent(in) :: psi_coef_ref_input(n_det_ref_input,N_states) + integer, intent(in) :: n_det_outer_input + integer(bit_kind), intent(in) :: psi_det_outer_input(N_int,2,n_det_outer_input) + double precision, intent(in) :: psi_coef_outer_input(n_det_outer_input,N_states) + + double precision, intent(inout) :: dressing_matrix(n_det_ref_input,n_det_ref_input) + + + integer :: i_pert, i_pert_count,i,j,k + double precision :: f,href,hka,lambda_i + double precision :: H_array(n_det_ref_input),accu + call i_h_j(psi_det_ref_input(1,1,1),psi_det_ref_input(1,1,1),N_int,href) + i_pert_count = 0 + do i = 1, n_det_outer_input + call i_h_j(psi_det_outer_input(1,1,i),psi_det_outer_input(1,1,i),N_int,hka) + f = 1.d0/(href - hka) + H_array = 0.d0 + accu = 0.d0 + do j=1,n_det_ref_input + call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka) + H_array(j) = hka + accu += psi_coef_ref_input(j,1) * hka + enddo + lambda_i = psi_coef_outer_input(i,1)/accu + i_pert = 1 + if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then + i_pert = 0 + endif + do j = 1, n_det_ref_input + if(dabs(H_array(j)*lambda_i).gt.0.3d0)then + i_pert = 1 + exit + endif + enddo + if(i_pert==1)then + lambda_i = f + i_pert_count +=1 + endif + do k=1,n_det_ref_input + double precision :: contrib + contrib = H_array(k) * H_array(k) * lambda_i + dressing_matrix(k, k) += contrib + do j=k+1,n_det_ref_input + contrib = H_array(k) * H_array(j) * lambda_i + dressing_matrix(k, j) += contrib + dressing_matrix(j, k) += contrib + enddo + enddo + enddo +end + + +subroutine diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_generators_input,dressing_matrix) + use bitmasks + implicit none + integer, intent(in) :: ndet_generators_input + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input) + double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input) + integer :: i,j + + double precision :: eigenvectors(Ndet_generators_input,Ndet_generators_input), eigenvalues(Ndet_generators_input) + + call lapack_diag(eigenvalues,eigenvectors,dressing_matrix,Ndet_generators_input,Ndet_generators_input) + print*,'Dressed eigenvalue, to be compared with the CI one' + print*,'E = ',eigenvalues(1)+nuclear_repulsion + print*,'Dressed matrix, to be compared to the intermediate Hamiltonian one' + do i = 1, Ndet_generators_input + write(*,'(100(F12.5,X))')dressing_matrix(i,:) + enddo + n_det = Ndet_generators_input + do i = 1, Ndet_generators_input + psi_coef(i,1) = eigenvectors(i,1) + do j = 1, N_int + psi_det(j,1,i) = psi_det_generators_input(j,1,i) + psi_det(j,2,i) = psi_det_generators_input(j,2,i) + enddo + enddo + + touch N_det psi_coef psi_det + +end + +subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) + use bitmasks + implicit none + integer, intent(out) :: n_det_1h1p, n_det_2h1p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + n_det_ref_restart_tmp = 0 + n_det_1h = 0 + n_det_1h1p = 0 + n_det_2h1p = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_restart_tmp +=1 + else if (n_h ==1 .and. n_p==0)then + n_det_1h +=1 + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p +=1 + else if (n_h ==2 .and. n_p==1)then + n_det_2h1p +=1 + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1h1p or 2h1p' + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo +! if(n_det_1h.ne.1)then +! print*,'PB !! You have more than one 1h' +! stop +! endif + if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then + print*,'PB !!!!' + print*,'You have forgotten something in your generators ... ' + stop + endif + + +end + +subroutine give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) + use bitmasks + implicit none + integer, intent(out) :: n_det_1h1p, n_det_1h2p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + n_det_ref_restart_tmp = 0 + n_det_1h = 0 + n_det_1h1p = 0 + n_det_1h2p = 0 + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + n_det_ref_restart_tmp +=1 + else if (n_h ==0 .and. n_p==1)then + n_det_1h +=1 + else if (n_h ==1 .and. n_p==1)then + n_det_1h1p +=1 + else if (n_h ==1 .and. n_p==2)then + n_det_1h2p +=1 + else + print*,'PB !!!!' + print*,'You have something else than a 1p, 1h1p or 1h2p' + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then + print*,'PB !!!!' + print*,'You have forgotten something in your generators ... ' + stop + endif + + +end + + +subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) + use bitmasks + implicit none + integer, intent(in) :: n_det_1h1p,n_det_2h1p + integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) + integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) + integer(bit_kind), intent(out) :: psi_2h1p(N_int,2,n_det_2h1p) + double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) + double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p, N_states) + + integer :: i,j + integer :: degree + integer :: number_of_holes,n_h, number_of_particles,n_p + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp + integer, allocatable :: index_generator(:) + integer, allocatable :: index_1h1p(:) + integer, allocatable :: index_2h1p(:) + + allocate(index_1h1p(n_det)) + allocate(index_2h1p(n_det)) + allocate(index_generator(N_det)) + + + n_det_generators_tmp = 0 + n_det_1h1p_tmp = 0 + n_det_2h1p_tmp = 0 + do i = 1, n_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if (n_h ==1 .and. n_p==1)then + n_det_1h1p_tmp +=1 + index_1h1p(n_det_1h1p_tmp) = i + else if (n_h ==2 .and. n_p==1)then + n_det_2h1p_tmp +=1 + index_2h1p(n_det_2h1p_tmp) = i + endif + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) + if(degree == 0)then + n_det_generators_tmp +=1 + index_generator(n_det_generators_tmp) = i + endif + enddo + enddo + if(n_det_1h1p_tmp.ne.n_det_1h1p)then + print*,'PB !!!' + print*,'n_det_1h1p_tmp.ne.n_det_1h1p)' + stop + endif + + + if(n_det_2h1p_tmp.ne.n_det_2h1p)then + print*,'PB !!!' + print*,'n_det_2h1p_tmp.ne.n_det_2h1p)' + stop + endif + + if(N_det_generators_tmp.ne.n_det_generators)then + print*,'PB !!!' + print*,'N_det_generators_tmp.ne.n_det_generators' + stop + endif + + do i = 1,N_det_generators + do j = 1, N_int + psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i)) + psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i)) + enddo + do j = 1, N_states + psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j) + enddo + enddo + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i)) + enddo + do j = 1, N_states + psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j) + enddo + enddo + + do i = 1, n_det_2h1p + do j = 1, N_int + psi_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i)) + psi_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i)) + enddo + do j = 1, N_states + psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j) + enddo + enddo + + + deallocate(index_generator) + deallocate(index_1h1p) + deallocate(index_2h1p) + +end + + +subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) + use bitmasks + implicit none + integer, intent(in) :: n_det_1h1p,n_det_1h2p + integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) + integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) + integer(bit_kind), intent(out) :: psi_1h2p(N_int,2,n_det_1h2p) + double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) + double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p, N_states) + + integer :: i,j + integer :: degree + integer :: number_of_holes,n_h, number_of_particles,n_p + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp + integer, allocatable :: index_generator(:) + integer, allocatable :: index_1h1p(:) + integer, allocatable :: index_1h2p(:) + + allocate(index_1h1p(n_det)) + allocate(index_1h2p(n_det)) + allocate(index_generator(N_det)) + + + n_det_generators_tmp = 0 + n_det_1h1p_tmp = 0 + n_det_1h2p_tmp = 0 + do i = 1, n_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if (n_h ==1 .and. n_p==1)then + n_det_1h1p_tmp +=1 + index_1h1p(n_det_1h1p_tmp) = i + else if (n_h ==1 .and. n_p==2)then + n_det_1h2p_tmp +=1 + index_1h2p(n_det_1h2p_tmp) = i + endif + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) + if(degree == 0)then + n_det_generators_tmp +=1 + index_generator(n_det_generators_tmp) = i + endif + enddo + enddo + if(n_det_1h1p_tmp.ne.n_det_1h1p)then + print*,'PB !!!' + print*,'n_det_1h1p_tmp.ne.n_det_1h1p)' + stop + endif + + + if(n_det_1h2p_tmp.ne.n_det_1h2p)then + print*,'PB !!!' + print*,'n_det_1h2p_tmp.ne.n_det_1h2p)' + stop + endif + + if(N_det_generators_tmp.ne.n_det_generators)then + print*,'PB !!!' + print*,'N_det_generators_tmp.ne.n_det_generators' + stop + endif + + do i = 1,N_det_generators + do j = 1, N_int + psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i)) + psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i)) + enddo + do j = 1, N_states + psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j) + enddo + enddo + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i)) + enddo + do j = 1, N_states + psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j) + enddo + enddo + + do i = 1, n_det_1h2p + do j = 1, N_int + psi_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i)) + psi_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i)) + enddo + do j = 1, N_states + psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j) + enddo + enddo + + + deallocate(index_generator) + deallocate(index_1h1p) + deallocate(index_1h2p) + +end + + diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f new file mode 100644 index 00000000..696011a9 --- /dev/null +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -0,0 +1,616 @@ +subroutine set_intermediate_normalization_lmct_old(norm,i_hole) + implicit none + integer, intent(in) :: i_hole + double precision, intent(out) :: norm(N_states) + integer :: i,j,degree,index_ref_generators_restart,k + integer:: number_of_holes,n_h, number_of_particles,n_p + integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:) + integer, allocatable :: index_one_p(:) + integer :: n_one_hole,n_one_hole_one_p,n_two_hole_one_p,n_two_hole,n_one_p + logical :: is_the_hole_in_det + double precision :: inv_coef_ref_generators_restart(N_states),hij,hii,accu + integer :: index_good_hole(1000) + integer :: n_good_hole + logical,allocatable :: is_a_ref_det(:) + allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det)) + + n_one_hole = 0 + n_one_hole_one_p = 0 + n_two_hole_one_p = 0 + n_two_hole = 0 + n_one_p = 0 + n_good_hole = 0 + ! Find the one holes and one hole one particle + is_a_ref_det = .False. + do i = 1, N_det + ! Find the reference determinant for intermediate normalization + call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart = i + do k = 1, N_states + inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k) + enddo +! cycle + endif + + ! Find all the determinants present in the reference wave function + do j = 1, N_det_generators_restart + call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int) + if(degree == 0)then + is_a_ref_det(i) = .True. + exit + endif + enddo + if(is_a_ref_det(i))cycle + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 1 .and. n_p == 0)then + if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then + n_good_hole +=1 + index_good_hole(n_good_hole) = i + else + do k = 1, N_states + psi_coef(i,k) = 0.d0 + enddo + endif + else + do k = 1, N_states + psi_coef(i,k) = 0.d0 + enddo + endif + enddo +!do k = 1, N_det +! call debug_det(psi_det(1,1,k),N_int) +! print*,'k,coef = ',k,psi_coef(k,1)/psi_coef(index_ref_generators_restart,1) +!enddo + print*,'' + print*,'n_good_hole = ',n_good_hole + do k = 1,N_states + print*,'state ',k + do i = 1, n_good_hole + print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart,k) + enddo + print*,'' + enddo + norm = 0.d0 + + ! Set the wave function to the intermediate normalization + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k) + enddo + enddo + do k = 1,N_states + print*,'state ',k + do i = 1, N_det +!! print*,'psi_coef(i_ref) = ',psi_coef(i,1) + if (is_a_ref_det(i))then + print*,'i,psi_coef_ref = ',psi_coef(i,k) + cycle + endif + norm(k) += psi_coef(i,k) * psi_coef(i,k) + enddo + print*,'norm = ',norm(k) + enddo + deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det) + soft_touch psi_coef +end + + +subroutine set_intermediate_normalization_mlct_old(norm,i_particl) + implicit none + integer, intent(in) :: i_particl + double precision, intent(out) :: norm(N_states) + integer :: i,j,degree,index_ref_generators_restart,k + integer:: number_of_holes,n_h, number_of_particles,n_p + integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:) + integer, allocatable :: index_one_p(:),index_one_hole_two_p(:) + integer :: n_one_hole,n_one_hole_one_p,n_two_hole_one_p,n_two_hole,n_one_p,n_one_hole_two_p + logical :: is_the_particl_in_det + double precision :: inv_coef_ref_generators_restart(N_states) + integer :: exc(0:2,2,2) + double precision :: phase,hij,hii,accu + integer :: h1,p1,h2,p2,s1,s2 + integer :: index_good_particl(1000) + integer :: n_good_particl + logical,allocatable :: is_a_ref_det(:) + integer :: i_count + allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det)) + allocate(index_one_hole_two_p(n_det)) + + n_one_hole = 0 + n_one_hole_one_p = 0 + n_two_hole_one_p = 0 + n_two_hole = 0 + n_one_p = 0 + n_one_hole_two_p = 0 + n_good_particl = 0 + ! Find the one holes and one hole one particle + i_count = 0 + is_a_ref_det = .False. + do i = 1, N_det + call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int) + if(degree == 0)then + index_ref_generators_restart = i + do k = 1, N_states + inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k) + enddo +! cycle + endif + + ! Find all the determinants present in the reference wave function + do j = 1, N_det_generators_restart + call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int) + if(degree == 0)then + is_a_ref_det(i) = .True. + exit + endif + enddo + if(is_a_ref_det(i))cycle + + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 1)then ! 1p + if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then + n_good_particl += 1 + index_good_particl(n_good_particl) = i + else + do k = 1, N_states + psi_coef(i,k) = 0.d0 + enddo + endif + else + do k = 1, N_states + psi_coef(i,k) = 0.d0 + enddo + endif + enddo + + norm = 0.d0 + print*,'' + print*,'n_good_particl = ',n_good_particl + do k = 1, N_states + print*,'state ',k + do i = 1, n_good_particl + print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart,k) + enddo + print*,'' + enddo + + + ! Set the wave function to the intermediate normalization + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k) + enddo + enddo + do k = 1, N_states + print*,'state ',k + do i = 1, N_det +!! print*,'i = ',i, psi_coef(i,1) + if (is_a_ref_det(i))then + print*,'i,psi_coef_ref = ',psi_coef(i,k) + cycle + endif + norm(k) += psi_coef(i,k) * psi_coef(i,k) + enddo + print*,'norm = ',norm + enddo + soft_touch psi_coef + deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det) +end + + +subroutine update_density_matrix_osoci + implicit none + BEGIN_DOC + ! one_body_dm_mo_alpha_osoci += Delta rho alpha + ! one_body_dm_mo_beta_osoci += Delta rho beta + END_DOC + integer :: i,j + integer :: iorb,jorb + do i = 1, mo_tot_num + do j = 1, mo_tot_num + one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha(i,j) - one_body_dm_mo_alpha_generators_restart(i,j)) + one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta(i,j) - one_body_dm_mo_beta_generators_restart(i,j)) + enddo + enddo + + +end + + +subroutine initialize_density_matrix_osoci + implicit none + one_body_dm_mo_alpha_osoci = one_body_dm_mo_alpha_generators_restart + one_body_dm_mo_beta_osoci = one_body_dm_mo_beta_generators_restart +end + +subroutine rescale_density_matrix_osoci(norm) + implicit none + double precision, intent(in) :: norm(N_states) + integer :: i,j + double precision :: norm_tmp + norm_tmp = 0.d0 + do i = 1, N_states + norm_tmp += norm(i) + enddo + print*,'norm = ',norm_tmp + + do i = 1, mo_tot_num + do j = 1,mo_tot_num + one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) * norm_tmp + one_body_dm_mo_beta_osoci(j,i) = one_body_dm_mo_beta_osoci(j,i) * norm_tmp + enddo + enddo +end + +subroutine save_osoci_natural_mos + + implicit none + BEGIN_DOC + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:),tmp_bis(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha,i,i_core,j_core,iorb,jorb,j,i_inact,j_inact,i_virt,j_virt + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + allocate(tmp_bis(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + allocate (occ(N_int*bit_kind_size,2)) + + ! Negation to have the occupied MOs first after the diagonalization + tmp_bis = -one_body_dm_mo_alpha_osoci - one_body_dm_mo_beta_osoci + ! Set to Zero the core-inact-act-virt part + do i = 1, n_core_orb + i_core = list_core(i) + tmp_bis(i_core,i_core) = -10.d0 + do j = i+1, n_core_orb + j_core = list_core(j) + tmp_bis(i_core,j_core) = 0.d0 + tmp_bis(j_core,i_core) = 0.d0 + enddo + do j = 1, n_inact_orb + iorb = list_inact(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + do j = 1, n_act_orb + iorb = list_act(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + do j = 1, n_virt_orb + iorb = list_virt(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + enddo + do i = 1, n_core_orb + print*,'dm core = ',list_core(i),tmp_bis(list_core(i),list_core(i)) + enddo + ! Set to Zero the inact-inact part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = i+1, n_inact_orb + j_inact = list_inact(j) + tmp_bis(i_inact,j_inact) = 0.d0 + tmp_bis(j_inact,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the inact-virt part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = 1, n_virt_orb + j_virt = list_virt(j) + tmp_bis(i_inact,j_virt) = 0.d0 + tmp_bis(j_virt,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the virt-virt part to avoid arbitrary rotations + do i = 1, n_virt_orb + i_virt = list_virt(i) + do j = i+1, n_virt_orb + j_virt = list_virt(j) + tmp_bis(i_virt,j_virt) = 0.d0 + tmp_bis(j_virt,i_virt) = 0.d0 + enddo + enddo + + double precision :: accu + ! Set to Zero the act-act part to avoid arbitrary rotations + do i = 1,n_act_orb + iorb = list_act(i) + do j = i+1,n_act_orb + jorb = list_act(j) + tmp_bis(iorb,jorb) = 0.d0 + tmp_bis(jorb,iorb) = 0.d0 + enddo + enddo + + tmp = tmp_bis +!! Symetrization act-virt + do j = 1, n_virt_orb + j_virt= list_virt(j) + accu = 0.d0 + do i = 1, n_act_orb + jorb = list_act(i) + accu += dabs(tmp_bis(j_virt,jorb)) + enddo + do i = 1, n_act_orb + iorb = list_act(i) + tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) + tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb)) + enddo + enddo + +!! Symetrization act-inact +!do j = 1, n_inact_orb +! j_inact = list_inact(j) +! accu = 0.d0 +! do i = 1, n_act_orb +! jorb = list_act(i) +! accu += dabs(tmp_bis(j_inact,jorb)) +! enddo +! do i = 1, n_act_orb +! iorb = list_act(i) +! tmp(j_inact,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_inact,iorb)) +! tmp(iorb,j_inact) = dsign(accu/dble(n_act_orb),tmp_bis(j_inact,iorb)) +! enddo +!enddo + +!!! Symetrization act-act +!!accu = 0.d0 +!!do i = 1, n_act_orb +!! iorb = list_act(i) +!! accu += tmp_bis(iorb,iorb) +!!enddo +!!do i = 1, n_act_orb +!! iorb = list_act(i) +!! tmp(iorb,iorb) = accu/dble(n_act_orb) +!!enddo + + call bitstring_to_list(reunion_of_bitmask(1,1), occ(1,1), n_occ_alpha, N_int) + double precision :: maxvaldm,imax,jmax + maxvaldm = 0.d0 + imax = 1 + jmax = 1 + print*,'' + print*,'Inactive-active Part of the One body DM' + print*,'' + do i = 1,n_act_orb + iorb = list_act(i) + print*,'' + print*,'ACTIVE ORBITAL ',iorb + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + print*,'INACTIVE ' + print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + endif + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + print*,'VIRT ' + print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + endif + enddo + enddo + do i = 1, mo_tot_num + do j = i+1, mo_tot_num + if(dabs(tmp(i,j)).le.threshold_fobo_dm)then + tmp(i,j) = 0.d0 + tmp(j,i) = 0.d0 + endif + enddo + enddo + + label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + soft_touch mo_coef + deallocate(tmp,occ) + + +end + +subroutine set_osoci_natural_mos + + implicit none + BEGIN_DOC + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:),tmp_bis(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha,i,i_core,j_core,iorb,jorb,j,i_inact,j_inact,i_virt,j_virt + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + allocate(tmp_bis(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + allocate (occ(N_int*bit_kind_size,2)) + + ! Negation to have the occupied MOs first after the diagonalization + tmp_bis = -one_body_dm_mo_alpha_osoci - one_body_dm_mo_beta_osoci + ! Set to Zero the core-inact-act-virt part + do i = 1, n_core_orb + i_core = list_core(i) + tmp_bis(i_core,i_core) = -10.d0 + do j = i+1, n_core_orb + j_core = list_core(j) + tmp_bis(i_core,j_core) = 0.d0 + tmp_bis(j_core,i_core) = 0.d0 + enddo + do j = 1, n_inact_orb + iorb = list_inact(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + do j = 1, n_act_orb + iorb = list_act(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + do j = 1, n_virt_orb + iorb = list_virt(j) + tmp_bis(i_core,iorb) = 0.d0 + tmp_bis(iorb,i_core) = 0.d0 + enddo + enddo + do i = 1, n_core_orb + print*,'dm core = ',list_core(i),tmp_bis(list_core(i),list_core(i)) + enddo + ! Set to Zero the inact-inact part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = i+1, n_inact_orb + j_inact = list_inact(j) + tmp_bis(i_inact,j_inact) = 0.d0 + tmp_bis(j_inact,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the inact-virt part to avoid arbitrary rotations + do i = 1, n_inact_orb + i_inact = list_inact(i) + do j = 1, n_virt_orb + j_virt = list_virt(j) + tmp_bis(i_inact,j_virt) = 0.d0 + tmp_bis(j_virt,i_inact) = 0.d0 + enddo + enddo + + ! Set to Zero the virt-virt part to avoid arbitrary rotations + do i = 1, n_virt_orb + i_virt = list_virt(i) + do j = i+1, n_virt_orb + j_virt = list_virt(j) + tmp_bis(i_virt,j_virt) = 0.d0 + tmp_bis(j_virt,i_virt) = 0.d0 + enddo + enddo + + double precision :: accu + ! Set to Zero the act-act part to avoid arbitrary rotations + do i = 1,n_act_orb + iorb = list_act(i) + do j = i+1,n_act_orb + jorb = list_act(j) + tmp_bis(iorb,jorb) = 0.d0 + tmp_bis(jorb,iorb) = 0.d0 + enddo + enddo + + tmp = tmp_bis + + call bitstring_to_list(reunion_of_bitmask(1,1), occ(1,1), n_occ_alpha, N_int) + double precision :: maxvaldm,imax,jmax + maxvaldm = 0.d0 + imax = 1 + jmax = 1 + print*,'' + print*,'Inactive-active Part of the One body DM' + print*,'' + do i = 1,n_act_orb + iorb = list_act(i) + print*,'' + print*,'ACTIVE ORBITAL ',iorb + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + print*,'INACTIVE ' + print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + endif + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + print*,'VIRT ' + print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + endif + enddo + enddo + do i = 1, mo_tot_num + do j = i+1, mo_tot_num + if(dabs(tmp(i,j)).le.threshold_fobo_dm)then + tmp(i,j) = 0.d0 + tmp(j,i) = 0.d0 + endif + enddo + enddo + + label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + soft_touch mo_coef + deallocate(tmp,occ) + + +end + +subroutine check_symetry(i_hole,thr,test) + implicit none + integer, intent(in) :: i_hole + double precision, intent(in) :: thr + logical, intent(out) :: test + integer :: i,j,k,l + double precision :: accu + accu = 0.d0 + do i = 1, n_act_orb + accu += dabs(mo_mono_elec_integral(i_hole,list_act(i))) + enddo + if(accu.gt.thr)then + test = .True. + else + test = .false. + endif +end + +subroutine check_symetry_1h1p(i_hole,i_part,thr,test) + implicit none + integer, intent(in) :: i_hole,i_part + double precision, intent(in) :: thr + logical, intent(out) :: test + integer :: i,j,k,l + double precision :: accu + accu = dabs(mo_mono_elec_integral(i_hole,i_part)) + if(accu.gt.thr)then + test = .True. + else + test = .false. + endif +end + + + subroutine update_one_body_dm_mo + implicit none + integer :: i + double precision :: accu_tot,accu_sd + print*,'touched the one_body_dm_mo_beta' + one_body_dm_mo_alpha = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta = one_body_dm_mo_beta_osoci + touch one_body_dm_mo_alpha one_body_dm_mo_beta + accu_tot = 0.d0 + accu_sd = 0.d0 + do i = 1, mo_tot_num + accu_tot += one_body_dm_mo_alpha(i,i) + one_body_dm_mo_beta(i,i) + accu_sd += one_body_dm_mo_alpha(i,i) - one_body_dm_mo_beta(i,i) + enddo + print*,'accu_tot = ',accu_tot + print*,'accu_sdt = ',accu_sd + end + + subroutine provide_properties + implicit none + integer :: i + double precision :: accu + if(.True.)then + accu= 0.d0 + do i = 1, nucl_num + accu += mulliken_spin_densities(i) + print*,i,nucl_charge(i),mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu + endif + end + diff --git a/plugins/FOBOCI/save_fock_diag_inactiv_virt.irp.f b/plugins/FOBOCI/save_fock_diag_inactiv_virt.irp.f new file mode 100644 index 00000000..dceb8546 --- /dev/null +++ b/plugins/FOBOCI/save_fock_diag_inactiv_virt.irp.f @@ -0,0 +1,7 @@ +program save_fock_inactiv_virt_mos + implicit none + call diag_inactive_virt_and_update_mos + call save_mos + + +end diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f index 56bf96d9..33966743 100644 --- a/plugins/Full_CI/full_ci.irp.f +++ b/plugins/Full_CI/full_ci.irp.f @@ -64,9 +64,9 @@ program full_ci print *, 'N_states = ', N_states do k = 1, N_states print*,'State ',k - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E(before)+PT2 = ', E_CI_before+pt2 + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) enddo print *, '-----' E_CI_before = CI_energy diff --git a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES index 5b16423e..7392852a 100644 --- a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES +++ b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full Psiref_Utils +Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS diff --git a/plugins/MRCC_Utils/mrcc_dummy.irp.f b/plugins/MRCC_Utils/mrcc_dummy.irp.f new file mode 100644 index 00000000..8f1deda8 --- /dev/null +++ b/plugins/MRCC_Utils/mrcc_dummy.irp.f @@ -0,0 +1,4 @@ +program pouet + + +end diff --git a/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES b/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..01b1f4b9 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Dressed_Ref_Hamiltonian OVB diff --git a/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f b/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f new file mode 100644 index 00000000..4ef79aeb --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f @@ -0,0 +1,59 @@ + BEGIN_PROVIDER [double precision, H_OVB_dressing, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of all the + END_DOC + implicit none + integer :: i,j,istate,k,l + double precision :: accu,hij + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + accu = 0.d0 + do istate = 1, N_states + do k = 1, ionic_index(i,0) + do l = 1, ionic_index(j,0) + hij = dressing_ref_hamiltonian(ionic_index(i,k),ionic_index(j,l),istate) +! accu += psi_ref_coef(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * & +! psi_ref_coef(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + accu += psi_ref_coef_dressed(ionic_index(i,k),istate) * normalization_factor_ionic_dressed(i,istate) * & + psi_ref_coef_dressed(ionic_index(j,l),istate) * normalization_factor_ionic_dressed(j,istate) * hij + enddo + enddo + H_OVB_dressing(i,j,istate) = accu + enddo + enddo + enddo + END_PROVIDER + + + + BEGIN_PROVIDER [double precision, H_OVB_total_dressed, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of all the + END_DOC + implicit none + integer :: i,j,istate + double precision :: accu,hij + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + do istate = 1, N_states + H_OVB_total_dressed(i,j,istate) = H_OVB_dressing(i,j,istate) + H_OVB_naked(i,j,istate) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, normalization_factor_ionic_dressed, (min_number_ionic:max_number_ionic, N_states) ] + implicit none + integer :: i,j,istate,k + double precision :: accu + do j = min_number_ionic, max_number_ionic + do istate = 1, N_states + accu = 0.d0 + do k = 1, ionic_index(j,0) + accu += psi_ref_coef_dressed(ionic_index(j,k),istate) **2 + enddo + normalization_factor_ionic_dressed(j,istate) = 1.d0/dsqrt(accu) + enddo + enddo + + END_PROVIDER diff --git a/plugins/OVB_effective_Hamiltonian/README.rst b/plugins/OVB_effective_Hamiltonian/README.rst new file mode 100644 index 00000000..a21ffcc6 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/README.rst @@ -0,0 +1,13 @@ +========================= +OVB_effective_Hamiltonian +========================= +Dressing of the OVB matrix by use of the Dressed_Ref_Hamiltonian dressing + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f b/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f new file mode 100644 index 00000000..bbd29b8e --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f @@ -0,0 +1,101 @@ +program print + read_wf = .True. + touch read_wf + call provide_all_stuffs +end +subroutine provide_all_stuffs + implicit none + provide ref_hamiltonian_matrix dressing_ref_hamiltonian + integer :: i,j,istate + double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + double precision, allocatable :: H_naked(:,:) + double precision, allocatable :: H_dressed(:,:) + double precision, allocatable :: H_print(:,:) + double precision :: accu_norm + allocate (H_dressed(max_number_ionic+1,max_number_ionic+1)) + allocate (H_print(min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)) + allocate (H_naked(max_number_ionic+1,max_number_ionic+1)) + allocate (psi_restart_ref_normalized(min_number_ionic:max_number_ionic)) + allocate (psi_ref_zeroth_order(min_number_ionic:max_number_ionic)) + print*,'# nuclear_repulsion = ',nuclear_repulsion + allocate (psi_ref_dressed(min_number_ionic:max_number_ionic)) + allocate (eigvalues(max_number_ionic+1)) + allocate (eigvectors(max_number_ionic+1,max_number_ionic+1)) + + do istate= 1, N_states + print*,'ISTATE = ',istate + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_print(i,j) = H_OVB_naked(j,i,istate) + enddo + enddo + do i = min_number_ionic,max_number_ionic + H_print(i,i) -= H_OVB_naked(min_number_ionic,min_number_ionic,istate) + enddo + + print*,'Ref Hamiltonian matrix emelent = ',H_OVB_naked(min_number_ionic,min_number_ionic,istate) + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:) + enddo + print*,'CAS MATRIX DRESSING' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_dressing(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + print*,'CAS MATRIX DRESSED ' + print*,'' + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_print(i,j) = H_OVB_total_dressed(j,i,istate) + enddo + enddo + do i = min_number_ionic,max_number_ionic + H_print(i,i) -= H_OVB_total_dressed(min_number_ionic,min_number_ionic,istate) + enddo + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:) + enddo + print*,'' + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + H_dressed(j+1,i+1) = H_OVB_total_dressed(i,j,istate) + H_naked(j+1,i+1) = H_OVB_naked(i,j,istate) + enddo + enddo + + call lapack_diagd(eigvalues,eigvectors,H_naked,max_number_ionic+1,max_number_ionic+1) + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = min_number_ionic,max_number_ionic + psi_ref_zeroth_order(i) = eigvectors(i+1,istate) + enddo + + + + call lapack_diagd(eigvalues,eigvectors,H_dressed,max_number_ionic+1,max_number_ionic+1) + do i = min_number_ionic,max_number_ionic + psi_ref_dressed(i) = eigvectors(i+1,istate) + enddo + print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion + do i = min_number_ionic,max_number_ionic + write(*,'(10(F10.7 ,4X))') psi_ref_dressed(i)/psi_ref_dressed(min_number_ionic) ,psi_ref_zeroth_order(i)/psi_ref_zeroth_order(min_number_ionic) + enddo + enddo + + deallocate (H_dressed) + deallocate (H_naked) + deallocate (psi_restart_ref_normalized) + deallocate (psi_ref_zeroth_order) + deallocate (psi_ref_dressed) + + deallocate (eigvalues) + deallocate (eigvectors) + +end diff --git a/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f b/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f new file mode 100644 index 00000000..4398b152 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f @@ -0,0 +1,90 @@ +program save_wf + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + use bitmasks + integer :: i,j,k,l + integer(bit_kind), allocatable :: psi_save_final(:,:,:) + double precision, allocatable :: psi_coef_save_final(:,:) + integer :: index_ref_determinants_save(psi_det_size) + integer :: n_det_ref_determinants_save + integer :: index_non_ref_determinants_save(psi_det_size) + integer :: n_det_non_ref_determinants_save + + integer :: n_det_save_final + integer :: number_of_particles + n_det_ref_determinants_save = 0 + integer :: ionic_level + ionic_level = 1 + do i = 1, ionic_index(ionic_level,0) ! number of determinants in the ref wf that are neutrals + n_det_ref_determinants_save +=1 + index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(ionic_level,i) + enddo + ! save all the 1p determinants in order to have the single excitations + ! on the top of the neutral structures + n_det_non_ref_determinants_save = 0 + do i = 1, N_det_non_ref + if(number_of_particles(psi_non_ref(1,1,i))==1)then + n_det_non_ref_determinants_save +=1 + index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i + endif + enddo + print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save + print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save + n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save + allocate (psi_save_final(N_int,2,n_det_save_final)) + allocate (psi_coef_save_final(n_det_save_final,1)) + integer :: n_det_tmp + n_det_tmp = 0 + do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i)) + enddo + psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1) + enddo + pause + do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i)) + enddo + accu = 0.d0 + double precision :: t_ik,hij + do j = 1, n_det_ref_determinants_save + call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij) + t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i)) + accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik + enddo + psi_coef_save_final(n_det_tmp,1) = accu + enddo + double precision :: accu + accu = 0.d0 + do i = 1, n_det_save_final + accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1) + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, n_det_save_final + psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1) + enddo + + do i = 1, n_det_save_final + print*,'' + print*,'Det' + call debug_det(psi_save_final(1,1,i),N_int) + print*,'coef = ',psi_coef_save_final(i,1) + enddo + + call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final) + deallocate (psi_save_final) + deallocate (psi_coef_save_final) + + +end diff --git a/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f b/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f new file mode 100644 index 00000000..a4b83d33 --- /dev/null +++ b/plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f @@ -0,0 +1,88 @@ +program save_wf + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + use bitmasks + integer :: i,j,k,l + integer(bit_kind), allocatable :: psi_save_final(:,:,:) + double precision, allocatable :: psi_coef_save_final(:,:) + integer :: index_ref_determinants_save(psi_det_size) + integer :: n_det_ref_determinants_save + integer :: index_non_ref_determinants_save(psi_det_size) + integer :: n_det_non_ref_determinants_save + + integer :: n_det_save_final + integer :: number_of_particles + n_det_ref_determinants_save = 0 + do i = 1, ionic_index(0,0) ! number of determinants in the ref wf that are neutrals + n_det_ref_determinants_save +=1 + index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(0,i) + enddo + ! save all the 1p determinants in order to have the single excitations + ! on the top of the neutral structures + n_det_non_ref_determinants_save = 0 + do i = 1, N_det_non_ref + if(number_of_particles(psi_non_ref(1,1,i))==1)then + n_det_non_ref_determinants_save +=1 + index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i + endif + enddo + print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save + print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save + n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save + allocate (psi_save_final(N_int,2,n_det_save_final)) + allocate (psi_coef_save_final(n_det_save_final,1)) + integer :: n_det_tmp + n_det_tmp = 0 + do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i)) + enddo + psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1) + enddo + pause + do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final + n_det_tmp +=1 + do j = 1, N_int + psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i)) + psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i)) + enddo + accu = 0.d0 + double precision :: t_ik,hij + do j = 1, n_det_ref_determinants_save + call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij) + t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i)) + accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik + enddo + psi_coef_save_final(n_det_tmp,1) = accu + enddo + double precision :: accu + accu = 0.d0 + do i = 1, n_det_save_final + accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1) + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, n_det_save_final + psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1) + enddo + + do i = 1, n_det_save_final + print*,'' + print*,'Det' + call debug_det(psi_save_final(1,1,i),N_int) + print*,'coef = ',psi_coef_save_final(i,1) + enddo + + call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final) + deallocate (psi_save_final) + deallocate (psi_coef_save_final) + + +end diff --git a/plugins/loc_cele/loc.f b/plugins/loc_cele/loc.f index e2439b7f..575932a3 100644 --- a/plugins/loc_cele/loc.f +++ b/plugins/loc_cele/loc.f @@ -17,7 +17,7 @@ C data small/1.d-6/ zprt=.true. - niter=100 + niter=500 conv=1.d-8 write (6,5) n,m,conv diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f index 12f90b64..e9c26f9d 100644 --- a/plugins/loc_cele/loc_cele.irp.f +++ b/plugins/loc_cele/loc_cele.irp.f @@ -9,7 +9,7 @@ ! id1=max is the number of MO in a given symmetry. END_DOC - integer id1 + integer id1,i_atom,shift,shift_h parameter (id1=300) @@ -92,7 +92,7 @@ - nrot(1) = 6 ! number of orbitals to be localized + nrot(1) = 6 ! number of orbitals to be localized integer :: index_rot(1000,1) @@ -101,12 +101,30 @@ cmoref = 0.d0 ! Definition of the index of the MO to be rotated - irot(1,1) = 20 ! the first mo to be rotated is the 19 th MO - irot(2,1) = 21 ! the first mo to be rotated is the 20 th MO - irot(3,1) = 22 ! etc.... - irot(4,1) = 23 ! - irot(5,1) = 24 ! - irot(6,1) = 25 ! +! irot(2,1) = 21 ! the first mo to be rotated is the 21 th MO +! irot(3,1) = 22 ! etc.... +! irot(4,1) = 23 ! +! irot(5,1) = 24 ! +! irot(6,1) = 25 ! +! do i = 1,12 +! irot(i,1) = i+6 +! enddo + irot(1,1) = 5 + irot(2,1) = 6 + irot(3,1) = 7 + irot(4,1) = 8 + irot(5,1) = 9 + irot(6,1) = 10 + do i = 1, nrot(1) + print*,'irot(i,1) = ',irot(i,1) + enddo + pause + cmoref(4,1,1) = 1.d0 ! 2S function + cmoref(5,2,1) = 1.d0 ! 2S function + cmoref(6,3,1) = 1.d0 ! 2S function + cmoref(19,4,1) = 1.d0 ! 2S function + cmoref(20,5,1) = 1.d0 ! 2S function + cmoref(21,6,1) = 1.d0 ! 2S function ! you define the guess vectors that you want ! the new MO to be close to @@ -120,23 +138,221 @@ ! own guess vectors for the MOs ! The new MOs are provided in output ! in the same order than the guess MOs - cmoref(3,1,1) = 1.d0 ! - cmoref(12,1,1) = 1.d0 ! + + ! C-C bonds + ! 1-2 +! i_atom = 1 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,1,1) = -0.012d0 ! 2S function +! cmoref(2+shift,1,1) = 0.18d0 ! +! cmoref(3+shift,1,1) = 0.1d0 ! - cmoref(21,2,1) = 1.d0 ! - cmoref(30,2,1) = 1.d0 ! +! cmoref(5+shift,1,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,1,1) = -0.1d0 ! 2pZ function - cmoref(39,3,1) = 1.d0 ! - cmoref(48,3,1) = 1.d0 ! +! i_atom = 2 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,1,1) = -0.012d0 ! 2S function +! cmoref(2+shift,1,1) = 0.18d0 ! +! cmoref(3+shift,1,1) = 0.1d0 ! - cmoref(3,4,1) = 1.d0 ! - cmoref(12,4,1) =-1.d0 ! +! cmoref(5+shift,1,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,1,1) = 0.1d0 ! 2pZ function - cmoref(21,5,1) = 1.d0 ! - cmoref(30,5,1) =-1.d0 ! - cmoref(39,6,1) = 1.d0 ! - cmoref(48,6,1) =-1.d0 ! +! ! 1-3 +! i_atom = 1 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,2,1) = -0.012d0 ! 2S function +! cmoref(2+shift,2,1) = 0.18d0 ! +! cmoref(3+shift,2,1) = 0.1d0 ! + +! cmoref(5+shift,2,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,2,1) = -0.1d0 ! 2pZ function + +! i_atom = 3 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,2,1) = -0.012d0 ! 2S function +! cmoref(2+shift,2,1) = 0.18d0 ! +! cmoref(3+shift,2,1) = 0.1d0 ! + +! cmoref(5+shift,2,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,2,1) = 0.1d0 ! 2pZ function + +! ! 4-6 +! i_atom = 4 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,3,1) = -0.012d0 ! 2S function +! cmoref(2+shift,3,1) = 0.18d0 ! +! cmoref(3+shift,3,1) = 0.1d0 ! + +! cmoref(5+shift,3,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,3,1) = -0.1d0 ! 2pZ function + +! i_atom = 6 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,3,1) = -0.012d0 ! 2S function +! cmoref(2+shift,3,1) = 0.18d0 ! +! cmoref(3+shift,3,1) = 0.1d0 ! + +! cmoref(5+shift,3,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,3,1) = 0.1d0 ! 2pZ function + + +! ! 6-5 +! i_atom = 6 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,4,1) = -0.012d0 ! 2S function +! cmoref(2+shift,4,1) = 0.18d0 ! +! cmoref(3+shift,4,1) = 0.1d0 ! + +! cmoref(5+shift,4,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,4,1) = 0.1d0 ! 2pZ function + +! i_atom = 5 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,4,1) = -0.012d0 ! 2S function +! cmoref(2+shift,4,1) = 0.18d0 ! +! cmoref(3+shift,4,1) = 0.1d0 ! + +! cmoref(5+shift,4,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,4,1) = -0.1d0 ! 2pZ function + + +! ! 2-4 +! i_atom = 2 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,5,1) = -0.012d0 ! 2S function +! cmoref(2+shift,5,1) = 0.18d0 ! +! cmoref(3+shift,5,1) = 0.1d0 ! + +! cmoref(6+shift,5,1) = 0.1d0 ! 2pZ function + +! i_atom = 4 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,5,1) = -0.012d0 ! 2S function +! cmoref(2+shift,5,1) = 0.18d0 ! +! cmoref(3+shift,5,1) = 0.1d0 ! + +! cmoref(6+shift,5,1) = -0.1d0 ! 2pZ function + + +! ! 3-5 +! i_atom = 3 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,6,1) = -0.012d0 ! 2S function +! cmoref(2+shift,6,1) = 0.18d0 ! +! cmoref(3+shift,6,1) = 0.1d0 ! + +! cmoref(6+shift,6,1) = 0.1d0 ! 2pZ function + +! i_atom = 5 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,6,1) = -0.012d0 ! 2S function +! cmoref(2+shift,6,1) = 0.18d0 ! +! cmoref(3+shift,6,1) = 0.1d0 ! + +! cmoref(6+shift,6,1) = -0.1d0 ! 2pZ function + +! ! C-H bonds +! ! 2-7 +! i_atom = 2 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,7,1) = -0.012d0 ! 2S function +! cmoref(2+shift,7,1) = 0.18d0 ! +! cmoref(3+shift,7,1) = 0.1d0 ! + +! cmoref(5+shift,7,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,7,1) = 0.1d0 ! 2pZ function +! +! i_atom = 7 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,7,1) = 0.12d0 ! 1S function + +! ! 4-10 +! i_atom = 4 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,8,1) = -0.012d0 ! 2S function +! cmoref(2+shift,8,1) = 0.18d0 ! +! cmoref(3+shift,8,1) = 0.1d0 ! + +! cmoref(5+shift,8,1) = -0.1d0 ! 2pX function +! cmoref(6+shift,8,1) = -0.1d0 ! 2pZ function +! +! i_atom = 10 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,8,1) = 0.12d0 ! 1S function + +! ! 5-11 +! i_atom = 5 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,9,1) = -0.012d0 ! 2S function +! cmoref(2+shift,9,1) = 0.18d0 ! +! cmoref(3+shift,9,1) = 0.1d0 ! + +! cmoref(5+shift,9,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,9,1) = -0.1d0 ! 2pZ function +! +! i_atom = 11 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,9,1) = 0.12d0 ! 1S function + +! ! 3-8 +! i_atom = 3 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,10,1) = -0.012d0 ! 2S function +! cmoref(2+shift,10,1) = 0.18d0 ! +! cmoref(3+shift,10,1) = 0.1d0 ! +! +! cmoref(5+shift,10,1) = 0.1d0 ! 2pX function +! cmoref(6+shift,10,1) = 0.1d0 ! 2pZ function +! +! i_atom = 8 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,10,1) = 0.12d0 ! 1S function + +! ! 1-9 +! i_atom = 1 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,11,1) = -0.012d0 ! 2S function +! cmoref(2+shift,11,1) = 0.18d0 ! +! cmoref(3+shift,11,1) = 0.1d0 ! +! +! cmoref(6+shift,11,1) = 0.1d0 ! 2pZ function + +! i_atom = 9 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,11,1) = 0.12d0 ! 1S function + +! +! ! 6-12 +! i_atom = 6 +! shift = (i_atom -1) * 15 +! cmoref(1+shift,12,1) = -0.012d0 ! 2S function +! cmoref(2+shift,12,1) = 0.18d0 ! +! cmoref(3+shift,12,1) = 0.1d0 ! +! +! cmoref(6+shift,12,1) = -0.1d0 ! 2pZ function + +! i_atom = 12 +! shift_h = (6-1) * 15 + (i_atom - 6)*5 +! cmoref(1+shift_h,12,1) = 0.12d0 ! 1S function +! cmoref(12,1,1) = 1.d0 ! + +! cmoref(21,2,1) = 1.d0 ! +! cmoref(30,2,1) = 1.d0 ! + +! cmoref(39,3,1) = 1.d0 ! +! cmoref(48,3,1) = 1.d0 ! + +! cmoref(3,4,1) = 1.d0 ! +! cmoref(12,4,1) =-1.d0 ! + +! cmoref(21,5,1) = 1.d0 ! +! cmoref(30,5,1) =-1.d0 ! + +! cmoref(39,6,1) = 1.d0 ! +! cmoref(48,6,1) =-1.d0 ! @@ -146,48 +362,11 @@ - do isym=1,nsym - - if (nrot(isym).eq.0) cycle - - do i=1,ao_num - - s(i,i,isym)=1.d0 - - do j=1,ao_num - - if (i.ne.j) s(i,j,isym)=0.d0 - - ddum(i,j)=0.d0 - - do k=1,nmo(isym) - - ddum(i,j)=ddum(i,j)+cmo(i,k,isym)*cmo(j,k,isym) - - enddo - + do i = 1, ao_num + do j = 1, ao_num + s(i,j,1) = ao_overlap(i,j) enddo - - enddo - - call dgesv(ao_num,ao_num,ddum,id1,ipiv,s(1,1,isym),id1,info) - - if (info.ne.0) then - - write (6,*) 'Something wrong in dgsev',isym - - stop - - endif - - - enddo - - - - - !Now big loop over symmetry diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index f83d8b41..ca2be5d6 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -23,7 +23,12 @@ deinit_thread skip init_main filter_integrals +filter2p filter2h2p +filter1h +filter1p +only_2p_single +only_2p_double filter_only_1h1p_single filter_only_1h1p_double filterhole @@ -44,7 +49,7 @@ class H_apply(object): self.selection_pt2 = None self.perturbation = None - + self.do_double_exc = do_double_exc #s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) & s["omp_parallel"] = """ PROVIDE elec_num_tab !$OMP PARALLEL DEFAULT(SHARED) & @@ -152,6 +157,32 @@ class H_apply(object): self["filterparticle"] = """ if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle """ + def filter_1h(self): + self["filter1h"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h(hole)) cycle + """ + def filter_2p(self): + self["filter2p"] = """ +! ! DIR$ FORCEINLINE + if (is_a_2p(hole)) cycle + """ + def filter_1p(self): + self["filter0p"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1p(hole)) cycle + """ + + def filter_only_2p(self): + self["only_2p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_2p(hole).eq..False.) cycle + """ + self["only_2p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_2p(key).eq..False.) cycle + """ + def filter_only_1h1p(self): self["filter_only_1h1p_single"] = """ @@ -216,10 +247,23 @@ class H_apply(object): self.data["initialization"] = """ PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit """ - self.data["keys_work"] = """ - call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) - """%(pert,) + if self.do_double_exc == True: + self.data["keys_work"] = """ +! if(check_double_excitation)then + call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) +! else +! call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & +! sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) +! endif + """%(pert,pert) + else: + self.data["keys_work"] = """ + call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) + """%(pert) + + self.data["finalization"] = """ """ self.data["copy_buffer"] = "" diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 50a95312..d6d8fcb0 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -399,18 +399,6 @@ END_PROVIDER enddo END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Bitmask of the core orbitals that are never excited in post CAS method - END_DOC - integer :: i,j - do i = 1, N_int - core_bitmask(i,1) = iand(ref_bitmask(i,1),reunion_of_bitmask(i,1)) - core_bitmask(i,2) = iand(ref_bitmask(i,2),reunion_of_bitmask(i,2)) - enddo - END_PROVIDER - BEGIN_PROVIDER [integer, list_core, (n_core_orb)] BEGIN_DOC ! List of the core orbitals that are never excited in post CAS method @@ -426,20 +414,21 @@ END_PROVIDER enddo END_PROVIDER - BEGIN_PROVIDER [ integer, n_core_orb ] + BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)] +&BEGIN_PROVIDER [ integer, n_core_orb] implicit none BEGIN_DOC - ! Number of core orbitals that are never excited in post CAS method + ! Core orbitals bitmask END_DOC - logical :: exists - integer :: j,i - integer :: i_hole,i_part,i_gen - + integer :: i,j n_core_orb = 0 - do j = 1, N_int - n_core_orb += popcnt(core_bitmask(j,1)) + do i = 1, N_int + core_bitmask(i,1) = xor(closed_shell_ref_bitmask(i,1),reunion_of_cas_inact_bitmask(i,1)) + core_bitmask(i,2) = xor(closed_shell_ref_bitmask(i,2),reunion_of_cas_inact_bitmask(i,2)) + n_core_orb += popcnt(core_bitmask(i,1)) enddo - END_PROVIDER + print*,'n_core_orb = ',n_core_orb + END_PROVIDER BEGIN_PROVIDER [ integer, i_bitmask_gen ] @@ -490,3 +479,27 @@ BEGIN_PROVIDER [integer, list_act, (n_act_orb)] enddo END_PROVIDER + + BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] + implicit none + integer :: i,j + do i = 1, N_int + closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),cas_bitmask(i,1,1)) + closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),cas_bitmask(i,2,1)) + enddo + END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_cas_inact_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + reunion_of_cas_inact_bitmask(i,1) = ior(cas_bitmask(i,1,1),inact_bitmask(i,1)) + reunion_of_cas_inact_bitmask(i,2) = ior(cas_bitmask(i,2,1),inact_bitmask(i,2)) + enddo + END_PROVIDER + + diff --git a/src/Bitmask/find_hole.irp.f b/src/Bitmask/find_hole.irp.f new file mode 100644 index 00000000..bc74ce52 --- /dev/null +++ b/src/Bitmask/find_hole.irp.f @@ -0,0 +1,55 @@ +logical function is_the_hole_in_det(key_in,ispin,i_hole) + use bitmasks + ! returns true if the electron ispin is absent from i_hole + implicit none + integer, intent(in) :: i_hole,ispin + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer(bit_kind) :: key_tmp(N_int) + integer(bit_kind) :: itest(N_int) + integer :: i,j,k + do i = 1, N_int + itest(i) = 0_bit_kind + enddo + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + itest(k) = ibset(itest(k),j) + j = 0 + do i = 1, N_int + key_tmp(i) = iand(itest(i),key_in(i,ispin)) + j += popcnt(key_tmp(i)) + enddo + if(j==0)then + is_the_hole_in_det = .True. + else + is_the_hole_in_det = .False. + endif + +end + +logical function is_the_particl_in_det(key_in,ispin,i_particl) + use bitmasks + ! returns true if the electron ispin is absent from i_particl + implicit none + integer, intent(in) :: i_particl,ispin + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer(bit_kind) :: key_tmp(N_int) + integer(bit_kind) :: itest(N_int) + integer :: i,j,k + do i = 1, N_int + itest(i) = 0_bit_kind + enddo + k = ishft(i_particl-1,-bit_kind_shift)+1 + j = i_particl-ishft(k-1,bit_kind_shift)-1 + itest(k) = ibset(itest(k),j) + j = 0 + do i = 1, N_int + key_tmp(i) = iand(itest(i),key_in(i,ispin)) + j += popcnt(key_tmp(i)) + enddo + if(j==0)then + is_the_particl_in_det = .False. + else + is_the_particl_in_det = .True. + endif + +end diff --git a/src/Bitmask/modify_bitmasks.irp.f b/src/Bitmask/modify_bitmasks.irp.f new file mode 100644 index 00000000..51d77af6 --- /dev/null +++ b/src/Bitmask/modify_bitmasks.irp.f @@ -0,0 +1,280 @@ + +use bitmasks +subroutine initialize_bitmask_to_restart_ones + implicit none + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC + ! Initialization of the generators_bitmask to the restart bitmask + END_DOC + do i = 1, N_int + do k=1,N_generators_bitmask + do ispin=1,2 + generators_bitmask(i,ispin,s_hole ,k) = generators_bitmask_restart(i,ispin,s_hole ,k) + generators_bitmask(i,ispin,s_part ,k) = generators_bitmask_restart(i,ispin,s_part ,k) + generators_bitmask(i,ispin,d_hole1,k) = generators_bitmask_restart(i,ispin,d_hole1,k) + generators_bitmask(i,ispin,d_part1,k) = generators_bitmask_restart(i,ispin,d_part1,k) + generators_bitmask(i,ispin,d_hole2,k) = generators_bitmask_restart(i,ispin,d_hole2,k) + generators_bitmask(i,ispin,d_part2,k) = generators_bitmask_restart(i,ispin,d_part2,k) + enddo + enddo + enddo +end + + +subroutine modify_bitmasks_for_hole(i_hole) + implicit none + integer, intent(in) :: i_hole + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC +! modify the generators_bitmask in order that one can only excite +! the electrons occupying i_hole + END_DOC + + ! Set to Zero the holes + do k=1,N_generators_bitmask + do l = 1, 3 + i = index_holes_bitmask(l) + do ispin=1,2 + do j = 1, N_int + generators_bitmask(j,ispin,i,k) = 0_bit_kind + enddo + enddo + enddo + enddo + + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + do m = 1, N_generators_bitmask + do l = 1, 3 + i = index_holes_bitmask(l) + do ispin=1,2 + generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + enddo + enddo + enddo + +end + +subroutine modify_bitmasks_for_hole_in_out(i_hole) + implicit none + integer, intent(in) :: i_hole + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC +! modify the generators_bitmask in order that one can only excite +! the electrons occupying i_hole + END_DOC + + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + do m = 1, N_generators_bitmask + do l = 1, 3 + i = index_holes_bitmask(l) + do ispin=1,2 + generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + enddo + enddo + enddo + +end + +subroutine modify_bitmasks_for_particl(i_part) + implicit none + integer, intent(in) :: i_part + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC +! modify the generators_bitmask in order that one can only excite +! the electrons to the orbital i_part + END_DOC + + ! Set to Zero the particles + do k=1,N_generators_bitmask + do l = 1, 3 + i = index_particl_bitmask(l) + do ispin=1,2 + do j = 1, N_int + generators_bitmask(j,ispin,i,k) = 0_bit_kind + enddo + enddo + enddo + enddo + + k = ishft(i_part-1,-bit_kind_shift)+1 + j = i_part-ishft(k-1,bit_kind_shift)-1 + do m = 1, N_generators_bitmask + do l = 1, 3 + i = index_particl_bitmask(l) + do ispin=1,2 + generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j) + enddo + enddo + enddo + +end + + +subroutine set_bitmask_particl_as_input(input_bimask) + implicit none + integer(bit_kind), intent(in) :: input_bimask(N_int,2) + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC +! set the generators_bitmask for the particles +! as the input_bimask + END_DOC + + do k=1,N_generators_bitmask + do l = 1, 3 + i = index_particl_bitmask(l) + do ispin=1,2 + do j = 1, N_int + generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin) + enddo + enddo + enddo + enddo + touch generators_bitmask + +end + + +subroutine set_bitmask_hole_as_input(input_bimask) + implicit none + integer(bit_kind), intent(in) :: input_bimask(N_int,2) + integer :: i,j,k,l,m + integer :: ispin + BEGIN_DOC +! set the generators_bitmask for the holes +! as the input_bimask + END_DOC + + do k=1,N_generators_bitmask + do l = 1, 3 + i = index_holes_bitmask(l) + do ispin=1,2 + do j = 1, N_int + generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin) + enddo + enddo + enddo + enddo + touch generators_bitmask + +end + + +subroutine print_generators_bitmasks_holes + implicit none + integer :: i,j,k,l + integer(bit_kind),allocatable :: key_tmp(:,:) + + allocate(key_tmp(N_int,2)) + do l = 1, 3 + k = 1 + i = index_holes_bitmask(l) + do j = 1, N_int + key_tmp(j,1) = generators_bitmask(j,1,i,k) + key_tmp(j,2) = generators_bitmask(j,2,i,k) + enddo + print*,'' + print*,'index hole = ',i + call print_det(key_tmp,N_int) + print*,'' + enddo + deallocate(key_tmp) + +end + +subroutine print_generators_bitmasks_particles + implicit none + integer :: i,j,k,l + integer(bit_kind),allocatable :: key_tmp(:,:) + + allocate(key_tmp(N_int,2)) + do l = 1, 3 + k = 1 + i = index_particl_bitmask(l) + do j = 1, N_int + key_tmp(j,1) = generators_bitmask(j,1,i,k) + key_tmp(j,2) = generators_bitmask(j,2,i,k) + enddo + print*,'' + print*,'index particl ',i + call print_det(key_tmp,N_int) + print*,'' + enddo + deallocate(key_tmp) + +end + +subroutine print_generators_bitmasks_holes_for_one_generator(i_gen) + implicit none + integer, intent(in) :: i_gen + integer :: i,j,k,l + integer(bit_kind),allocatable :: key_tmp(:,:) + + allocate(key_tmp(N_int,2)) + do l = 1, 3 + k = i_gen + i = index_holes_bitmask(l) + do j = 1, N_int + key_tmp(j,1) = generators_bitmask(j,1,i,k) + key_tmp(j,2) = generators_bitmask(j,2,i,k) + enddo + print*,'' + print*,'index hole = ',i + call print_det(key_tmp,N_int) + print*,'' + enddo + deallocate(key_tmp) + +end + +subroutine print_generators_bitmasks_particles_for_one_generator(i_gen) + implicit none + integer, intent(in) :: i_gen + integer :: i,j,k,l + integer(bit_kind),allocatable :: key_tmp(:,:) + + allocate(key_tmp(N_int,2)) + do l = 1, 3 + k = i_gen + i = index_particl_bitmask(l) + do j = 1, N_int + key_tmp(j,1) = generators_bitmask(j,1,i,k) + key_tmp(j,2) = generators_bitmask(j,2,i,k) + enddo + print*,'' + print*,'index particl ',i + call print_det(key_tmp,N_int) + print*,'' + enddo + deallocate(key_tmp) + +end + + + BEGIN_PROVIDER [integer, index_holes_bitmask, (3)] + implicit none + BEGIN_DOC +! Index of the holes in the generators_bitmasks + END_DOC + index_holes_bitmask(1) = d_hole1 + index_holes_bitmask(2) = d_hole2 + index_holes_bitmask(3) = s_hole + + END_PROVIDER + + BEGIN_PROVIDER [integer, index_particl_bitmask, (3)] + implicit none + BEGIN_DOC +! Index of the holes in the generators_bitmasks + END_DOC + index_particl_bitmask(1) = d_part1 + index_particl_bitmask(2) = d_part2 + index_particl_bitmask(3) = s_part + + END_PROVIDER diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 3f99f705..86780430 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -166,6 +166,9 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl logical :: check_double_excitation logical :: is_a_1h1p + logical :: is_a_1h + logical :: is_a_1p + logical :: is_a_2p logical :: b_cycle check_double_excitation = .True. iproc = iproc_in @@ -300,6 +303,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl key(k,other_spin) = ibset(key(k,other_spin),l) $filter2h2p $filter_only_1h1p_double + $only_2p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -351,6 +355,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl key(k,ispin) = ibset(key(k,ispin),l) $filter2h2p $filter_only_1h1p_double + $only_2p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -422,6 +427,9 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato logical :: check_double_excitation logical :: is_a_1h1p + logical :: is_a_1h + logical :: is_a_1p + logical :: is_a_2p key_mask(:,:) = 0_bit_kind @@ -493,6 +501,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 $filterparticle hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + $only_2p_single + $filter1h + $filter1p + $filter2p $filter2h2p $filter_only_1h1p_single key_idx += 1 diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index b93b18b6..7a54bdbc 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -189,6 +189,39 @@ logical function is_connected_to(key,keys,Nint,Ndet) enddo end +logical function is_connected_to_by_mono(key,keys,Nint,Ndet) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: i, l + integer :: degree_x2 + + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + is_connected_to_by_mono = .false. + + do i=1,Ndet + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +& + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if (degree_x2 > 2) then + cycle + else + is_connected_to_by_mono = .true. + return + endif + enddo +end + integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) use bitmasks diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 5fe18c49..63ed7a92 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -642,6 +642,14 @@ subroutine read_dets(det,Nint,Ndet) end +subroutine save_ref_determinant + implicit none + use bitmasks + call save_wavefunction_general(1,1,ref_bitmask,1,1.d0) +end + + + subroutine save_wavefunction implicit none diff --git a/src/Determinants/save_HF_determinant.irp.f b/src/Determinants/save_HF_determinant.irp.f new file mode 100644 index 00000000..0b81f136 --- /dev/null +++ b/src/Determinants/save_HF_determinant.irp.f @@ -0,0 +1,5 @@ +program save_HF + implicit none + call save_ref_determinant + +end