diff --git a/config/ifort.cfg b/config/ifort.cfg index cc848cba..2efa9eac 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -51,7 +51,7 @@ FCFLAGS : -xSSE4.2 -O2 -ip -ftz # [DEBUG] FC : -g -traceback -FCFLAGS : -xSSE2 -C -fpe0 +FCFLAGS : -xSSE2 -C IRPF90_FLAGS : --openmp # OpenMP flags diff --git a/plugins/All_singles/H_apply.irp.f b/plugins/All_singles/H_apply.irp.f index d0a41f90..f34f003c 100644 --- a/plugins/All_singles/H_apply.irp.f +++ b/plugins/All_singles/H_apply.irp.f @@ -8,10 +8,9 @@ s.unset_skip() s.filter_only_1h1p() print s -s = H_apply("just_mono") +s = H_apply("just_mono",do_double_exc=False) s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() -s.unset_double_excitations() print s END_SHELL diff --git a/plugins/All_singles/all_singles.irp.f b/plugins/All_singles/all_singles.irp.f index 3b5c5cce..ad8648c7 100644 --- a/plugins/All_singles/all_singles.irp.f +++ b/plugins/All_singles/all_singles.irp.f @@ -15,7 +15,7 @@ subroutine routine integer :: N_st, degree double precision,allocatable :: E_before(:) integer :: n_det_before - N_st = N_states + N_st = N_states_diag allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) i = 0 print*,'N_det = ',N_det diff --git a/plugins/CAS_SD/H_apply.irp.f b/plugins/CAS_SD/H_apply.irp.f index aa393bc7..35c45fb6 100644 --- a/plugins/CAS_SD/H_apply.irp.f +++ b/plugins/CAS_SD/H_apply.irp.f @@ -20,22 +20,18 @@ 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/DDCI_selected/ddci.irp.f b/plugins/DDCI_selected/ddci.irp.f index 3fcb443b..08e17cdd 100644 --- a/plugins/DDCI_selected/ddci.irp.f +++ b/plugins/DDCI_selected/ddci.irp.f @@ -3,10 +3,10 @@ program ddci integer :: i,k - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) integer :: N_st, degree - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + N_st = N_states_diag + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) character*(64) :: perturbation pt2 = 1.d0 @@ -27,6 +27,8 @@ program ddci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_DDCI_selection(pt2, norm_pert, H_pert_diag, N_st) @@ -47,8 +49,21 @@ program ddci print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 + print *, 'E+PT2 = ', E_before+pt2 print *, '-----' + 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 call ezfio_set_ddci_selected_energy(CI_energy) if (abort_all) then exit diff --git a/plugins/FOBOCI/H_apply.irp.f b/plugins/FOBOCI/H_apply.irp.f index 0a488753..d8ab02f1 100644 --- a/plugins/FOBOCI/H_apply.irp.f +++ b/plugins/FOBOCI/H_apply.irp.f @@ -18,8 +18,22 @@ print s -s = H_apply("standard") +s = H_apply("only_1h2p") s.set_selection_pt2("epstein_nesbet") +s.filter_only_1h2p() +s.unset_skip() +print s + +s = H_apply("only_2h2p") +s.set_selection_pt2("epstein_nesbet") +s.filter_only_2h2p() +s.unset_skip() +print s + + +s = H_apply("only_2p") +s.set_selection_pt2("epstein_nesbet") +s.filter_only_2p() s.unset_skip() print s diff --git a/plugins/FOBOCI/H_apply_dressed_autonom.irp.f b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f index c5b0aa5c..5f09390e 100644 --- a/plugins/FOBOCI/H_apply_dressed_autonom.irp.f +++ b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f @@ -449,8 +449,8 @@ subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_g 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(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref integer :: i_generator, nmax diff --git a/plugins/FOBOCI/H_apply_dressed_autonom_bis.irp.f b/plugins/FOBOCI/H_apply_dressed_autonom_bis.irp.f new file mode 100644 index 00000000..a9b05fc7 --- /dev/null +++ b/plugins/FOBOCI/H_apply_dressed_autonom_bis.irp.f @@ -0,0 +1,385 @@ +subroutine H_apply_dressed_pert_monoexc_bis(key_in, hole_1,particl_1,i_generator,iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input ) + 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,n_det_input + double precision, intent(inout) :: delta_ij_generators_(n_det_input,n_det_input),E_ref + + integer(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input) + 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_bis(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input) + key_idx = 0 + endif + enddo ! ii +! !$OMP ENDDO NOWAIT + enddo ! ispin + call standard_dress_bis(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input) + + 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_pertk_single(delta_ij_, Ndet_generators,psi_det_generators_input,E_ref,psi_det_input,n_det_input) + 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,n_det_input + integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) + integer(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input) + double precision, intent(inout) :: delta_ij_(n_det_input,n_det_input),E_ref + + + + 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 + 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_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input) + 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(.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_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input) + 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 + + +subroutine standard_dress_bis(delta_ij_generators_,size_buffer,Ndet_generators,i_generator,n_selected,det_buffer,Nint,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc,n_det_input + integer, intent(in) :: Ndet_generators,size_buffer + double precision, intent(inout) :: delta_ij_generators_(n_det_input,n_det_input),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(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input) + integer :: i,j,k,m + integer :: new_size + integer :: degree(n_det_input) + integer :: idx(0:n_det_input) + 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(n_det_input) + double precision :: contrib,lambda_i,accu + integer :: number_of_holes,n_h, number_of_particles,n_p + + 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 + print* + n_h = number_of_holes(det_buffer(1,1,i)) + n_p = number_of_particles(det_buffer(1,1,i)) + print*,'n_h,n_p = ',n_h,n_p + call get_excitation_degree_vector(psi_det_input,det_buffer(1,1,i),degree,N_int,n_det_input,idx) + H_array = 0.d0 + do k=1,idx(0) + call i_h_j(det_buffer(1,1,i),psi_det_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) + + lambda_i = f + 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 + enddo +end + diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES index adeefe99..f6c0c1c4 100644 --- a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES +++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Generators_restart Selectors_no_sorted +Perturbation Selectors_no_sorted Hartree_Fock diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f index 924e1c72..c4f0b7ae 100644 --- a/plugins/FOBOCI/all_singles.irp.f +++ b/plugins/FOBOCI/all_singles.irp.f @@ -8,7 +8,7 @@ subroutine all_single 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 + threshold_davidson = 1.d-9 soft_touch threshold_davidson davidson_criterion i = 0 print*,'Doing all the mono excitations !' @@ -52,10 +52,12 @@ subroutine all_single enddo endif E_before = CI_energy + !!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO + exit enddo - threshold_davidson = 1.d-10 - soft_touch threshold_davidson davidson_criterion - call diagonalize_CI +! threshold_davidson = 1.d-8 +! soft_touch threshold_davidson davidson_criterion +! call diagonalize_CI print*,'Final Step ' print*,'N_det = ',N_det do i = 1, N_states_diag @@ -67,10 +69,250 @@ subroutine all_single do i = 1, 2 print*,'psi_coef = ',psi_coef(i,1) enddo -! call save_wavefunction deallocate(pt2,norm_pert,E_before) end +subroutine all_1h2p + 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 1h2P 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_only_1h2p(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 + + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo + deallocate(pt2,norm_pert,E_before) +end + +subroutine all_2h2p + 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 2h2P 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_only_2h2p(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 + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo + 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_only_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 + 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 + deallocate(pt2,norm_pert,E_before) + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo +end + subroutine all_single_no_1h_or_1p implicit none integer :: i,k @@ -126,7 +368,7 @@ subroutine all_single_no_1h_or_1p endif E_before = CI_energy enddo - threshold_davidson = 1.d-10 + threshold_davidson = 1.d-16 soft_touch threshold_davidson davidson_criterion call diagonalize_CI print*,'Final Step ' @@ -217,85 +459,6 @@ subroutine all_single_no_1h_or_1p_or_2p 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 diff --git a/plugins/FOBOCI/all_singles_split.irp.f b/plugins/FOBOCI/all_singles_split.irp.f index 57f6ebde..9ddf369a 100644 --- a/plugins/FOBOCI/all_singles_split.irp.f +++ b/plugins/FOBOCI/all_singles_split.irp.f @@ -5,7 +5,7 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N 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 + integer :: i,i_hole,j n_det_max_jacobi = 50 soft_touch n_det_max_jacobi do i = 1, n_inact_orb @@ -22,16 +22,170 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N 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 diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) enddo + + do i = 1, n_act_orb + i_hole = list_act(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 +! call diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) + call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input) + enddo + + do i = 1, n_virt_orb + i_hole = list_virt(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 +! call diagonalize_CI_SC2 +! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2) + 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_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) + implicit none + use bitmasks + integer, intent(in) :: i_particl + 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) + double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) + integer :: i,j + n_det_max_jacobi = 50 + soft_touch n_det_max_jacobi + + call all_single + + threshold_davidson = 1.d-12 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + + + + double precision, allocatable :: matrix_ref_1h_1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h2p(:,:) + double precision, allocatable :: psi_coef_ref_1h_1p(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_1h2p(:,:) + integer(bit_kind), allocatable :: psi_det_1h2p(:,:,:) + integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:) + integer :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + double precision :: hka + double precision,allocatable :: eigenvectors(:,:), eigenvalues(:) + + + call give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p) + + allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h2p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)) + allocate(psi_det_1h2p(N_int,2,n_det_1h2p), psi_coef_1h2p(n_det_1h2p,N_states)) + allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states)) + + call give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p) + + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka) + matrix_ref_1h_1p(i,j) = hka + enddo + enddo + matrix_ref_1h_1p_dressing_1h1p = 0.d0 + matrix_ref_1h_1p_dressing_1h2p = 0.d0 + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h2p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h2p,psi_coef_1h2p,n_det_1h2p) + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h1p,psi_coef_1h1p,n_det_1h1p) + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_1h2p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j) + enddo + enddo + + allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p)) + call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p) +!do j = 1, n_det_ref_1h_1p +! print*,'coef = ',eigenvectors(j,1) +!enddo + print*,'' + print*,'-----------------------' + print*,'-----------------------' + print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion + print*,'-----------------------' + ! Extract the + integer, allocatable :: index_generator(:) + integer :: n_det_generators_tmp,degree + n_det_generators_tmp = 0 + allocate(index_generator(n_det_ref_1h_1p)) + do i = 1, n_det_ref_1h_1p + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(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_generators_tmp .ne. n_det_generators)then + print*,'PB !!!' + print*,'if(n_det_generators_tmp .ne. n_det_genrators)then' + stop + endif + do i = 1, N_det_generators + print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1) + do j = 1, N_det_generators + dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j)) + dressing_matrix_1h2p(i,j) += matrix_ref_1h_1p_dressing_1h2p(index_generator(i),index_generator(j)) + enddo + enddo + print*,'-----------------------' + print*,'-----------------------' + + + deallocate(matrix_ref_1h_1p) + deallocate(matrix_ref_1h_1p_dressing_1h1p) + deallocate(matrix_ref_1h_1p_dressing_1h2p) + deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p) + deallocate(psi_det_1h2p, psi_coef_1h2p) + deallocate(psi_det_1h1p, psi_coef_1h1p) + deallocate(eigenvectors,eigenvalues) + deallocate(index_generator) + + +end + subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p) implicit none use bitmasks @@ -39,50 +193,168 @@ subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dr 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) double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) - integer :: i + integer :: i,j n_det_max_jacobi = 50 soft_touch n_det_max_jacobi - integer :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p - integer(bit_kind), allocatable :: psi_ref_out(:,:,:) - integer(bit_kind), allocatable :: psi_1h1p(:,:,:) - integer(bit_kind), allocatable :: psi_2h1p(:,:,:) - integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) - double precision, allocatable :: psi_ref_coef_out(:,:) - double precision, allocatable :: psi_coef_1h1p(:,:) - double precision, allocatable :: psi_coef_2h1p(:,:) - double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) -!call all_single_no_1h_or_1p call all_single threshold_davidson = 1.d-12 soft_touch threshold_davidson davidson_criterion call diagonalize_CI - call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) - 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)) - allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) - call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) - 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) - call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) - deallocate(psi_ref_out) - deallocate(psi_1h1p) - deallocate(psi_2h1p) - deallocate(psi_extra_1h_or_1p) - deallocate(psi_ref_coef_out) - deallocate(psi_coef_1h1p) - deallocate(psi_coef_2h1p) - deallocate(psi_coef_extra_1h_or_1p) + + + double precision, allocatable :: matrix_ref_1h_1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:) + double precision, allocatable :: matrix_ref_1h_1p_dressing_2h1p(:,:) + double precision, allocatable :: psi_coef_ref_1h_1p(:,:) + double precision, allocatable :: psi_coef_1h1p(:,:) + double precision, allocatable :: psi_coef_2h1p(:,:) + integer(bit_kind), allocatable :: psi_det_2h1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:) + integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:) + integer :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p + double precision :: hka + double precision,allocatable :: eigenvectors(:,:), eigenvalues(:) + + + call give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p) + + allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(matrix_ref_1h_1p_dressing_2h1p(n_det_ref_1h_1p,n_det_ref_1h_1p)) + allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)) + allocate(psi_det_2h1p(N_int,2,n_det_2h1p), psi_coef_2h1p(n_det_2h1p,N_states)) + allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states)) + + call give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p) + + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka) + matrix_ref_1h_1p(i,j) = hka + enddo + enddo + matrix_ref_1h_1p_dressing_1h1p = 0.d0 + matrix_ref_1h_1p_dressing_2h1p = 0.d0 + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_2h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_2h1p,psi_coef_2h1p,n_det_2h1p) + call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, & + psi_det_1h1p,psi_coef_1h1p,n_det_1h1p) + do i = 1, n_det_ref_1h_1p + do j = 1, n_det_ref_1h_1p + matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j) + enddo + enddo + + allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p)) + call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p) +!do j = 1, n_det_ref_1h_1p +! print*,'coef = ',eigenvectors(j,1) +!enddo + print*,'' + print*,'-----------------------' + print*,'-----------------------' + print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion + print*,'-----------------------' + ! Extract the + integer, allocatable :: index_generator(:) + integer :: n_det_generators_tmp,degree + n_det_generators_tmp = 0 + allocate(index_generator(n_det_ref_1h_1p)) + do i = 1, n_det_ref_1h_1p + do j = 1, N_det_generators + call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(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_generators_tmp .ne. n_det_generators)then + print*,'PB !!!' + print*,'if(n_det_generators_tmp .ne. n_det_genrators)then' + stop + endif + do i = 1, N_det_generators + print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1) + do j = 1, N_det_generators + dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j)) + dressing_matrix_2h1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(index_generator(i),index_generator(j)) + enddo + enddo + print*,'-----------------------' + print*,'-----------------------' + + + deallocate(matrix_ref_1h_1p) + deallocate(matrix_ref_1h_1p_dressing_1h1p) + deallocate(matrix_ref_1h_1p_dressing_2h1p) + deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p) + deallocate(psi_det_2h1p, psi_coef_2h1p) + deallocate(psi_det_1h1p, psi_coef_1h1p) + deallocate(eigenvectors,eigenvalues) + deallocate(index_generator) +!return +! + +!integer(bit_kind), allocatable :: psi_ref_out(:,:,:) +!integer(bit_kind), allocatable :: psi_1h1p(:,:,:) +!integer(bit_kind), allocatable :: psi_2h1p(:,:,:) +!integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) +!double precision, allocatable :: psi_ref_coef_out(:,:) +!double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) + +!call all_single_no_1h_or_1p + +!call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) +!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)) +!allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) +!call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) +!do i = 1, n_det_extra_1h_or_1p +! print*,'----' +! print*,'c = ',psi_coef_extra_1h_or_1p(i,1) +! call debug_det(psi_extra_1h_or_1p(1,1,i),N_int) +! print*,'----' +!enddo +!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) +!print*,'Dressing 1h1p ' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_1h1p(j,:) +!enddo + +!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) +!print*,'Dressing 2h1p ' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_2h1p(j,:) +!enddo + +!call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) +!print*,',dressing_matrix_extra_1h_or_1p' +!do j =1, N_det_generators +! print*,' dressing ',dressing_matrix_extra_1h_or_1p(j,:) +!enddo + + +!deallocate(psi_ref_out) +!deallocate(psi_1h1p) +!deallocate(psi_2h1p) +!deallocate(psi_extra_1h_or_1p) +!deallocate(psi_ref_coef_out) +!deallocate(psi_coef_1h1p) +!deallocate(psi_coef_2h1p) +!deallocate(psi_coef_extra_1h_or_1p) end @@ -208,56 +480,56 @@ subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) soft_touch n_det_max_jacobi end -subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) - implicit none - use bitmasks - integer, intent(in ) :: i_particl - 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) - double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) - integer :: i - n_det_max_jacobi = 50 - soft_touch n_det_max_jacobi - - integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p - integer(bit_kind), allocatable :: psi_ref_out(:,:,:) - integer(bit_kind), allocatable :: psi_1h1p(:,:,:) - integer(bit_kind), allocatable :: psi_1h2p(:,:,:) - integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) - double precision, allocatable :: psi_ref_coef_out(:,:) - double precision, allocatable :: psi_coef_1h1p(:,:) - double precision, allocatable :: psi_coef_1h2p(:,:) - double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) -!call all_single_no_1h_or_1p_or_2p - call all_single - - threshold_davidson = 1.d-12 - soft_touch threshold_davidson davidson_criterion - call diagonalize_CI - call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) - 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)) - allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) - call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) - 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) - call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & - psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) - - deallocate(psi_ref_out) - deallocate(psi_1h1p) - deallocate(psi_1h2p) - deallocate(psi_ref_coef_out) - deallocate(psi_coef_1h1p) - deallocate(psi_coef_1h2p) - -end +! subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) +! implicit none +! use bitmasks +! integer, intent(in ) :: i_particl +! 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) +! double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) +! integer :: i +! n_det_max_jacobi = 50 +! soft_touch n_det_max_jacobi +! +! integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p +! integer(bit_kind), allocatable :: psi_ref_out(:,:,:) +! integer(bit_kind), allocatable :: psi_1h1p(:,:,:) +! integer(bit_kind), allocatable :: psi_1h2p(:,:,:) +! integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) +! double precision, allocatable :: psi_ref_coef_out(:,:) +! double precision, allocatable :: psi_coef_1h1p(:,:) +! double precision, allocatable :: psi_coef_1h2p(:,:) +! double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) +!!!!call all_single_no_1h_or_1p_or_2p +! call all_single +! +! threshold_davidson = 1.d-12 +! soft_touch threshold_davidson davidson_criterion +! call diagonalize_CI +! call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) +! 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)) +! allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) +! call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) +! 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) +! call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & +! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) +! +! 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/collect_all_lmct.irp.f b/plugins/FOBOCI/collect_all_lmct.irp.f new file mode 100644 index 00000000..ebece3ed --- /dev/null +++ b/plugins/FOBOCI/collect_all_lmct.irp.f @@ -0,0 +1,436 @@ +use bitmasks + +subroutine collect_lmct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + print*,'COLLECTING THE PERTINENT LMCT (1h)' + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-2)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_mlct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + print*,'COLLECTING THE PERTINENT MLCT (1p)' + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-3)then + n_couples +=1 + hole_particle(n_couples,1) = iorb + hole_particle(n_couples,2) = jorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_lmct_mlct(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + double precision, allocatable :: tmp(:,:) + print*,'COLLECTING THE PERTINENT LMCT (1h)' + print*,'AND THE PERTINENT MLCT (1p)' + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + n_couples +=1 + hole_particle(n_couples,1) = iorb + hole_particle(n_couples,2) = jorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + +subroutine collect_1h1p(hole_particle,n_couples) + implicit none + integer, intent(out) :: hole_particle(1000,2), n_couples + BEGIN_DOC + ! Collect all the couple holes/particles of the important LMCT + ! hole_particle(i,1) = ith hole + ! hole_particle(i,2) = ith particle + ! n_couples is the number of important excitations + END_DOC + double precision, allocatable :: tmp(:,:) + print*,'COLLECTING THE PERTINENT 1h1p' + allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2))) + tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci + integer :: i,j,iorb,jorb + n_couples = 0 + do i = 1,n_virt_orb + iorb = list_virt(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + if(dabs(tmp(iorb,jorb)).gt.1.d-2)then + n_couples +=1 + hole_particle(n_couples,1) = jorb + hole_particle(n_couples,2) = iorb + print*,'DM' + print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb) + endif + enddo + enddo + deallocate(tmp) + print*,'number of meaning full couples of holes/particles ' + print*,'n_couples = ',n_couples + + +end + + + +subroutine set_lmct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + print*,'i_hole,i_particle 2 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + print*,'i_hole,i_particle 1 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_mlct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_mlct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + print*,'i_hole,i_particle 2 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + print*,'i_hole,i_particle 1 = ',i_hole,i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + print*,'i_ok = ',i_ok + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_lmct_mlct_to_generators_restart + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct_mlct(hole_particle,n_couples) + call set_generators_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_cas + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_cas(n,1,m) + key_tmp(n,2) = psi_cas(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det_generators(n,1,N_det_total) = key_tmp(n,1) + psi_det_generators(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + N_det_generators = N_det_total + do i = 1, N_det_generators + psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total)) + enddo + print*,'number of generators in total = ',N_det_generators + touch N_det_generators psi_coef_generators psi_det_generators +end + +subroutine set_lmct_mlct_to_psi_det + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_lmct_mlct(hole_particle,n_couples) + call set_psi_det_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_generators_restart + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + + N_det = N_det_total + integer :: k + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total)) + enddo + enddo + SOFT_TOUCH N_det psi_det psi_coef + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) +end + +subroutine set_1h1p_to_psi_det + implicit none + integer :: i,j,m,n,i_hole,i_particle + integer :: hole_particle(1000,2), n_couples + integer(bit_kind) :: key_tmp(N_int,2) + integer :: N_det_total,i_ok + + call collect_1h1p(hole_particle,n_couples) + call set_psi_det_to_generators_restart + N_det_total = N_det_generators_restart + do i = 1, n_couples + i_hole = hole_particle(i,1) + i_particle = hole_particle(i,2) + do m = 1, N_det_generators_restart + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + ! You excite the beta electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + + do n = 1, N_int + key_tmp(n,1) = psi_det_generators_restart(n,1,m) + key_tmp(n,2) = psi_det_generators_restart(n,2,m) + enddo + + ! You excite the alpha electron from i_hole to i_particle + call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok) + if(i_ok==1)then + N_det_total +=1 + do n = 1, N_int + psi_det(n,1,N_det_total) = key_tmp(n,1) + psi_det(n,2,N_det_total) = key_tmp(n,2) + enddo + endif + enddo + enddo + + N_det = N_det_total + integer :: k + do k = 1, N_states + do i = 1, N_det + psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total)) + 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/corr_energy_2h2p.irp.f b/plugins/FOBOCI/corr_energy_2h2p.irp.f new file mode 100644 index 00000000..8c4f2fe3 --- /dev/null +++ b/plugins/FOBOCI/corr_energy_2h2p.irp.f @@ -0,0 +1,363 @@ + BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_2h2p] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_2h2p = 0.d0 + corr_energy_2h2p_per_orb_ab = 0.d0 + corr_energy_2h2p_per_orb_aa = 0.d0 + corr_energy_2h2p_per_orb_bb = 0.d0 + do i = 1, n_inact_orb ! beta + i_hole = list_inact(i) + do k = 1, n_virt_orb ! beta + k_part = list_virt(k) + do j = 1, n_inact_orb ! alpha + j_hole = list_inact(j) + do l = 1, n_virt_orb ! alpha + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = (ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = hij*hij/delta_e + total_corr_e_2h2p += contrib + corr_energy_2h2p_per_orb_ab(i_hole) += contrib + corr_energy_2h2p_per_orb_ab(k_part) += contrib + enddo + enddo + enddo + enddo + + ! alpha alpha correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_virt_orb + k_part = list_virt(k) + do l = k+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h2p += contrib + corr_energy_2h2p_per_orb_aa(i_hole) += contrib + corr_energy_2h2p_per_orb_aa(k_part) += contrib + enddo + enddo + enddo + enddo + + ! beta beta correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_virt_orb + k_part = list_virt(k) + do l = k+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h2p += contrib + corr_energy_2h2p_per_orb_bb(i_hole) += contrib + corr_energy_2h2p_per_orb_bb(k_part) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_2h1p] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_2h1p = 0.d0 + corr_energy_2h1p_per_orb_ab = 0.d0 + corr_energy_2h1p_per_orb_aa = 0.d0 + corr_energy_2h1p_per_orb_bb = 0.d0 + do i = 1, n_inact_orb + i_hole = list_inact(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_2h1p += contrib + corr_energy_2h1p_per_orb_ab(i_hole) += contrib + corr_energy_2h1p_per_orb_ab(l_part) += contrib + enddo + enddo + enddo + enddo + + ! Alpha Alpha spin correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = 1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_2h1p += contrib + corr_energy_2h1p_per_orb_aa(i_hole) += contrib + corr_energy_2h1p_per_orb_aa(l_part) += contrib + enddo + enddo + enddo + enddo + + ! Beta Beta correlation energy + do i = 1, n_inact_orb + i_hole = list_inact(i) + do j = i+1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = 1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_2h1p += contrib + corr_energy_2h1p_per_orb_bb(i_hole) += contrib + corr_energy_2h1p_per_orb_aa(l_part) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_1h2p] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_1h2p = 0.d0 + corr_energy_1h2p_per_orb_ab = 0.d0 + corr_energy_1h2p_per_orb_aa = 0.d0 + corr_energy_1h2p_per_orb_bb = 0.d0 + do i = 1, n_virt_orb + i_hole = list_virt(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_ab(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + enddo + enddo + enddo + enddo + + ! Alpha Alpha correlation energy + do i = 1, n_virt_orb + i_hole = list_virt(i) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = i+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + + key_tmp = ref_bitmask + ispin = 1 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_aa(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + enddo + enddo + enddo + enddo + + ! Beta Beta correlation energy + do i = 1, n_virt_orb + i_hole = list_virt(i) + do j = 1, n_inact_orb + j_hole = list_inact(j) + do k = 1, n_act_orb + k_part = list_act(k) + do l = i+1,n_virt_orb + l_part = list_virt(l) + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 2 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + hij = hij - exc + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + total_corr_e_1h2p += contrib + corr_energy_1h2p_per_orb_bb(i_hole) += contrib + corr_energy_1h2p_per_orb_ab(j_hole) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, corr_energy_1h1p_spin_flip_per_orb, (mo_tot_num)] +&BEGIN_PROVIDER [ double precision, total_corr_e_1h1p_spin_flip] + use bitmasks + implicit none + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i,j,k,l + integer :: i_hole,j_hole,k_part,l_part + double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib + double precision :: diag_H_mat_elem + integer :: i_ok,ispin + ! Alpha - Beta correlation energy + total_corr_e_1h1p_spin_flip = 0.d0 + corr_energy_1h1p_spin_flip_per_orb = 0.d0 + do i = 1, n_inact_orb + i_hole = list_inact(i) + do k = 1, n_act_orb + k_part = list_act(k) + do j = 1, n_act_orb + j_hole = list_act(j) + do l = 1, n_virt_orb + l_part = list_virt(l) + + key_tmp = ref_bitmask + ispin = 2 + call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok) + if(i_ok .ne.1)cycle + ispin = 1 + call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok) + if(i_ok .ne.1)cycle + delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) + + hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) + contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + + total_corr_e_1h1p_spin_flip += contrib + corr_energy_1h1p_spin_flip_per_orb(i_hole) += contrib + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f index a4c6b652..83955e61 100644 --- a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f +++ b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f @@ -3,6 +3,7 @@ subroutine diag_inactive_virt_and_update_mos integer :: i,j,i_inact,j_inact,i_virt,j_virt double precision :: tmp(mo_tot_num_align,mo_tot_num) character*(64) :: label + print*,'Diagonalizing the occ and virt Fock operator' tmp = 0.d0 do i = 1, mo_tot_num tmp(i,i) = Fock_matrix_mo(i,i) @@ -33,3 +34,50 @@ subroutine diag_inactive_virt_and_update_mos end + +subroutine diag_inactive_virt_new_and_update_mos + implicit none + integer :: i,j,i_inact,j_inact,i_virt,j_virt,k,k_act + double precision :: tmp(mo_tot_num_align,mo_tot_num),accu,get_mo_bielec_integral_schwartz + 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) + accu =0.d0 + do k = 1, n_act_orb + k_act = list_act(k) + accu += get_mo_bielec_integral_schwartz(i_inact,k_act,j_inact,k_act,mo_integrals_map) + accu -= get_mo_bielec_integral_schwartz(i_inact,k_act,k_act,j_inact,mo_integrals_map) + enddo + tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact) + accu + tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact) + accu + 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) + accu =0.d0 + do k = 1, n_act_orb + k_act = list_act(k) + accu += get_mo_bielec_integral_schwartz(i_virt,k_act,j_virt,k_act,mo_integrals_map) + enddo + tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt) - accu + tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt) - accu + 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 index 2f662f4d..9df94140 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -85,33 +85,6 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen 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 diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f new file mode 100644 index 00000000..1b134733 --- /dev/null +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -0,0 +1,35 @@ +program foboscf + implicit none + no_oa_or_av_opt = .True. + touch no_oa_or_av_opt + call run_prepare + call routine_fobo_scf + call save_mos + +end + +subroutine run_prepare + implicit none + call damping_SCF + call diag_inactive_virt_and_update_mos +end + +subroutine routine_fobo_scf + implicit none + integer :: i,j + print*,'' + print*,'' + character*(64) :: label + label = "Natural" + do i = 1, 5 + call FOBOCI_lmct_mlct_old_thr + call save_osoci_natural_mos + call damping_SCF + call diag_inactive_virt_and_update_mos + call clear_mo_map + call provide_properties + enddo + + + +end diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index 087f791b..b406284f 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -36,6 +36,10 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'' print*,'' print*,'DOING FIRST LMCT !!' + integer(bit_kind) , allocatable :: zero_bitmask(:,:) + integer(bit_kind) , allocatable :: psi_singles(:,:,:) + double precision, allocatable :: psi_singles_coef(:,:) + allocate( zero_bitmask(N_int,2) ) do i = 1, n_inact_orb integer :: i_hole_osoci i_hole_osoci = list_inact(i) @@ -54,24 +58,85 @@ subroutine FOBOCI_lmct_mlct_old_thr 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)) + dressing_matrix = 0.d0 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) + hkl = dressing_matrix(1,1) + do k = 1, N_det_generators + dressing_matrix(k,k) = dressing_matrix(k,k) - hkl + enddo + print*,'Naked matrix' + do k = 1, N_det_generators + write(*,'(100(F12.5,X))')dressing_matrix(k,:) + enddo + + ! Do all the single excitations on top of the CAS and 1h determinants + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single + +! ! Change the mask of the holes and particles to perform all the +! ! double excitations that starts from the active space in order +! ! to introduce the Coulomb hole in the active space +! ! These are the 1h2p excitations that have the i_hole_osoci hole in common +! ! and the 2p if there is more than one electron in the active space +! do k = 1, N_int +! zero_bitmask(k,1) = 0_bit_kind +! zero_bitmask(k,2) = 0_bit_kind +! enddo +! ! hole is possible only in the orbital i_hole_osoci +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int) +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int) +! ! and in the active space +! do k = 1, n_act_orb +! call set_bit_to_integer(list_act(k),zero_bitmask(1,1),N_int) +! call set_bit_to_integer(list_act(k),zero_bitmask(1,2),N_int) +! enddo +! call set_bitmask_hole_as_input(zero_bitmask) + +! call set_bitmask_particl_as_input(reunion_of_bitmask) + +! call all_1h2p +! call diagonalize_CI_SC2 +! call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) + +! ! Change the mask of the holes and particles to perform all the +! ! double excitations that from the orbital i_hole_osoci +! do k = 1, N_int +! zero_bitmask(k,1) = 0_bit_kind +! zero_bitmask(k,2) = 0_bit_kind +! enddo +! ! hole is possible only in the orbital i_hole_osoci +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int) +! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int) +! call set_bitmask_hole_as_input(zero_bitmask) + +! call set_bitmask_particl_as_input(reunion_of_bitmask) + +! call set_psi_det_to_generators +! call all_2h2p +! call diagonalize_CI_SC2 + double precision :: hkl + call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) + hkl = dressing_matrix(1,1) + do k = 1, N_det_generators + dressing_matrix(k,k) = dressing_matrix(k,k) - hkl + enddo + print*,'Dressed matrix' + do k = 1, N_det_generators + write(*,'(100(F12.5,X))')dressing_matrix(k,:) + enddo +! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) 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) @@ -132,24 +197,6 @@ subroutine FOBOCI_lmct_mlct_old_thr 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 diff --git a/plugins/FOBOCI/foboci_reunion.irp.f b/plugins/FOBOCI/foboci_reunion.irp.f new file mode 100644 index 00000000..fcfaff58 --- /dev/null +++ b/plugins/FOBOCI/foboci_reunion.irp.f @@ -0,0 +1,18 @@ +program osoci_program +implicit none + do_it_perturbative = .True. + touch do_it_perturbative + 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 set_lmct_mlct_to_psi_det + call diagonalize_CI + call save_wavefunction + + +end + diff --git a/plugins/FOBOCI/generators_restart_save.irp.f b/plugins/FOBOCI/generators_restart_save.irp.f index dca4c901..09d4aa2b 100644 --- a/plugins/FOBOCI/generators_restart_save.irp.f +++ b/plugins/FOBOCI/generators_restart_save.irp.f @@ -1,126 +1,74 @@ -use bitmasks +use bitmasks + BEGIN_PROVIDER [ integer, N_det_generators_restart ] implicit none BEGIN_DOC - ! Number of determinants in the wave function + ! Read the wave function END_DOC - logical :: exists - character*64 :: label + integer :: i 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) + double precision :: norm + if(ifirst == 0)then + call ezfio_get_determinants_n_det(N_det_generators_restart) ifirst = 1 -!endif + else + print*,'PB in generators_restart restart !!!' + endif + call write_int(output_determinants,N_det_generators_restart,'Number of generators_restart') END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,N_det_generators_restart) ] &BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ] implicit none BEGIN_DOC - ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file - ! is empty + ! read wf + ! END_DOC - integer :: i - logical :: exists - character*64 :: label - + integer :: i, k 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 + double precision, allocatable :: psi_coef_read(:,:) + if(ifirst == 0)then + call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) + do k = 1, N_int + ref_generators_restart(k,1) = psi_det_generators_restart(k,1,1) + ref_generators_restart(k,2) = psi_det_generators_restart(k,2,1) + enddo + 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 ifirst = 1 -!endif + deallocate(psi_coef_read) + else + print*,'PB in generators_restart restart !!!' + 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 - - - +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 END_PROVIDER +BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] + implicit none + BEGIN_DOC + ! Memo to skip useless selectors + END_DOC + select_max = huge(1.d0) +END_PROVIDER + + BEGIN_PROVIDER [ integer, N_det_generators ] +&BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,10000) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (10000,N_states) ] + +END_PROVIDER diff --git a/plugins/FOBOCI/hcc_1h1p.irp.f b/plugins/FOBOCI/hcc_1h1p.irp.f new file mode 100644 index 00000000..66cf2fd4 --- /dev/null +++ b/plugins/FOBOCI/hcc_1h1p.irp.f @@ -0,0 +1,83 @@ +program test_sc2 + implicit none + read_wf = .True. + touch read_wf + call routine + + +end + +subroutine routine + implicit none + double precision, allocatable :: energies(:),diag_H_elements(:) + double precision, allocatable :: H_matrix(:,:) + allocate(energies(N_states),diag_H_elements(N_det)) + call diagonalize_CI + call test_hcc + call test_mulliken +! call SC2_1h1p(psi_det,psi_coef,energies, & +! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + allocate(H_matrix(N_det,N_det)) + call SC2_1h1p_full(psi_det,psi_coef,energies, & + H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + deallocate(H_matrix) + integer :: i,j + double precision :: accu,coef_hf +! coef_hf = 1.d0/psi_coef(1,1) +! do i = 1, N_det +! psi_coef(i,1) *= coef_hf +! enddo + touch psi_coef + call pouet +end + +subroutine pouet + implicit none + double precision :: accu,coef_hf +! provide one_body_dm_mo_alpha one_body_dm_mo_beta +! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int) +! touch one_body_dm_mo_alpha one_body_dm_mo_beta + call test_hcc + call test_mulliken +! call save_wavefunction + +end + +subroutine test_hcc + implicit none + double precision :: accu + integer :: i,j + print*,'Z AU GAUSS MHZ cm^-1' + do i = 1, nucl_num + write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) + enddo + +end + +subroutine test_mulliken + double precision :: accu + integer :: i + integer :: j + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),mulliken_spin_densities(i) + accu += mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu +!print*,'AO SPIN POPULATIONS' + accu = 0.d0 +!do i = 1, ao_num +! accu += spin_gross_orbital_product(i) +! write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) +!enddo +!print*,'sum = ',accu +!accu = 0.d0 +!print*,'Angular momentum analysis' +!do i = 0, ao_l_max +! accu += spin_population_angular_momentum(i) +! print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) +!print*,'sum = ',accu +!enddo + +end + diff --git a/plugins/FOBOCI/modify_generators.irp.f b/plugins/FOBOCI/modify_generators.irp.f index c756f0c2..359b6405 100644 --- a/plugins/FOBOCI/modify_generators.irp.f +++ b/plugins/FOBOCI/modify_generators.irp.f @@ -6,6 +6,7 @@ subroutine set_generators_to_psi_det END_DOC N_det_generators = N_det integer :: i,k + print*,'N_det = ',N_det do i=1,N_det_generators do k=1,N_int psi_det_generators(k,1,i) = psi_det(k,1,i) diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 90b29faa..2e551dcd 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -213,6 +213,8 @@ subroutine new_approach deallocate(dressing_matrix_1h2p) deallocate(dressing_matrix_extra_1h_or_1p) enddo + + double precision, allocatable :: H_matrix_total(:,:) integer :: n_det_total n_det_total = N_det_generators_restart + n_good_det @@ -251,25 +253,79 @@ subroutine new_approach 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 + + ! Adding the correlation energy + logical :: orb_taken_good_det(mo_tot_num) + double precision :: phase + integer :: n_h,n_p,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + integer :: h1,h2,p1,p2,s1,s2 + logical, allocatable :: one_hole_or_one_p(:) + integer, allocatable :: holes_or_particle(:) + allocate(one_hole_or_one_p(n_good_det), holes_or_particle(n_good_det)) + orb_taken_good_det = .False. + do i = 1, n_good_det + n_h = number_of_holes(psi_good_det(1,1,i)) + n_p = number_of_particles(psi_good_det(1,1,i)) + call get_excitation(ref_bitmask,psi_good_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(n_h == 0 .and. n_p == 1)then + orb_taken_good_det(h1) = .True. + one_hole_or_one_p(i) = .True. + holes_or_particle(i) = h1 + endif + if(n_h == 1 .and. n_p == 0)then + orb_taken_good_det(p1) = .True. + one_hole_or_one_p(i) = .False. + holes_or_particle(i) = p1 + endif 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 + + do i = 1, N_det_generators_restart + ! Add the 2h2p, 2h1p and 1h2p correlation energy + H_matrix_total(i,i) += total_corr_e_2h2p + total_corr_e_2h1p + total_corr_e_1h2p + total_corr_e_1h1p_spin_flip + ! Substract the 2h1p part that have already been taken into account + do j = 1, n_inact_orb + iorb = list_inact(j) + if(.not.orb_taken_good_det(iorb))cycle + H_matrix_total(i,i) -= corr_energy_2h1p_per_orb_ab(iorb) - corr_energy_2h1p_per_orb_bb(iorb) - corr_energy_1h1p_spin_flip_per_orb(iorb) + enddo + ! Substract the 1h2p part that have already been taken into account + do j = 1, n_virt_orb + iorb = list_virt(j) + if(.not.orb_taken_good_det(iorb))cycle + H_matrix_total(i,i) -= corr_energy_1h2p_per_orb_ab(iorb) - corr_energy_1h2p_per_orb_aa(iorb) + enddo + enddo + + do i = 1, N_good_det + ! Repeat the 2h2p correlation energy + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_2h2p + ! Substract the part that can not be repeated + ! If it is a 1h + if(one_hole_or_one_p(i))then + ! 2h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h2p_per_orb_bb(holes_or_particle(i)) + ! You can repeat a certain part of the 1h2p correlation energy + ! that is everything except the part that involves the hole of the 1h + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_1h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_1h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_1h2p_per_orb_bb(holes_or_particle(i)) + + else + ! 2h2p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h2p_per_orb_aa(holes_or_particle(i)) + ! You can repeat a certain part of the 2h1p correlation energy + ! that is everything except the part that involves the hole of the 1p + ! 2h1p + H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h1p_per_orb_ab(holes_or_particle(i)) & + -corr_energy_2h1p_per_orb_aa(holes_or_particle(i)) + endif + enddo + 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 @@ -284,22 +340,222 @@ subroutine new_approach psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i) enddo enddo - norm = 0.d0 + + + double precision :: href + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + integer(bit_kind), allocatable :: psi_det_final(:,:,:) + double precision, allocatable :: psi_coef_final(:,:) + double precision :: norm + 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*,'' + print*,'' + print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1) + print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion 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) + print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1) enddo - print*,'norm = ',norm + + integer(bit_kind), allocatable :: psi_det_remaining_1h_or_1p(:,:,:) + integer(bit_kind), allocatable :: key_tmp(:,:) + integer :: n_det_remaining_1h_or_1p + integer :: ispin,i_ok + allocate(key_tmp(N_int,2),psi_det_remaining_1h_or_1p(N_int,2,n_inact_orb*n_act_orb+n_virt_orb*n_act_orb)) + logical :: is_already_present + logical, allocatable :: one_hole_or_one_p_bis(:) + integer, allocatable :: holes_or_particle_bis(:) + double precision,allocatable :: H_array(:) + allocate(one_hole_or_one_p_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb), holes_or_particle_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb)) + allocate(H_array(n_det_total)) + ! Dressing with the remaining 1h determinants + print*,'' + print*,'' + print*,'Dressing with the remaining 1h determinants' + n_det_remaining_1h_or_1p = 0 + do i = 1, n_inact_orb + iorb = list_inact(i) + if(orb_taken_good_det(iorb))cycle + do j = 1, n_act_orb + jorb = list_act(j) + ispin = 2 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,iorb,jorb,ispin,i_ok) + if(i_ok .ne.1)cycle + is_already_present = .False. + H_array = 0.d0 + call i_h_j(key_tmp,key_tmp,N_int,hij) + href = ref_bitmask_energy - hij + href = 1.d0/href + do k = 1, n_det_total + call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int) + if(degree == 0)then + is_already_present = .True. + exit + endif + enddo + if(is_already_present)cycle + n_det_remaining_1h_or_1p +=1 + one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .True. + holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb + do k = 1, N_int + psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1) + psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2) + enddo + ! do k = 1, n_det_total + ! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij) + ! H_array(k) = hij + ! enddo + ! do k = 1, n_det_total + ! do l = 1, n_det_total + ! H_matrix_total(k,l) += H_array(k) * H_array(l) * href + ! enddo + ! enddo + enddo + enddo + ! Dressing with the remaining 1p determinants + print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p + print*,'Dressing with the remaining 1p determinants' + do i = 1, n_virt_orb + iorb = list_virt(i) + if(orb_taken_good_det(iorb))cycle + do j = 1, n_act_orb + jorb = list_act(j) + ispin = 1 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,jorb,iorb,ispin,i_ok) + if(i_ok .ne.1)cycle + is_already_present = .False. + H_array = 0.d0 + call i_h_j(key_tmp,key_tmp,N_int,hij) + href = ref_bitmask_energy - hij + href = 1.d0/href + do k = 1, n_det_total + call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int) + if(degree == 0)then + is_already_present = .True. + exit + endif + enddo + if(is_already_present)cycle + n_det_remaining_1h_or_1p +=1 + one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .False. + holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb + do k = 1, N_int + psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1) + psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2) + enddo +! do k = 1, n_det_total +! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij) +! H_array(k) = hij +! enddo +! do k = 1, n_det_total +! do l = 1, n_det_total +! H_matrix_total(k,l) += H_array(k) * H_array(l) * href +! enddo +! enddo + enddo + enddo + print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p + deallocate(key_tmp,H_array) + + double precision, allocatable :: eigvalues_bis(:),eigvectors_bis(:,:),H_matrix_total_bis(:,:) + integer :: n_det_final + n_det_final = n_det_total + n_det_remaining_1h_or_1p + allocate(eigvalues_bis(n_det_final),eigvectors_bis(n_det_final,n_det_final),H_matrix_total_bis(n_det_final,n_det_final)) + print*,'passed the allocate, building the big matrix' + do i = 1, n_det_total + do j = 1, n_det_total + H_matrix_total_bis(i,j) = H_matrix_total(i,j) + enddo + enddo + do i = 1, n_det_remaining_1h_or_1p + do j = 1, n_det_remaining_1h_or_1p + call i_h_j(psi_det_remaining_1h_or_1p(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij) + H_matrix_total_bis(n_det_total+i,n_det_total+j) = hij + enddo + enddo + do i = 1, n_det_total + do j = 1, n_det_remaining_1h_or_1p + call i_h_j(psi_det_final(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij) + H_matrix_total_bis(i,n_det_total+j) = hij + H_matrix_total_bis(n_det_total+j,i) = hij + enddo + enddo + print*,'passed the matrix' + do i = 1, n_det_remaining_1h_or_1p + if(one_hole_or_one_p_bis(i))then + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h2p_per_orb_bb(holes_or_particle_bis(i)) + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_1h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_1h2p_per_orb_bb(holes_or_particle_bis(i)) + else + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h2p_per_orb_aa(holes_or_particle_bis(i)) + H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_2h1p_per_orb_ab(holes_or_particle_bis(i)) & + -corr_energy_2h1p_per_orb_aa(holes_or_particle_bis(i)) + + endif + enddo + do i = 2, n_det_final + do j = i+1, n_det_final + H_matrix_total_bis(i,j) = 0.d0 + H_matrix_total_bis(j,i) = 0.d0 + enddo + enddo + do i = 1, n_det_final + write(*,'(500(F10.5,X))')H_matrix_total_bis(i,:) + enddo + call lapack_diag(eigvalues_bis,eigvectors_bis,H_matrix_total_bis,n_det_final,n_det_final) + print*,'e_dressed = ',eigvalues_bis(1) + nuclear_repulsion + do i = 1, n_det_final + print*,'coef = ',eigvectors_bis(i,1),H_matrix_total_bis(i,i) - H_matrix_total_bis(1,1) + enddo + do j = 1, N_states + do i = 1, n_det_total + psi_coef_final(i,j) = eigvectors_bis(i,j) + norm += psi_coef_final(i,j)**2 + enddo + norm = 1.d0/dsqrt(norm) + do i = 1, n_det_total + psi_coef_final(i,j) = psi_coef_final(i,j) * norm + enddo + enddo + + + deallocate(eigvalues_bis,eigvectors_bis,H_matrix_total_bis) + + +!print*,'H matrix to diagonalize' +!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 +!call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total) +!print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1) +!print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion +!do i = 1, n_det_total +! print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1) +!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 +!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 + + 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 diff --git a/plugins/FOBOCI/new_new_approach.irp.f b/plugins/FOBOCI/new_new_approach.irp.f new file mode 100644 index 00000000..b904a5b3 --- /dev/null +++ b/plugins/FOBOCI/new_new_approach.irp.f @@ -0,0 +1,132 @@ +program test_new_new + implicit none + read_wf = .True. + touch read_wf + call test +end + + +subroutine test + implicit none + integer :: i,j,k,l + call diagonalize_CI + call set_generators_to_psi_det + print*,'Initial coefficients' + do i = 1, N_det + print*,'' + call debug_det(psi_det(1,1,i),N_int) + print*,'psi_coef = ',psi_coef(i,1) + print*,'' + enddo + double precision, allocatable :: dressing_matrix(:,:) + double precision :: hij + double precision :: phase + integer :: n_h,n_p,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + integer :: h1,h2,p1,p2,s1,s2 + allocate(dressing_matrix(N_det_generators,N_det_generators)) + do i = 1, N_det_generators + do j = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij) + dressing_matrix(i,j) = hij + enddo + enddo + href = dressing_matrix(1,1) + print*,'Diagonal part of the dressing' + do i = 1, N_det_generators + print*,'delta e = ',dressing_matrix(i,i) - href + enddo + call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) + double precision :: href + print*,'' + ! One considers that the following excitation classes are not repeatable on the 1h and 1p determinants : + ! + 1h1p spin flip + ! + 2h1p + ! + 1h2p + ! But the 2h2p are correctly taken into account +!dressing_matrix(1,1) += total_corr_e_1h2p + total_corr_e_2h1p + total_corr_e_1h1p_spin_flip +!do i = 1, N_det_generators +! dressing_matrix(i,i) += total_corr_e_2h2p +! 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 +! +! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! print*,'' +! print*,' 1h det ' +! print*,'' +! call debug_det(psi_det_generators(1,1,i),N_int) +! print*,'h1,p1 = ',h1,p1 +! print*,'total_corr_e_2h2p ',total_corr_e_2h2p +! print*,'corr_energy_2h2p_per_orb_ab(h1)',corr_energy_2h2p_per_orb_ab(h1) +! print*,'corr_energy_2h2p_per_orb_bb(h1)',corr_energy_2h2p_per_orb_bb(h1) +! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1) +! dressing_matrix(1,1) += -corr_energy_2h1p_per_orb_aa(h1) - corr_energy_2h1p_per_orb_ab(h1) -corr_energy_2h1p_per_orb_bb(h1) & +! -corr_energy_1h1p_spin_flip_per_orb(h1) +! endif +! if(n_h == 0 .and. n_p ==1)then +! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! print*,'' +! print*,' 1p det ' +! print*,'' +! call debug_det(psi_det_generators(1,1,i),N_int) +! print*,'h1,p1 = ',h1,p1 +! print*,'total_corr_e_2h2p ',total_corr_e_2h2p +! print*,'corr_energy_2h2p_per_orb_ab(p1)',corr_energy_2h2p_per_orb_ab(p1) +! print*,'corr_energy_2h2p_per_orb_aa(p1)',corr_energy_2h2p_per_orb_aa(p1) +! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(p1) - corr_energy_2h2p_per_orb_aa(p1) +! dressing_matrix(1,1) += -corr_energy_1h2p_per_orb_aa(p1) - corr_energy_1h2p_per_orb_ab(p1) -corr_energy_1h2p_per_orb_bb(p1) +! endif +!enddo +!href = dressing_matrix(1,1) +!print*,'Diagonal part of the dressing' +!do i = 1, N_det_generators +! print*,'delta e = ',dressing_matrix(i,i) - href +!enddo + call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) + print*,'After dressing matrix' + print*,'' + print*,'' + do i = 1, N_det + print*,'psi_coef = ',psi_coef(i,1) + enddo +!print*,'' +!print*,'' +!print*,'Canceling the dressing part of the interaction between 1h and 1p' +!do i = 2, N_det_generators +! do j = i+1, N_det_generators +! call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij) +! dressing_matrix(i,j) = hij +! dressing_matrix(j,i) = hij +! enddo +!enddo +!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) +!print*,'' +!print*,'' +!do i = 1, N_det +! print*,'psi_coef = ',psi_coef(i,1) +!enddo +!print*,'' +!print*,'' +!print*,'Canceling the interaction between 1h and 1p' + +!print*,'' +!print*,'' +!do i = 2, N_det_generators +! do j = i+1, N_det_generators +! dressing_matrix(i,j) = 0.d0 +! dressing_matrix(j,i) = 0.d0 +! enddo +!enddo +!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) +!do i = 1, N_det +! print*,'psi_coef = ',psi_coef(i,1) +!enddo + call save_natural_mos + deallocate(dressing_matrix) + + +end diff --git a/plugins/FOBOCI/routines_dressing.irp.f b/plugins/FOBOCI/routines_dressing.irp.f index b0edd949..125143da 100644 --- a/plugins/FOBOCI/routines_dressing.irp.f +++ b/plugins/FOBOCI/routines_dressing.irp.f @@ -55,15 +55,11 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det i_pert = 0 endif do j = 1, ndet_generators_input - if(dabs(H_array(j)*lambda_i).gt.0.5d0)then + if(dabs(H_array(j)*lambda_i).gt.0.1d0)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 @@ -79,8 +75,52 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det enddo enddo enddo + href = dressing_matrix(1,1) + print*,'Diagonal part of the dressing' + do i = 1, ndet_generators_input + print*,'delta e = ',dressing_matrix(i,i) - href + enddo !print*,'i_pert_count = ',i_pert_count end + +subroutine update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,H_jj_in) + 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(in) :: H_jj_in(N_det) + double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input) + integer :: i,j,n_det_ref_tmp,degree + double precision :: href + 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 + dressing_matrix(j,j) += H_jj_in(i) + n_det_ref_tmp +=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 + + href = dressing_matrix(1,1) + print*,'' + print*,'Update with the SC2 dressing' + print*,'' + print*,'Diagonal part of the dressing' + do i = 1, ndet_generators_input + print*,'delta e = ',dressing_matrix(i,i) - href + enddo +end + subroutine provide_matrix_dressing_for_extra_1h_or_1p(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 @@ -125,11 +165,14 @@ subroutine provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix,psi_det_re i_pert = 0 endif do j = 1, n_det_ref_input - if(dabs(H_array(j)*lambda_i).gt.0.3d0)then + if(dabs(H_array(j)*lambda_i).gt.0.5d0)then i_pert = 1 exit endif enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! i_pert = 0 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(i_pert==1)then lambda_i = f i_pert_count +=1 @@ -178,16 +221,17 @@ subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi accu += psi_coef_ref_input(j,1) * hka enddo lambda_i = psi_coef_outer_input(i,1)/accu - i_pert = 1 + i_pert = 0 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 + if(dabs(H_array(j)*lambda_i).gt.0.5d0)then i_pert = 1 exit endif enddo +! i_pert = 0 if(i_pert==1)then lambda_i = f i_pert_count +=1 @@ -275,19 +319,257 @@ subroutine give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1 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 - + if(n_det_2h1p + n_det_1h1p + n_det_extra_1h_or_1p + n_det_generators .ne. N_det)then + print*,'PB !!!!' + print*,'You have forgotten something in your generators ... ' + stop + endif end +subroutine give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p) + use bitmasks + implicit none + integer, intent(out) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + n_det_ref_1h_1p = 0 + n_det_2h1p = 0 + n_det_1h1p = 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_1h_1p +=1 + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p +=1 + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p +=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, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + +end + +subroutine give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p) + use bitmasks + implicit none + integer, intent(out) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + integer :: i + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + n_det_ref_1h_1p = 0 + n_det_1h2p = 0 + n_det_1h1p = 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_1h_1p +=1 + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p +=1 + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p +=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 1h, 1p, 1h1p or 1h2p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + +end + +subroutine give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p + integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p) + integer(bit_kind), intent(out) :: psi_det_2h1p(N_int,2,n_det_2h1p) + integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p) + double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states) + double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states) + integer :: n_det_ref_1h_1p_tmp,n_det_2h1p_tmp,n_det_1h1p_tmp + integer :: i,j + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + integer, allocatable :: index_ref_1h_1p(:) + integer, allocatable :: index_2h1p(:) + integer, allocatable :: index_1h1p(:) + allocate(index_ref_1h_1p(n_det)) + allocate(index_2h1p(n_det)) + allocate(index_1h1p(n_det)) + n_det_ref_1h_1p_tmp = 0 + n_det_2h1p_tmp = 0 + n_det_1h1p_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 == 0 .and. n_p == 0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else 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 + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + do i = 1, n_det_2h1p + do j = 1, N_int + psi_det_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i)) + psi_det_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 + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_det_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_ref_1h_1p + do j = 1, N_int + psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i)) + psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i)) + enddo + do j = 1, N_states + psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j) + enddo + enddo + +end + +subroutine give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,& + psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p + integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p) + integer(bit_kind), intent(out) :: psi_det_1h2p(N_int,2,n_det_1h2p) + integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p) + double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states) + double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p,N_states) + double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states) + integer :: n_det_ref_1h_1p_tmp,n_det_1h2p_tmp,n_det_1h1p_tmp + integer :: i,j + integer :: n_det_ref_restart_tmp,n_det_1h + integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det + integer, allocatable :: index_ref_1h_1p(:) + integer, allocatable :: index_1h2p(:) + integer, allocatable :: index_1h1p(:) + allocate(index_ref_1h_1p(n_det)) + allocate(index_1h2p(n_det)) + allocate(index_1h1p(n_det)) + n_det_ref_1h_1p_tmp = 0 + n_det_1h2p_tmp = 0 + n_det_1h1p_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 == 0 .and. n_p == 0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_ref_1h_1p_tmp +=1 + index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i + else 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 + else + print*,'PB !!!!' + print*,'You have something else than a 1h, 1p, 1h1p or 1h2p' + print*,'n_h,n_p = ',n_h,n_p + call debug_det(psi_det(1,1,i),N_int) + stop + endif + enddo + do i = 1, n_det_1h2p + do j = 1, N_int + psi_det_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i)) + psi_det_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 + + do i = 1, n_det_1h1p + do j = 1, N_int + psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i)) + psi_det_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_ref_1h_1p + do j = 1, N_int + psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i)) + psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i)) + enddo + do j = 1, N_states + psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j) + enddo + enddo + +end + + + subroutine give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p) use bitmasks implicit none @@ -353,7 +635,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d 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,n_det_extra_1h_or_1p_tmp - integer :: n_det_extra_1h_tmp + integer :: n_det_1h_tmp integer, allocatable :: index_generator(:) integer, allocatable :: index_1h1p(:) integer, allocatable :: index_2h1p(:) @@ -370,7 +652,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d n_det_1h1p_tmp = 0 n_det_2h1p_tmp = 0 n_det_extra_1h_or_1p_tmp = 0 - n_det_extra_1h_tmp = 0 + n_det_1h_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)) @@ -385,7 +667,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i else 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_det_extra_1h_tmp +=1 + n_det_1h_tmp +=1 else n_det_extra_1h_or_1p_tmp +=1 index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 696011a9..4def99e2 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -332,20 +332,20 @@ subroutine save_osoci_natural_mos 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-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 @@ -389,14 +389,14 @@ subroutine save_osoci_natural_mos jorb = list_inact(j) if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then print*,'INACTIVE ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(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)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo enddo @@ -410,8 +410,9 @@ subroutine save_osoci_natural_mos enddo label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) - soft_touch mo_coef +!soft_touch mo_coef deallocate(tmp,occ) @@ -520,14 +521,14 @@ subroutine set_osoci_natural_mos jorb = list_inact(j) if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then print*,'INACTIVE ' - print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb)) + print*,'DM ',iorb,jorb,(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)) + print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo enddo @@ -602,15 +603,7 @@ 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 + call print_mulliken_sd + call print_hcc end diff --git a/plugins/Generators_restart/generators.irp.f b/plugins/Generators_restart/generators.irp.f index 0a82e6f9..17854330 100644 --- a/plugins/Generators_restart/generators.irp.f +++ b/plugins/Generators_restart/generators.irp.f @@ -1,5 +1,5 @@ use bitmasks - + BEGIN_PROVIDER [ integer, N_det_generators ] implicit none BEGIN_DOC @@ -8,17 +8,18 @@ BEGIN_PROVIDER [ integer, N_det_generators ] integer :: i integer, save :: ifirst = 0 double precision :: norm - read_wf = .True. if(ifirst == 0)then - N_det_generators = N_det + call ezfio_get_determinants_n_det(N_det_generators) ifirst = 1 + else + print*,'PB in generators restart !!!' endif call write_int(output_determinants,N_det_generators,'Number of generators') END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ] implicit none BEGIN_DOC ! read wf @@ -26,17 +27,20 @@ END_PROVIDER END_DOC integer :: i, k integer, save :: ifirst = 0 + double precision, allocatable :: psi_coef_read(:,:) if(ifirst == 0)then - 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 + call read_dets(psi_det_generators,N_int,N_det_generators) + allocate (psi_coef_read(N_det_generators,N_states)) + call ezfio_get_determinants_psi_coef(psi_coef_read) do k = 1, N_states - psi_coef_generators(i,k) = psi_coef(i,k) + do i = 1, N_det_generators + psi_coef_generators(i,k) = psi_coef_read(i,k) + enddo enddo - enddo ifirst = 1 + deallocate(psi_coef_read) + else + print*,'PB in generators restart !!!' endif END_PROVIDER diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index d77c91c5..5c626db9 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,9 @@ subroutine damping_SCF write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' write(output_hartree_fock,*) - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + if(.not.no_oa_or_av_opt)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + endif call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') call ezfio_set_hartree_fock_energy(E_min) diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 72d03808..0086c67e 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -126,6 +126,8 @@ subroutine pt2_moller_plesset ($arguments) delta_e = (Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)) + & (Fock_matrix_diag_mo(h2) - Fock_matrix_diag_mo(p2)) delta_e = 1.d0/delta_e +! print*,'h1,p1',h1,p1 +! print*,'h2,p2',h2,p2 else if (degree == 1) then call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1) diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f index c1d88d2c..e31b3ba4 100644 --- a/plugins/Properties/hyperfine_constants.irp.f +++ b/plugins/Properties/hyperfine_constants.irp.f @@ -133,3 +133,16 @@ END_PROVIDER enddo END_PROVIDER + + +subroutine print_hcc + implicit none + double precision :: accu + integer :: i,j + print*,'Z AU GAUSS MHZ cm^-1' + do i = 1, nucl_num + write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) + enddo + +end + diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f index d56c9a44..cc0a2f8e 100644 --- a/plugins/Properties/mulliken.irp.f +++ b/plugins/Properties/mulliken.irp.f @@ -105,3 +105,34 @@ END_PROVIDER enddo END_PROVIDER + + +subroutine print_mulliken_sd + implicit none + double precision :: accu + integer :: i + integer :: j + print*,'Mulliken spin densities' + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),mulliken_spin_densities(i) + accu += mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu + print*,'AO SPIN POPULATIONS' + accu = 0.d0 + do i = 1, ao_num + accu += spin_gross_orbital_product(i) + write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) + enddo + print*,'sum = ',accu + accu = 0.d0 + print*,'Angular momentum analysis' + do i = 0, ao_l_max + accu += spin_population_angular_momentum(i) + print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) + print*,'sum = ',accu + enddo + +end + diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f index f0091e1e..45bca5e6 100644 --- a/plugins/Properties/print_hcc.irp.f +++ b/plugins/Properties/print_hcc.irp.f @@ -1,17 +1,6 @@ -program print_hcc +program print_hcc_main implicit none read_wf = .True. touch read_wf - call test + call print_hcc end -subroutine test - implicit none - double precision :: accu - integer :: i,j - print*,'Z AU GAUSS MHZ cm^-1' - do i = 1, nucl_num - write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) - enddo - -end - diff --git a/plugins/Properties/print_mulliken.irp.f b/plugins/Properties/print_mulliken.irp.f index 100c8556..d4be534a 100644 --- a/plugins/Properties/print_mulliken.irp.f +++ b/plugins/Properties/print_mulliken.irp.f @@ -2,34 +2,5 @@ program print_mulliken implicit none read_wf = .True. touch read_wf - print*,'Mulliken spin densities' - - call test + call print_mulliken_sd end -subroutine test - double precision :: accu - integer :: i - integer :: j - accu= 0.d0 - do i = 1, nucl_num - print*,i,nucl_charge(i),mulliken_spin_densities(i) - accu += mulliken_spin_densities(i) - enddo - print*,'Sum of Mulliken SD = ',accu - print*,'AO SPIN POPULATIONS' - accu = 0.d0 - do i = 1, ao_num - accu += spin_gross_orbital_product(i) - write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) - enddo - print*,'sum = ',accu - accu = 0.d0 - print*,'Angular momentum analysis' - do i = 0, ao_l_max - accu += spin_population_angular_momentum(i) - print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) - print*,'sum = ',accu - enddo - -end - diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index f9be4fa6..84150c1b 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -24,13 +24,18 @@ skip init_main filter_integrals filter2p -filter2h2p +filter2h2p_double +filter2h2p_single filter1h filter1p only_2p_single only_2p_double filter_only_1h1p_single filter_only_1h1p_double +filter_only_1h2p_single +filter_only_1h2p_double +filter_only_2h2p_single +filter_only_2h2p_double filterhole filterparticle do_double_excitations @@ -194,6 +199,27 @@ class H_apply(object): if (is_a_1h1p(key).eqv..False.) cycle """ + def filter_only_2h2p(self): + self["filter_only_2h2p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_two_holes_two_particles(hole).eqv..False.) cycle + """ + self["filter_only_1h1p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_two_holes_two_particles(key).eqv..False.) cycle + """ + + + def filter_only_1h2p(self): + self["filter_only_1h2p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h2p(hole).eqv..False.) cycle + """ + self["filter_only_1h2p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h2p(key).eqv..False.) cycle + """ + def unset_skip(self): self["skip"] = """ @@ -201,9 +227,12 @@ class H_apply(object): def set_filter_2h_2p(self): - self["filter2h2p"] = """ + self["filter2h2p_double"] = """ if (is_a_two_holes_two_particles(key)) cycle """ + self["filter2h2p_single"] = """ + if (is_a_two_holes_two_particles(hole)) cycle + """ def set_perturbation(self,pert): diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 4441fb22..4984d9a8 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -212,6 +212,12 @@ logical function is_a_two_holes_two_particles(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) integer :: i,i_diff + integer :: number_of_holes, number_of_particles + is_a_two_holes_two_particles = .False. + if(number_of_holes(key_in) == 2 .and. number_of_particles(key_in) == 2)then + is_a_two_holes_two_particles = .True. + return + endif i_diff = 0 if(N_int == 1)then i_diff = i_diff & @@ -456,6 +462,17 @@ logical function is_a_1h1p(key_in) end +logical function is_a_1h2p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h2p = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.2)then + is_a_1h2p = .True. + endif + +end + logical function is_a_1h(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index d6d8fcb0..3f8299f6 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -98,9 +98,40 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ] END_PROVIDER +BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ] + implicit none + BEGIN_DOC + ! Number of bitmasks for generators + END_DOC + logical :: exists + PROVIDE ezfio_filename + + call ezfio_has_bitmasks_N_mask_gen(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_generators_bitmask_restart = 1 + endif + ASSERT (N_generators_bitmask_restart > 0) + +END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask) ] + + +BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ] implicit none BEGIN_DOC ! Bitmasks for generator determinants. @@ -258,7 +289,7 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] call ezfio_get_bitmasks_cas(cas_bitmask) print*,'---------------------' else - if(N_generators_bitmask == 1)then + if(N_generators_bitmask_restart == 1)then do i=1,N_cas_bitmask cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) enddo @@ -302,7 +333,7 @@ END_PROVIDER n_inact_orb = 0 n_virt_orb = 0 - if(N_generators_bitmask == 1)then + if(N_generators_bitmask_restart == 1)then do j = 1, N_int inact_bitmask(j,1) = xor(generators_bitmask_restart(j,1,1,1),cas_bitmask(j,1,1)) inact_bitmask(j,2) = xor(generators_bitmask_restart(j,2,1,1),cas_bitmask(j,2,1)) @@ -315,15 +346,15 @@ END_PROVIDER i_hole = 1 i_gen = 1 do i = 1, N_int - inact_bitmask(i,1) = generators_bitmask(i,1,i_hole,i_gen) - inact_bitmask(i,2) = generators_bitmask(i,2,i_hole,i_gen) + inact_bitmask(i,1) = generators_bitmask_restart(i,1,i_hole,i_gen) + inact_bitmask(i,2) = generators_bitmask_restart(i,2,i_hole,i_gen) n_inact_orb += popcnt(inact_bitmask(i,1)) enddo i_part = 2 i_gen = 3 do i = 1, N_int - virt_bitmask(i,1) = generators_bitmask(i,1,i_part,i_gen) - virt_bitmask(i,2) = generators_bitmask(i,2,i_part,i_gen) + virt_bitmask(i,1) = generators_bitmask_restart(i,1,i_part,i_gen) + virt_bitmask(i,2) = generators_bitmask_restart(i,2,i_part,i_gen) n_virt_orb += popcnt(virt_bitmask(i,1)) enddo endif diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 86780430..fdf2ec80 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -166,6 +166,7 @@ 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_1h2p logical :: is_a_1h logical :: is_a_1p logical :: is_a_2p @@ -301,8 +302,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl 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) - $filter2h2p + $filter2h2p_double $filter_only_1h1p_double + $filter_only_1h2p_double + $filter_only_2h2p_double $only_2p_double key_idx += 1 do k=1,N_int @@ -353,8 +356,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl 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) - $filter2h2p + $filter2h2p_double $filter_only_1h1p_double + $filter_only_1h2p_double + $filter_only_2h2p_double $only_2p_double key_idx += 1 do k=1,N_int @@ -427,6 +432,7 @@ 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_1h2p logical :: is_a_1h logical :: is_a_1p logical :: is_a_2p @@ -505,8 +511,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato $filter1h $filter1p $filter2p - $filter2h2p + $filter2h2p_single $filter_only_1h1p_single + $filter_only_1h2p_single + $filter_only_2h2p_single key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) @@ -566,7 +574,6 @@ subroutine $subroutine($params_main) iproc = 0 allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) ) do i_generator=1,nmax - progress_bar(1) = i_generator if (abort_here) then @@ -600,6 +607,16 @@ subroutine $subroutine($params_main) not(psi_det_generators(k,ispin,i_generator)) ) enddo enddo +! print*,'generator in ' +! call debug_det(psi_det_generators(1,1,i_generator),N_int) +! print*,'hole 1' +! call debug_det(mask(1,1,d_hole1),N_int) +! print*,'hole 2' +! call debug_det(mask(1,1,d_hole2),N_int) +! print*,'part 1' +! call debug_det(mask(1,1,d_part1),N_int) +! print*,'part 2' +! call debug_det(mask(1,1,d_part2),N_int) if($do_double_excitations)then call $subroutine_diexc(psi_det_generators(1,1,i_generator), & psi_det_generators(1,1,1), & diff --git a/src/Determinants/SC2.irp.f b/src/Determinants/SC2.irp.f index 440b2870..4f613abc 100644 --- a/src/Determinants/SC2.irp.f +++ b/src/Determinants/SC2.irp.f @@ -1,4 +1,4 @@ -subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) +subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence) use bitmasks implicit none BEGIN_DOC @@ -21,6 +21,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(inout) :: u_in(dim_in,N_st) double precision, intent(out) :: energies(N_st) + double precision, intent(out) :: diag_H_elements(dim_in) double precision, intent(in) :: convergence ASSERT (N_st > 0) ASSERT (sze > 0) @@ -200,6 +201,9 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) converged = dabs(e_corr_double - e_corr_double_before) < convergence converged = converged .or. abort_here if (converged) then + do i = 1, dim_in + diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i) + enddo exit endif e_corr_double_before = e_corr_double diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 63ed7a92..f453a8a3 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -58,7 +58,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ] else psi_det_size = 1 endif - psi_det_size = max(psi_det_size,10000) + psi_det_size = max(psi_det_size,100000) call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays') END_PROVIDER diff --git a/src/Determinants/diagonalize_CI_SC2.irp.f b/src/Determinants/diagonalize_CI_SC2.irp.f index 97161ad3..498792d9 100644 --- a/src/Determinants/diagonalize_CI_SC2.irp.f +++ b/src/Determinants/diagonalize_CI_SC2.irp.f @@ -23,8 +23,10 @@ END_PROVIDER threshold_convergence_SC2 = 1.d-10 END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ] implicit none BEGIN_DOC ! Eigenvectors/values of the CI matrix @@ -39,7 +41,8 @@ END_PROVIDER enddo call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & - size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) +! size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) END_PROVIDER subroutine diagonalize_CI_SC2 @@ -54,5 +57,6 @@ subroutine diagonalize_CI_SC2 psi_coef(i,j) = CI_SC2_eigenvectors(i,j) enddo enddo - SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors + SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2 +! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors end diff --git a/src/Determinants/save_natorb.irp.f b/src/Determinants/save_natorb.irp.f index e56f9821..674ba32e 100644 --- a/src/Determinants/save_natorb.irp.f +++ b/src/Determinants/save_natorb.irp.f @@ -2,5 +2,6 @@ program save_natorb read_wf = .True. touch read_wf call save_natural_mos + call save_ref_determinant end diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 84b08715..dc35f278 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -230,7 +230,6 @@ subroutine clear_ao_map end - !! MO Map !! ====== diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 0ff14168..000fcb4d 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -72,7 +72,7 @@ subroutine add_integrals_to_map(mask_ijkl) integer :: i2,i3,i4 double precision,parameter :: thr_coef = 1.d-10 - PROVIDE ao_bielec_integrals_in_map + PROVIDE ao_bielec_integrals_in_map mo_coef !Get list of MOs for i,j,k and l !------------------------------- @@ -341,7 +341,7 @@ end double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) if (.not.do_direct_integrals) then - PROVIDE ao_bielec_integrals_in_map + PROVIDE ao_bielec_integrals_in_map mo_coef endif mo_bielec_integral_jj_from_ao = 0.d0 @@ -513,4 +513,13 @@ subroutine clear_mo_map call map_deinit(mo_integrals_map) FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + + +end + +subroutine provide_all_mo_integrals + implicit none + provide mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti + provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map + end diff --git a/src/Integrals_Monoelec/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f index 714222ec..5bae9868 100644 --- a/src/Integrals_Monoelec/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -5,6 +5,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to ! array of the mono electronic hamiltonian on the MOs basis ! : sum of the kinetic and nuclear electronic potential END_DOC + print*,'Providing the mono electronic integrals' do j = 1, mo_tot_num do i = 1, mo_tot_num mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)