diff --git a/plugins/FOBOCI/EZFIO.cfg b/plugins/FOBOCI/EZFIO.cfg index 88189608..9b9f7d71 100644 --- a/plugins/FOBOCI/EZFIO.cfg +++ b/plugins/FOBOCI/EZFIO.cfg @@ -19,10 +19,15 @@ default: 0.00001 [do_it_perturbative] type: logical -doc: if true, you do the FOBOCI calculation perturbatively +doc: if true, when a given 1h or 1p determinant is not selected because of its perturbation estimate, then if its coefficient is lower than threshold_perturbative, it is acounted in the FOBOCI differential density matrices interface: ezfio,provider,ocaml default: .False. +[threshold_perturbative] +type: double precision +doc: when do_it_perturbative is True, threshold_perturbative select if a given determinant ia selected or not for beign taken into account in the FOBO-SCF treatment. In practive, if the coefficient is larger then threshold_perturbative it means that it not selected as the perturbation should not be too importan. A value of 0.01 is in general OK. +interface: ezfio,provider,ocaml +default: 0.001 [speed_up_convergence_foboscf] type: logical @@ -49,3 +54,9 @@ doc: if true, you do all 2p type excitation on the LMCT interface: ezfio,provider,ocaml default: .True. +[selected_fobo_ci] +type: logical +doc: if true, for each CI step you will run a CIPSI calculation that stops at pt2_max +interface: ezfio,provider,ocaml +default: .False. + diff --git a/plugins/FOBOCI/SC2_1h1p.irp.f b/plugins/FOBOCI/SC2_1h1p.irp.f index d347c6e5..6f6156f4 100644 --- a/plugins/FOBOCI/SC2_1h1p.irp.f +++ b/plugins/FOBOCI/SC2_1h1p.irp.f @@ -158,6 +158,7 @@ subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,conve ! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma} ! diag_H_elements(index_hf) += total_corr_e_2h2p + return c_ref = c_ref * c_ref print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf) do i = 1, n_singles @@ -186,6 +187,186 @@ subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,conve end + +subroutine dressing_1h1p_by_2h2p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: diag_H_elements(dim_in) + double precision, intent(in) :: convergence + + integer :: i,j,k,l + integer :: r,s,i0,j0,r0,s0 + integer :: n_singles + integer :: index_singles(sze),hole_particles_singles(sze,3) + integer :: n_doubles + integer :: index_doubles(sze),hole_particles_doubles(sze,2) + integer :: index_hf + double precision :: e_corr_singles(mo_tot_num,2) + double precision :: e_corr_doubles(mo_tot_num) + double precision :: e_corr_singles_total(2) + double precision :: e_corr_doubles_1h1p + + integer :: exc(0:2,2,2),degree + integer :: h1,h2,p1,p2,s1,s2 + integer :: other_spin(2) + double precision :: phase + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i_ok + double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz + double precision :: hij,c_ref,contrib + integer :: iorb + + other_spin(1) = 2 + other_spin(2) = 1 + + n_singles = 0 + n_doubles = 0 + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij) + diag_H_elements(i) = hij + if(degree == 0)then + index_hf = i + else if (degree == 1)then + n_singles +=1 + index_singles(n_singles) = i + ! h1 = inactive orbital of the hole + hole_particles_singles(n_singles,1) = h1 + ! p1 = virtual orbital of the particle + hole_particles_singles(n_singles,2) = p1 + ! s1 = spin of the electron excited + hole_particles_singles(n_singles,3) = s1 + else if (degree == 2)then + n_doubles +=1 + index_doubles(n_doubles) = i + ! h1 = inactive orbital of the hole (beta of course) + hole_particles_doubles(n_doubles,1) = h1 + ! p1 = virtual orbital of the particle (alpha of course) + hole_particles_doubles(n_doubles,2) = p2 + else + print*,'PB !! found out other thing than a single or double' + print*,'stopping ..' + stop + endif + enddo + double precision :: delta_e + double precision :: coef_ijrs + diag_H_elements = 0.d0 + do i0 = 1, n_core_inact_orb + i= list_core_inact(i0) + do j0 = i0+1, n_core_inact_orb + j = list_core_inact(j0) + print*, i,j + do r0 = 1, n_virt_orb + r = list_virt(r0) + do s0 = r0+1, n_virt_orb + s = list_virt(s0) + !!! alpha (i-->r) / beta (j-->s) + s1 = 1 + s2 = 2 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,i,r,s1,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call do_mono_excitation(key_tmp,j,s,s2,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call i_H_j(ref_bitmask, key_tmp, N_int,hij) + delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s) + coef_ijrs = hij/delta_e + do k = 1, n_singles + l = index_singles(k) + call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij) + diag_H_elements(l) += coef_ijrs * hij + enddo + !if(i>j.and.r>s)then + !! alpha (i-->r) / alpha (j-->s) + s1 = 1 + s2 = 1 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,i,r,s1,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call do_mono_excitation(key_tmp,j,s,s2,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call i_H_j(ref_bitmask, key_tmp, N_int,hij) + delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s) + coef_ijrs = hij/delta_e + do k = 1, n_singles + l = index_singles(k) + call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij) + diag_H_elements(l) += coef_ijrs * hij + enddo + !! beta (i-->r) / beta (j-->s) + s1 = 2 + s2 = 2 + key_tmp = ref_bitmask + call do_mono_excitation(key_tmp,i,r,s1,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call do_mono_excitation(key_tmp,j,s,s2,i_ok) + if(i_ok .ne.1)then + print*, 'pb !!' + stop + endif + call i_H_j(ref_bitmask, key_tmp, N_int,hij) + delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s) + coef_ijrs = hij/delta_e + do k = 1, n_singles + l = index_singles(k) + call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij) + diag_H_elements(l) += coef_ijrs * hij + enddo + !endif + enddo + enddo + enddo + enddo + c_ref = 1.d0/u_in(index_hf,1) + do k = 1, n_singles + l = index_singles(k) + diag_H_elements(0) -= diag_H_elements(l) + enddo +! do k = 1, n_doubles +! l = index_doubles(k) +! diag_H_elements(0) += diag_H_elements(l) +! enddo + + +end + + subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence) use bitmasks implicit none @@ -478,11 +659,13 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c 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 :: extra_diag_H_elements(dim_in) double precision, intent(in) :: convergence integer :: i,j,iter do iter = 1, 1 - call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) - if(sze<=N_det_max_jacobi)then +! call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) + call dressing_1h1p_by_2h2p(dets_in,u_in,extra_diag_H_elements,dim_in,sze,N_st,Nint,convergence) +! if(sze<=N_det_max_jacobi)then double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) do j=1,sze @@ -490,9 +673,14 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) enddo enddo - do i = 1,sze - H_matrix_tmp(i,i) = diag_H_elements(i) + H_matrix_tmp(1,1) += extra_diag_H_elements(1) + do i = 2,sze + H_matrix_tmp(1,i) += extra_diag_H_elements(i) + H_matrix_tmp(i,1) += extra_diag_H_elements(i) enddo + !do i = 1,sze + ! H_matrix_tmp(i,i) = diag_H_elements(i) + !enddo call lapack_diag(eigenvalues,eigenvectors, & H_matrix_tmp,size(H_matrix_all_dets,1),sze) do j=1,min(N_states_diag,sze) @@ -502,9 +690,9 @@ subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,c energies(j) = eigenvalues(j) enddo deallocate (H_matrix_tmp, eigenvalues, eigenvectors) - else - call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants) - endif +! else +! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants) +! endif print*,'E = ',energies(1) + nuclear_repulsion enddo diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f index 0594e56e..2968ab90 100644 --- a/plugins/FOBOCI/all_singles.irp.f +++ b/plugins/FOBOCI/all_singles.irp.f @@ -1,13 +1,25 @@ -subroutine all_single +subroutine all_single(e_pt2) implicit none + double precision, intent(in) :: e_pt2 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 + if(.not.selected_fobo_ci)then + selection_criterion = 0.d0 + soft_touch selection_criterion + else + selection_criterion = 0.1d0 + selection_criterion_factor = 0.01d0 + selection_criterion_min = selection_criterion + soft_touch selection_criterion + endif + print*, 'e_pt2 = ',e_pt2 + pt2_max = 0.15d0 * e_pt2 + soft_touch pt2_max + print*, 'pt2_max = ',pt2_max threshold_davidson = 1.d-9 soft_touch threshold_davidson davidson_criterion i = 0 @@ -17,6 +29,8 @@ subroutine all_single print*,'pt2_max = ',pt2_max print*,'N_det_generators = ',N_det_generators pt2=-1.d0 + print*, 'ref_bitmask_energy =',ref_bitmask_energy + print*, 'CI_expectation_value =',CI_expectation_value(1) E_before = ref_bitmask_energy print*,'Initial Step ' @@ -29,7 +43,7 @@ subroutine all_single print*,'S^2 = ',CI_eigenvectors_s2(i) enddo n_det_max = 100000 - do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > dabs(pt2_max)) i += 1 print*,'-----------------------' print*,'i = ',i @@ -39,6 +53,8 @@ subroutine all_single print*,'E = ',CI_energy(1) print*,'pt2 = ',pt2(1) print*,'E+PT2 = ',E_before + pt2(1) + print*,'pt2_max = ',pt2_max + print*, maxval(abs(pt2(1:N_st))) > dabs(pt2_max) if(N_states_diag.gt.1)then print*,'Variational Energy difference' do i = 2, N_st @@ -53,7 +69,6 @@ subroutine all_single 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-8 ! soft_touch threshold_davidson davidson_criterion diff --git a/plugins/FOBOCI/create_1h_or_1p.irp.f b/plugins/FOBOCI/create_1h_or_1p.irp.f index 140ed504..41ec7b6c 100644 --- a/plugins/FOBOCI/create_1h_or_1p.irp.f +++ b/plugins/FOBOCI/create_1h_or_1p.irp.f @@ -68,7 +68,9 @@ subroutine create_restart_and_1h(i_hole) SOFT_TOUCH N_det psi_det psi_coef logical :: found_duplicates + if(n_act_orb.gt.1)then call remove_duplicates_in_psi_det(found_duplicates) + endif end subroutine create_restart_and_1p(i_particle) @@ -213,6 +215,8 @@ subroutine create_restart_1h_1p(i_hole,i_part) SOFT_TOUCH N_det psi_det psi_coef logical :: found_duplicates + if(n_act_orb.gt.1)then call remove_duplicates_in_psi_det(found_duplicates) + endif end diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index 74759362..e6521c76 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -72,20 +72,21 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen end -subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) +subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) use bitmasks implicit none double precision, intent(in) :: threshold - logical, intent(out) :: is_ok,exit_loop + double precision, intent(out):: e_pt2 + logical, intent(out) :: is_ok,exit_loop,is_ok_perturbative logical, intent(in) :: verbose integer :: l,k,m double precision,allocatable :: dressed_H_matrix(:,:) - double precision,allocatable :: psi_coef_diagonalized_tmp(:,:) + double precision, allocatable :: psi_coef_diagonalized_tmp(:,:) integer(bit_kind), allocatable :: psi_det_generators_input(:,:,:) + double precision :: hij - allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators)) - allocate(psi_coef_diagonalized_tmp(N_det_generators,N_states)) + allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators),psi_coef_diagonalized_tmp(N_det_generators,N_states)) dressed_H_matrix = 0.d0 do k = 1, N_det_generators do l = 1, N_int @@ -94,9 +95,20 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) enddo enddo !call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input) - call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) - if(do_it_perturbative)then - if(is_ok)then + call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop,is_ok_perturbative) +!do m = 1, N_states +! do k = 1, N_det_generators +! do l = 1, N_int +! psi_selectors(l,1,k) = psi_det_generators_input(l,1,k) +! psi_selectors(l,2,k) = psi_det_generators_input(l,2,k) +! enddo +! psi_selectors_coef(k,m) = psi_coef_diagonalized_tmp(k,m) +! enddo +!enddo +!soft_touch psi_selectors psi_selectors_coef +!if(do_it_perturbative)then + print*, 'is_ok_perturbative',is_ok_perturbative + if(is_ok.or.is_ok_perturbative)then N_det = N_det_generators do m = 1, N_states do k = 1, N_det_generators @@ -105,11 +117,19 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) psi_det(l,2,k) = psi_det_generators_input(l,2,k) enddo psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m) + print*, 'psi_coef(k,m)',psi_coef(k,m) + enddo + enddo + soft_touch psi_det psi_coef N_det + e_pt2 = 0.d0 + do m =1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators_input(1,1,m),psi_det_generators_input(1,1,l),N_int,hij) ! Fill the zeroth order H matrix + e_pt2 += (dressed_H_matrix(m,l) - hij)* psi_coef_diagonalized_tmp(m,1)* psi_coef_diagonalized_tmp(l,1) enddo enddo - touch psi_coef psi_det N_det endif - endif +!endif deallocate(psi_det_generators_input,dressed_H_matrix,psi_coef_diagonalized_tmp) @@ -118,14 +138,14 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) end -subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) +subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop,is_ok_perturbative) use bitmasks implicit none integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) integer, intent(in) :: Ndet_generators double precision, intent(in) :: threshold logical, intent(in) :: verbose - logical, intent(out) :: is_ok,exit_loop + logical, intent(out) :: is_ok,exit_loop,is_ok_perturbative double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states) double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators) @@ -309,10 +329,124 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener exit endif enddo + if(.not.is_ok)then + is_ok_perturbative = .True. + do i = 1, Ndet_generators + if(is_a_ref_det(i))cycle + do k = 1, N_states + print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative + if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold_perturbative)then + is_ok_perturbative = .False. + exit + endif + enddo + if(.not.is_ok_perturbative)then + exit + endif + enddo + endif if(verbose)then - print*,'is_ok = ',is_ok + print*,'is_ok = ',is_ok + print*,'is_ok_perturbative = ',is_ok_perturbative endif end +subroutine fill_H_apply_buffer_no_selection_first_order_coef(n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + BEGIN_DOC + ! Fill the H_apply buffer with determiants for CISD + END_DOC + + integer, intent(in) :: n_selected, Nint, iproc + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k + integer :: new_size + PROVIDE H_apply_buffer_allocated + call omp_set_lock(H_apply_buffer_lock(1,iproc)) + new_size = H_apply_buffer(iproc)%N_det + n_selected + if (new_size > H_apply_buffer(iproc)%sze) then + call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc) + endif + do i=1,H_apply_buffer(iproc)%N_det + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) + enddo + do i=1,n_selected + do j=1,N_int + H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i) + H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i) + enddo + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num) + enddo + double precision :: i_H_psi_array(N_states),h,diag_H_mat_elem_fock,delta_e + do i=1,N_selected + call i_H_psi(det_buffer(1,1,i),psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_psi_array) + call i_H_j(det_buffer(1,1,i),det_buffer(1,1,i),N_int,h) + do j=1,N_states + delta_e = -1.d0 /(h - CI_expectation_value(j)) + H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = i_H_psi_array(j) * delta_e + enddo + enddo + H_apply_buffer(iproc)%N_det = new_size + do i=1,H_apply_buffer(iproc)%N_det + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) + enddo + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) +end + + + +subroutine make_s2_eigenfunction_first_order + implicit none + integer :: i,j,k + integer :: smax, s + integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:) + integer :: N_det_new + integer, parameter :: bufsze = 1000 + logical, external :: is_in_wavefunction + + allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) ) + smax = 1 + N_det_new = 0 + + do i=1,N_occ_pattern + call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int) + s += 1 + if (s > smax) then + deallocate(d) + allocate ( d(N_int,2,s) ) + smax = s + endif + call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) + do j=1,s + if (.not. is_in_wavefunction(d(1,1,j), N_int) ) then + N_det_new += 1 + do k=1,N_int + det_buffer(k,1,N_det_new) = d(k,1,j) + det_buffer(k,2,N_det_new) = d(k,2,j) + enddo + if (N_det_new == bufsze) then + call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0) + N_det_new = 0 + endif + endif + enddo + enddo + + if (N_det_new > 0) then + call fill_H_apply_buffer_no_selection_first_order_coef(N_det_new,det_buffer,N_int,0) + call copy_H_apply_buffer_to_wf + SOFT_TOUCH N_det psi_coef psi_det + endif + + deallocate(d,det_buffer) + + call write_int(output_determinants,N_det_new, 'Added deteminants for S^2') + +end + diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f index 8656b633..8a709154 100644 --- a/plugins/FOBOCI/fobo_scf.irp.f +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -1,8 +1,13 @@ program foboscf implicit none - call run_prepare +!if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Read" )then +! disk_access_ao_integrals = "Write" +! touch disk_access_ao_integrals +!endif +!print*, 'disk_access_ao_integrals',disk_access_ao_integrals no_oa_or_av_opt = .True. touch no_oa_or_av_opt + call run_prepare call routine_fobo_scf call save_mos @@ -10,8 +15,8 @@ end subroutine run_prepare implicit none - no_oa_or_av_opt = .False. - touch no_oa_or_av_opt +! no_oa_or_av_opt = .False. +! touch no_oa_or_av_opt call damping_SCF call diag_inactive_virt_and_update_mos end @@ -27,6 +32,7 @@ subroutine routine_fobo_scf print*,'*******************************************************************************' print*,'*******************************************************************************' print*,'FOBO-SCF Iteration ',i + print*, 'ao_bielec_integrals_in_map = ',ao_bielec_integrals_in_map print*,'*******************************************************************************' print*,'*******************************************************************************' if(speed_up_convergence_foboscf)then @@ -46,7 +52,7 @@ subroutine routine_fobo_scf soft_touch threshold_lmct threshold_mlct endif endif - call FOBOCI_lmct_mlct_old_thr + call FOBOCI_lmct_mlct_old_thr(i) call save_osoci_natural_mos call damping_SCF call diag_inactive_virt_and_update_mos diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index a072a918..46ca9662 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -1,7 +1,8 @@ -subroutine FOBOCI_lmct_mlct_old_thr +subroutine FOBOCI_lmct_mlct_old_thr(iter) use bitmasks implicit none + integer, intent(in) :: iter integer :: i,j,k,l integer(bit_kind),allocatable :: unpaired_bitmask(:,:) integer, allocatable :: occ(:,:) @@ -10,7 +11,7 @@ subroutine FOBOCI_lmct_mlct_old_thr logical :: test_sym double precision :: thr,hij double precision, allocatable :: dressing_matrix(:,:) - logical :: verbose,is_ok + logical :: verbose,is_ok,is_ok_perturbative verbose = .True. thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) @@ -46,89 +47,45 @@ subroutine FOBOCI_lmct_mlct_old_thr i_hole_osoci = list_inact(i) print*,'--------------------------' ! First set the current generators to the one of restart - call set_generators_to_generators_restart - call set_psi_det_to_generators call check_symetry(i_hole_osoci,thr,test_sym) if(.not.test_sym)cycle + call set_generators_to_generators_restart + call set_psi_det_to_generators print*,'i_hole_osoci = ',i_hole_osoci call create_restart_and_1h(i_hole_osoci) call set_generators_to_psi_det print*,'Passed set generators' call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) - call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop) + double precision :: e_pt2 + call is_a_good_candidate(threshold_lmct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) print*,'is_ok = ',is_ok - if(.not.is_ok)cycle - allocate(dressing_matrix(N_det_generators,N_det_generators)) - dressing_matrix = 0.d0 - if(.not.do_it_perturbative)then - - 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 + if(is_ok)then + allocate(dressing_matrix(N_det_generators,N_det_generators)) + 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 - enddo - 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 - call make_s2_eigenfunction - call diagonalize_ci -! if(dressing_2h2p)then -! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct) -! endif - -! ! 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 + 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(e_pt2) + call make_s2_eigenfunction_first_order + threshold_davidson = 1.d-6 + soft_touch threshold_davidson davidson_criterion + call diagonalize_ci double precision :: hkl call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators) hkl = dressing_matrix(1,1) @@ -139,7 +96,10 @@ subroutine FOBOCI_lmct_mlct_old_thr 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) + deallocate(dressing_matrix) + else + if(.not.do_it_perturbative)cycle + if(.not. is_ok_perturbative)cycle endif call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) @@ -148,7 +108,6 @@ subroutine FOBOCI_lmct_mlct_old_thr norm_total(k) += norm_tmp(k) enddo call update_density_matrix_osoci - deallocate(dressing_matrix) enddo if(.True.)then @@ -163,9 +122,9 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'--------------------------' ! First set the current generators to the one of restart call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle call set_generators_to_generators_restart call set_psi_det_to_generators - if(.not.test_sym)cycle print*,'i_particl_osoci= ',i_particl_osoci ! Initialize the bitmask to the restart ones call initialize_bitmask_to_restart_ones @@ -181,32 +140,33 @@ subroutine FOBOCI_lmct_mlct_old_thr call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) !! ! so all the mono excitation on the new generators - call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop) + call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) print*,'is_ok = ',is_ok - if(.not. is_ok)then + if(is_ok)then + allocate(dressing_matrix(N_det_generators,N_det_generators)) + dressing_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) + dressing_matrix(k,l) = hkl + enddo + enddo + call all_single(e_pt2) + call make_s2_eigenfunction_first_order + threshold_davidson = 1.d-6 + soft_touch threshold_davidson davidson_criterion + + call diagonalize_ci + deallocate(dressing_matrix) + else if(exit_loop)then + call set_generators_to_generators_restart + call set_psi_det_to_generators exit else - cycle + if(.not.do_it_perturbative)cycle + if(.not. is_ok_perturbative)cycle endif - endif - allocate(dressing_matrix(N_det_generators,N_det_generators)) - if(.not.do_it_perturbative)then - dressing_matrix = 0.d0 - do k = 1, N_det_generators - do l = 1, N_det_generators - call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) - dressing_matrix(k,l) = hkl - enddo - enddo - ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) - ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) - call all_single - call make_s2_eigenfunction - call diagonalize_ci -! if(dressing_2h2p)then -! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct) -! endif endif call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) do k = 1, N_states @@ -214,7 +174,6 @@ subroutine FOBOCI_lmct_mlct_old_thr norm_total(k) += norm_tmp(k) enddo call update_density_matrix_osoci - deallocate(dressing_matrix) enddo endif @@ -376,3 +335,303 @@ subroutine FOBOCI_lmct_old enddo print*,'accu = ',accu end + +subroutine FOBOCI_lmct_mlct_old_thr_restart(iter) + use bitmasks + implicit none + integer, intent(in) :: iter + integer :: i,j,k,l + integer(bit_kind),allocatable :: unpaired_bitmask(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + double precision :: norm_tmp(N_states),norm_total(N_states) + logical :: test_sym + double precision :: thr,hij + double precision, allocatable :: dressing_matrix(:,:) + logical :: verbose,is_ok,is_ok_perturbative + verbose = .True. + thr = 1.d-12 + allocate(unpaired_bitmask(N_int,2)) + allocate (occ(N_int*bit_kind_size,2)) + do i = 1, N_int + unpaired_bitmask(i,1) = unpaired_alpha_electrons(i) + unpaired_bitmask(i,2) = unpaired_alpha_electrons(i) + enddo + norm_total = 0.d0 + call initialize_density_matrix_osoci + call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int) + print*,'' + print*,'' + print*,'mulliken spin population analysis' + accu =0.d0 + do i = 1, nucl_num + accu += mulliken_spin_densities(i) + print*,i,nucl_charge(i),mulliken_spin_densities(i) + enddo + print*,'' + print*,'' + print*,'DOING FIRST LMCT !!' + print*,'Threshold_lmct = ',threshold_lmct + integer(bit_kind) , allocatable :: zero_bitmask(:,:) + integer(bit_kind) , allocatable :: psi_singles(:,:,:) + logical :: lmct + double precision, allocatable :: psi_singles_coef(:,:) + logical :: exit_loop + allocate( zero_bitmask(N_int,2) ) + if(iter.ne.1)then + do i = 1, n_inact_orb + lmct = .True. + integer :: i_hole_osoci + i_hole_osoci = list_inact(i) + print*,'--------------------------' + ! First set the current generators to the one of restart + call check_symetry(i_hole_osoci,thr,test_sym) + if(.not.test_sym)cycle + call set_generators_to_generators_restart + call set_psi_det_to_generators + print*,'i_hole_osoci = ',i_hole_osoci + call create_restart_and_1h(i_hole_osoci) + call set_generators_to_psi_det + print*,'Passed set generators' + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) + double precision :: e_pt2 + call is_a_good_candidate(threshold_lmct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) + print*,'is_ok = ',is_ok + if(is_ok)then + allocate(dressing_matrix(N_det_generators,N_det_generators)) + 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 + 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(e_pt2) + call make_s2_eigenfunction_first_order + threshold_davidson = 1.d-6 + soft_touch threshold_davidson davidson_criterion + call diagonalize_ci + 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 + deallocate(dressing_matrix) + else + if(.not.do_it_perturbative)cycle + if(.not. is_ok_perturbative)cycle + endif + call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci) + + do k = 1, N_states + print*,'norm_tmp = ',norm_tmp(k) + norm_total(k) += norm_tmp(k) + enddo + call update_density_matrix_osoci + enddo + else + double precision :: array_dm(mo_tot_num) + call read_dm_from_lmct(array_dm) + call update_density_matrix_beta_osoci_read(array_dm) + endif + + if(iter.ne.1)then + if(.True.)then + print*,'' + print*,'DOING THEN THE MLCT !!' + print*,'Threshold_mlct = ',threshold_mlct + lmct = .False. + do i = 1, n_virt_orb + integer :: i_particl_osoci + i_particl_osoci = list_virt(i) + + print*,'--------------------------' + ! First set the current generators to the one of restart + call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle + call set_generators_to_generators_restart + call set_psi_det_to_generators + print*,'i_particl_osoci= ',i_particl_osoci + ! Initialize the bitmask to the restart ones + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call modify_bitmasks_for_particl(i_particl_osoci) + call print_generators_bitmasks_holes + ! Impose that only the active part can be reached + call set_bitmask_hole_as_input(unpaired_bitmask) +!!! call all_single_h_core + call create_restart_and_1p(i_particl_osoci) +!!! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) +!!! ! so all the mono excitation on the new generators + call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) + print*,'is_ok = ',is_ok + if(is_ok)then + allocate(dressing_matrix(N_det_generators,N_det_generators)) + dressing_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) + dressing_matrix(k,l) = hkl + enddo + enddo + call all_single(e_pt2) + call make_s2_eigenfunction_first_order + threshold_davidson = 1.d-6 + soft_touch threshold_davidson davidson_criterion + + call diagonalize_ci + deallocate(dressing_matrix) + else + if(exit_loop)then + call set_generators_to_generators_restart + call set_psi_det_to_generators + exit + else + if(.not.do_it_perturbative)cycle + if(.not. is_ok_perturbative)cycle + endif + endif + call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) + do k = 1, N_states + print*,'norm_tmp = ',norm_tmp(k) + norm_total(k) += norm_tmp(k) + enddo + call update_density_matrix_osoci + enddo + endif + else + integer :: norb + call read_dm_from_mlct(array_dm,norb) + call update_density_matrix_alpha_osoci_read(array_dm) + do i = norb+1, n_virt_orb + i_particl_osoci = list_virt(i) + + print*,'--------------------------' + ! First set the current generators to the one of restart + call check_symetry(i_particl_osoci,thr,test_sym) + if(.not.test_sym)cycle + call set_generators_to_generators_restart + call set_psi_det_to_generators + print*,'i_particl_osoci= ',i_particl_osoci + ! Initialize the bitmask to the restart ones + call initialize_bitmask_to_restart_ones + ! Impose that only the hole i_hole_osoci can be done + call modify_bitmasks_for_particl(i_particl_osoci) + call print_generators_bitmasks_holes + ! Impose that only the active part can be reached + call set_bitmask_hole_as_input(unpaired_bitmask) +!!! call all_single_h_core + call create_restart_and_1p(i_particl_osoci) +!!! ! Update the generators + call set_generators_to_psi_det + call set_bitmask_particl_as_input(reunion_of_bitmask) + call set_bitmask_hole_as_input(reunion_of_bitmask) +!!! ! so all the mono excitation on the new generators + call is_a_good_candidate(threshold_mlct,is_ok,e_pt2,verbose,exit_loop,is_ok_perturbative) + print*,'is_ok = ',is_ok + if(is_ok)then + allocate(dressing_matrix(N_det_generators,N_det_generators)) + dressing_matrix = 0.d0 + do k = 1, N_det_generators + do l = 1, N_det_generators + call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl) + dressing_matrix(k,l) = hkl + enddo + enddo + call all_single(e_pt2) + call make_s2_eigenfunction_first_order + threshold_davidson = 1.d-6 + soft_touch threshold_davidson davidson_criterion + + call diagonalize_ci + deallocate(dressing_matrix) + else + if(exit_loop)then + call set_generators_to_generators_restart + call set_psi_det_to_generators + exit + else + if(.not.do_it_perturbative)cycle + if(.not. is_ok_perturbative)cycle + endif + endif + call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) + do k = 1, N_states + print*,'norm_tmp = ',norm_tmp(k) + norm_total(k) += norm_tmp(k) + enddo + call update_density_matrix_osoci + enddo + endif + + print*,'norm_total = ',norm_total + norm_total = norm_generators_restart + norm_total = 1.d0/norm_total +! call rescale_density_matrix_osoci(norm_total) + double precision :: accu + accu = 0.d0 + do i = 1, mo_tot_num + accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i) + enddo + print*,'accu = ',accu +end + +subroutine read_dm_from_lmct(array) + implicit none + integer :: i,iunit ,getUnitAndOpen + double precision :: stuff + double precision, intent(out) :: array(mo_tot_num) + character*(128) :: input + input=trim("fort.33") + iunit= getUnitAndOpen(input,'r') + print*, iunit + array = 0.d0 + do i = 1, n_inact_orb + read(iunit,*) stuff + print*, list_inact(i),stuff + array(list_inact(i)) = stuff + enddo +end + +subroutine read_dm_from_mlct(array,norb) + implicit none + integer :: i,iunit ,getUnitAndOpen + double precision :: stuff + double precision, intent(out) :: array(mo_tot_num) + character*(128) :: input + input=trim("fort.35") + iunit= getUnitAndOpen(input,'r') + integer,intent(out) :: norb + read(iunit,*)norb + print*, iunit + input=trim("fort.34") + iunit= getUnitAndOpen(input,'r') + array = 0.d0 + print*, 'norb = ',norb + do i = 1, norb + read(iunit,*) stuff + print*, list_virt(i),stuff + array(list_virt(i)) = stuff + enddo +end diff --git a/plugins/FOBOCI/generators_restart_save.irp.f b/plugins/FOBOCI/generators_restart_save.irp.f index 09d4aa2b..eba9f0ad 100644 --- a/plugins/FOBOCI/generators_restart_save.irp.f +++ b/plugins/FOBOCI/generators_restart_save.irp.f @@ -9,6 +9,7 @@ BEGIN_PROVIDER [ integer, N_det_generators_restart ] integer :: i integer, save :: ifirst = 0 double precision :: norm + print*, ' Providing N_det_generators_restart' if(ifirst == 0)then call ezfio_get_determinants_n_det(N_det_generators_restart) ifirst = 1 @@ -30,6 +31,7 @@ END_PROVIDER integer :: i, k integer, save :: ifirst = 0 double precision, allocatable :: psi_coef_read(:,:) + print*, ' Providing psi_det_generators_restart' if(ifirst == 0)then call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart) do k = 1, N_int diff --git a/plugins/FOBOCI/hcc_1h1p.irp.f b/plugins/FOBOCI/hcc_1h1p.irp.f index 66cf2fd4..ffad686f 100644 --- a/plugins/FOBOCI/hcc_1h1p.irp.f +++ b/plugins/FOBOCI/hcc_1h1p.irp.f @@ -15,12 +15,12 @@ subroutine routine 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) + 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) @@ -34,12 +34,12 @@ end subroutine pouet implicit none double precision :: accu,coef_hf -! provide one_body_dm_mo_alpha one_body_dm_mo_beta + 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 + call save_wavefunction end diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 3ecd7977..b3dfca52 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -218,6 +218,44 @@ subroutine update_density_matrix_osoci enddo +end + +subroutine update_density_matrix_beta_osoci_read(array) + implicit none + BEGIN_DOC + ! one_body_dm_mo_alpha_osoci += Delta rho alpha + ! one_body_dm_mo_beta_osoci += Delta rho beta + END_DOC + integer :: i,j + integer :: iorb,jorb + double precision :: array(mo_tot_num) + do i = 1, mo_tot_num + j = list_act(1) + one_body_dm_mo_beta_osoci(i,j) += array(i) + one_body_dm_mo_beta_osoci(j,i) += array(i) + one_body_dm_mo_beta_osoci(i,i) += array(i) * array(i) + enddo + + +end + +subroutine update_density_matrix_alpha_osoci_read(array) + implicit none + BEGIN_DOC + ! one_body_dm_mo_alpha_osoci += Delta rho alpha + ! one_body_dm_mo_beta_osoci += Delta rho beta + END_DOC + integer :: i,j + integer :: iorb,jorb + double precision :: array(mo_tot_num) + do i = 1, mo_tot_num + j = list_act(1) + one_body_dm_mo_alpha_osoci(i,j) += array(i) + one_body_dm_mo_alpha_osoci(j,i) += array(i) + one_body_dm_mo_alpha_osoci(i,i) += array(i) * array(i) + enddo + + end @@ -387,14 +425,14 @@ subroutine save_osoci_natural_mos print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then + if(dabs(tmp(iorb,jorb)).gt.0.0001d0)then print*,'INACTIVE ' 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_mlct)then + if(dabs(tmp(iorb,jorb)).gt.0.0001d0)then print*,'VIRT ' print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif @@ -412,6 +450,10 @@ subroutine save_osoci_natural_mos label = "Natural" call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) +!if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then +! disk_access_ao_integrals = "Read" +! touch disk_access_ao_integrals +!endif !soft_touch mo_coef deallocate(tmp,occ) diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index c1da3670..80739aa2 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -73,21 +73,21 @@ print*, '1h1p = ',accu ! 1h1p third order - delta_ij_tmp = 0.d0 - call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) -!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) -!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + if(do_third_order_1h1p)then + delta_ij_tmp = 0.d0 + call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det + do j = 1, N_det + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + enddo enddo - enddo - second_order_pt_new_1h1p(i_state) = accu(i_state) - enddo - print*, '1h1p(3)',accu + second_order_pt_new_1h1p(i_state) = accu(i_state) + enddo + print*, '1h1p(3)',accu + endif ! 2h delta_ij_tmp = 0.d0 diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index b4c7e6f4..f08af1d5 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -287,50 +287,46 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) integer :: ispin,jspin,kspin if (n_holes_act == 0 .and. n_particles_act == 1) then -! i_particle_act = particles_active_list_spin_traced(1) -! delta_e_act += one_creat_spin_trace(i_particle_act ) ispin = particle_list_practical(1,1) i_particle_act = particle_list_practical(2,1) - call get_excitation_degree(det_1,det_2,degree,N_int) - if(degree == 1)then - call get_excitation(det_1,det_2,exc,degree,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - i_hole = list_inact_reverse(h1) - i_part = list_act_reverse(p1) - do i_state = 1, N_states - delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) - enddo - else if (degree == 2)then +! call get_excitation_degree(det_1,det_2,degree,N_int) +! if(degree == 1)then +! call get_excitation(det_1,det_2,exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! i_hole = list_inact_reverse(h1) +! i_part = list_act_reverse(p1) +! do i_state = 1, N_states +! delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) +! enddo +! else if (degree == 2)then do i_state = 1, N_states delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) enddo - endif +! endif else if (n_holes_act == 1 .and. n_particles_act == 0) then -! i_hole_act = holes_active_list_spin_traced(1) -! delta_e_act += one_anhil_spin_trace(i_hole_act ) ispin = hole_list_practical(1,1) i_hole_act = hole_list_practical(2,1) - call get_excitation_degree(det_1,det_2,degree,N_int) - if(degree == 1)then - call get_excitation(det_1,det_2,exc,degree,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - i_hole = list_act_reverse(h1) - i_part = list_virt_reverse(p1) - do i_state = 1, N_states - if(isnan(one_creat_virt(i_hole,i_part,i_state)))then - print*, i_hole,i_part,i_state - call debug_det(det_1,N_int) - call debug_det(det_2,N_int) - stop - endif - delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state) - enddo - else if (degree == 2)then +! call get_excitation_degree(det_1,det_2,degree,N_int) +! if(degree == 1)then +! call get_excitation(det_1,det_2,exc,degree,phase,N_int) +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! i_hole = list_act_reverse(h1) +! i_part = list_virt_reverse(p1) +! do i_state = 1, N_states +! if(isnan(one_creat_virt(i_hole,i_part,i_state)))then +! print*, i_hole,i_part,i_state +! call debug_det(det_1,N_int) +! call debug_det(det_2,N_int) +! stop +! endif +! delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state) +! enddo +! else if (degree == 2)then do i_state = 1, N_states delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) enddo - endif +! endif else if (n_holes_act == 1 .and. n_particles_act == 1) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -460,21 +456,12 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) double precision :: phase call get_excitation_degree(det_1,det_2,degree,N_int) if(degree == 1)then -! call debug_det(det_1,N_int) call get_excitation(det_1,det_2,exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) i_hole = list_inact_reverse(h1) i_part = list_virt_reverse(p1) do i_state = 1, N_states -! if(one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1).gt.1.d-10)then -! print*, hij, one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) -! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) & -! * coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) -! print*, coef_array(i_state)* hij*coef_array(i_state)* hij,one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1), & -! coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) -! else - delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) -! endif +! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) enddo endif diff --git a/plugins/Perturbation/EZFIO.cfg b/plugins/Perturbation/EZFIO.cfg index ad26cfe5..c5d6379f 100644 --- a/plugins/Perturbation/EZFIO.cfg +++ b/plugins/Perturbation/EZFIO.cfg @@ -16,4 +16,4 @@ type: Normalized_float doc: The selection process stops when the energy ratio variational/(variational+PT2) is equal to var_pt2_ratio interface: ezfio,provider,ocaml -default: 0.75 \ No newline at end of file +default: 0.75 diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index c284e01e..29aed8d3 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -142,6 +142,60 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) end + + +subroutine pt2_epstein_nesbet_2x2_no_ci_diag($arguments) + use bitmasks + implicit none + $declarations + + BEGIN_DOC + ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + ! + ! for the various N_st states. + ! + ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + ! + ! c_pert(i) = e_2_pert(i)/ + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem_fock,delta_e, h + double precision :: i_H_psi_array(N_st) + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + PROVIDE CI_electronic_energy + + call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + !call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + do i =1,N_st + if (i_H_psi_array(i) /= 0.d0) then + delta_e = h - CI_expectation_value(i) + if (delta_e > 0.d0) then + e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + else + e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + endif + if (dabs(i_H_psi_array(i)) > 1.d-6) then + c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + else + c_pert(i) = 0.d0 + endif + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + else + e_2_pert(i) = 0.d0 + c_pert(i) = 0.d0 + H_pert_diag(i) = 0.d0 + endif + enddo + +end + + + subroutine pt2_moller_plesset ($arguments) use bitmasks implicit none diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f index de14da9f..63ad545b 100644 --- a/plugins/Properties/hyperfine_constants.irp.f +++ b/plugins/Properties/hyperfine_constants.irp.f @@ -121,6 +121,11 @@ END_PROVIDER conversion_factor_mhz_hcc(8) = -606.1958551736545d0 conversion_factor_gauss_hcc(8) = -216.30574771560407d0 conversion_factor_cm_1_hcc(8) = -202.20517197179822d0 + + ! Phosphore + conversion_factor_mhz_hcc(15) = 1811.0967763744873d0 + conversion_factor_gauss_hcc(15) = 646.2445276897648d0 + conversion_factor_cm_1_hcc(15) = 604.1170297381395d0 END_PROVIDER diff --git a/plugins/Selectors_no_sorted/selectors.irp.f b/plugins/Selectors_no_sorted/selectors.irp.f index 9273c7bb..83a8d472 100644 --- a/plugins/Selectors_no_sorted/selectors.irp.f +++ b/plugins/Selectors_no_sorted/selectors.irp.f @@ -40,6 +40,7 @@ END_PROVIDER do k=1,N_states do i=1,N_det_selectors psi_selectors_coef(i,k) = psi_coef(i,k) +! print*, 'psi_selectors_coef(i,k) == ',psi_selectors_coef(i,k) enddo enddo END_PROVIDER diff --git a/plugins/loc_cele/loc_exchange_int.irp.f b/plugins/loc_cele/loc_exchange_int.irp.f new file mode 100644 index 00000000..d7cc5c65 --- /dev/null +++ b/plugins/loc_cele/loc_exchange_int.irp.f @@ -0,0 +1,110 @@ +program loc_int + implicit none + integer :: i,j,k,l,iorb,jorb + double precision :: exchange_int(mo_tot_num) + integer :: iorder(mo_tot_num) + integer :: indices(mo_tot_num,2) + logical :: list_core_inact_check(mo_tot_num) + integer :: n_rot + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_core_inact_orb + jorb = list_core_inact(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_act_orb + iorb = list_act(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_act_orb + jorb = list_act(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_virt_orb + iorb = list_virt(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_virt_orb + jorb = list_virt(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + + + call save_mos + + +end diff --git a/plugins/loc_cele/loc_exchange_int_act.irp.f b/plugins/loc_cele/loc_exchange_int_act.irp.f new file mode 100644 index 00000000..b9bbeb82 --- /dev/null +++ b/plugins/loc_cele/loc_exchange_int_act.irp.f @@ -0,0 +1,45 @@ +program loc_int + implicit none + integer :: i,j,k,l,iorb,jorb + double precision :: exchange_int(mo_tot_num) + integer :: iorder(mo_tot_num) + integer :: indices(mo_tot_num,2) + logical :: list_core_inact_check(mo_tot_num) + integer :: n_rot + + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_act_orb + iorb = list_act(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_act_orb + jorb = list_act(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + call save_mos + + +end diff --git a/plugins/loc_cele/loc_exchange_int_inact.irp.f b/plugins/loc_cele/loc_exchange_int_inact.irp.f new file mode 100644 index 00000000..2ff3c85f --- /dev/null +++ b/plugins/loc_cele/loc_exchange_int_inact.irp.f @@ -0,0 +1,45 @@ +program loc_int + implicit none + integer :: i,j,k,l,iorb,jorb + double precision :: exchange_int(mo_tot_num) + integer :: iorder(mo_tot_num) + integer :: indices(mo_tot_num,2) + logical :: list_core_inact_check(mo_tot_num) + integer :: n_rot + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_core_inact_orb + jorb = list_core_inact(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + + call save_mos + + +end diff --git a/plugins/loc_cele/loc_exchange_int_virt.irp.f b/plugins/loc_cele/loc_exchange_int_virt.irp.f new file mode 100644 index 00000000..333f189b --- /dev/null +++ b/plugins/loc_cele/loc_exchange_int_virt.irp.f @@ -0,0 +1,47 @@ +program loc_int + implicit none + integer :: i,j,k,l,iorb,jorb + double precision :: exchange_int(mo_tot_num) + integer :: iorder(mo_tot_num) + integer :: indices(mo_tot_num,2) + logical :: list_core_inact_check(mo_tot_num) + integer :: n_rot + + indices = 0 + list_core_inact_check = .True. + n_rot = 0 + do i = 1, n_virt_orb + iorb = list_virt(i) + exchange_int = 0.d0 + iorder = 0 + print*,'' + if(list_core_inact_check(iorb) == .False.)cycle + do j = i+1, n_virt_orb + jorb = list_virt(j) + iorder(jorb) = jorb + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + enddo + n_rot += 1 + call dsort(exchange_int,iorder,mo_tot_num) + indices(n_rot,1) = iorb + indices(n_rot,2) = iorder(1) + list_core_inact_check(iorder(1)) = .False. + print*,indices(n_rot,1),indices(n_rot,2) + print*,'' + print*,'' + enddo + print*,'****************************' + print*,'-+++++++++++++++++++++++++' + do i = 1, n_rot + iorb = indices(i,1) + jorb = indices(i,2) + print*,iorb,jorb + call mix_mo_jk(iorb,jorb) + enddo + + + + call save_mos + + +end diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index b1c459ba..324005aa 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -124,3 +124,16 @@ interface: ezfio,provider,ocaml doc: Energy that should be obtained when truncating the wave function (optional) type: Energy default: 0. + +[store_full_H_mat] +type: logical +doc: If True, the Davidson diagonalization is performed by storring the full H matrix up to n_det_max_stored. Be carefull, it can cost a lot of memory but can also save a lot of CPU time +interface: ezfio,provider,ocaml +default: False + +[n_det_max_stored] +type: Det_number_max +doc: Maximum number of determinants for which the full H matrix is stored. Be carefull, the memory requested scales as 10*n_det_max_stored**2. For instance, 90000 determinants represent a matrix of size 60 Gb. +interface: ezfio,provider,ocaml +default: 90000 + diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index 28513597..d5b972e4 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -306,7 +306,6 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) call omp_unset_lock(H_apply_buffer_lock(1,iproc)) end - subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) use f77_zmq implicit none diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index deba43c5..130bd56d 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -334,6 +334,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun double precision :: to_print(2,N_st) double precision :: cpu, wall + if(store_full_H_mat.and.sze.le.n_det_max_stored)then + provide H_matrix_all_dets + endif call write_time(iunit) @@ -439,7 +442,11 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! ---------------------- do k=1,N_st - call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + if(store_full_H_mat.and.sze.le.n_det_max_stored)then + call H_u_0_stored(W(1,k,iter),U(1,k,iter),H_matrix_all_dets,sze) + else + call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + endif enddo diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f index 11ec6db5..49714082 100644 --- a/src/Determinants/diagonalize_CI.irp.f +++ b/src/Determinants/diagonalize_CI.irp.f @@ -33,6 +33,14 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_expectation_value, (N_states_diag) ] + implicit none + integer :: i + do i = 1, N_states + call u0_H_u_0(CI_expectation_value(i),psi_coef(1,i),n_det,psi_det,N_int) + enddo + END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ] @@ -69,10 +77,14 @@ END_PROVIDER if (diag_algorithm == "Davidson") then + print*, '------------- In Davidson ' call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, & size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_determinants) + print*, '------------- Out Davidson ' do j=1,N_states_diag + print*, '------------- In S^2' call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),CI_eigenvectors_s2(j)) + print*, '------------- Out S^2' enddo @@ -233,16 +245,20 @@ END_PROVIDER else - ! Sorting the N_states_diag by energy, whatever the S^2 value is + !! Sorting the N_states_diag by energy, whatever the S^2 value is allocate(e_array(n_states_diag),iorder(n_states_diag)) - do j = 1, N_states_diag - call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) + do j = 2, N_states_diag + if(store_full_H_mat.and.n_det.le.n_det_max_stored)then + call u_0_H_u_0_stored(e_0,CI_eigenvectors(1,j),H_matrix_all_dets,n_det) + else + call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) + endif e_array(j) = e_0 iorder(j) = j enddo call dsort(e_array,iorder,n_states_diag) - do j = 1, N_states_diag + do j = 2, N_states_diag CI_electronic_energy(j) = e_array(j) do i = 1, N_det CI_eigenvectors(i,j) = psi_coef(i,iorder(j)) @@ -253,6 +269,7 @@ END_PROVIDER endif deallocate(s2_eigvalues) endif + print*, 'out provider' END_PROVIDER diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index e2e12974..c58d1f82 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -256,27 +256,6 @@ subroutine make_s2_eigenfunction integer :: N_det_new integer, parameter :: bufsze = 1000 logical, external :: is_in_wavefunction -! return - -! !TODO DEBUG -! do i=1,N_det -! do j=i+1,N_det -! s = 0 -! do k=1,N_int -! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & -! (psi_det(k,2,j) /= psi_det(k,2,i))) then -! s=1 -! exit -! endif -! enddo -! if ( s == 0 ) then -! print *, 'Error0: det ', j, 'already in wf' -! call debug_det(psi_det(1,1,j),N_int) -! stop -! endif -! enddo -! enddo -! !TODO DEBUG allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) ) smax = 1 @@ -308,33 +287,15 @@ subroutine make_s2_eigenfunction if (N_det_new > 0) then call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0) +! call fill_H_apply_buffer_no_selection_first_order_coef(N_det_new,det_buffer,N_int,0) call copy_H_apply_buffer_to_wf SOFT_TOUCH N_det psi_coef psi_det endif deallocate(d,det_buffer) - -! !TODO DEBUG -! do i=1,N_det -! do j=i+1,N_det -! s = 0 -! do k=1,N_int -! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & -! (psi_det(k,2,j) /= psi_det(k,2,i))) then -! s=1 -! exit -! endif -! enddo -! if ( s == 0 ) then -! print *, 'Error : det ', j, 'already in wf at ', i -! call debug_det(psi_det(1,1,j),N_int) -! stop -! endif -! enddo -! enddo -! !TODO DEBUG call write_int(output_determinants,N_det_new, 'Added deteminants for S^2') end + diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 7f0e7e57..5cbed15e 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -431,7 +431,7 @@ end -subroutine i_H_j(key_i,key_j,Nint,hij) +subroutine i_H_j_new(key_i,key_j,Nint,hij) use bitmasks implicit none BEGIN_DOC @@ -463,6 +463,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) hij = 0.d0 !DIR$ FORCEINLINE call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin select case (degree) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) @@ -507,6 +508,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) ! Mono alpha m = exc(1,1,1) p = exc(1,2,1) + spin = 1 do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then @@ -534,6 +536,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij) ! Mono beta m = exc(1,1,2) p = exc(1,2,2) + spin = 2 + do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then @@ -559,6 +563,154 @@ subroutine i_H_j(key_i,key_j,Nint,hij) endif hij = phase*(hij + mo_mono_elec_integral(m,p)) + + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + +subroutine i_H_j(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral_schwartz + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 +! do k = 1, elec_alpha_num +! i = occ(k,1) +! if (.not.has_mipi(i)) then +! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) +! miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) +! has_mipi(i) = .True. +! endif +! enddo +! do k = 1, elec_beta_num +! i = occ(k,2) +! if (.not.has_mipi(i)) then +! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) +! has_mipi(i) = .True. +! endif +! enddo +! +! do k = 1, elec_alpha_num +! hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) +! enddo +! do k = 1, elec_beta_num +! hij = hij + mipi(occ(k,2)) +! enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + +! do k = 1, elec_beta_num +! i = occ(k,2) +! if (.not.has_mipi(i)) then +! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) +! miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) +! has_mipi(i) = .True. +! endif +! enddo +! do k = 1, elec_alpha_num +! i = occ(k,1) +! if (.not.has_mipi(i)) then +! mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) +! has_mipi(i) = .True. +! endif +! enddo +! +! do k = 1, elec_alpha_num +! hij = hij + mipi(occ(k,1)) +! enddo +! do k = 1, elec_beta_num +! hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) +! enddo + + endif +! hij = phase*(hij + mo_mono_elec_integral(m,p)) + + call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) case (0) hij = diag_H_mat_elem(key_i,Nint) @@ -2182,3 +2334,43 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) deallocate (shortcut, sort_idx, sorted, version) end +subroutine H_u_0_stored(v_0,u_0,hmatrix,sze) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! uses the big_matrix_stored array + END_DOC + integer, intent(in) :: sze + double precision, intent(in) :: hmatrix(sze,sze) + double precision, intent(out) :: v_0(sze) + double precision, intent(in) :: u_0(sze) + v_0 = 0.d0 + call matrix_vector_product(u_0,v_0,hmatrix,sze,sze) + +end + +subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze) + use bitmasks + implicit none + BEGIN_DOC + ! Computes e_0 = + ! + ! n : number of determinants + ! + ! uses the big_matrix_stored array + END_DOC + integer, intent(in) :: sze + double precision, intent(in) :: hmatrix(sze,sze) + double precision, intent(out) :: e_0 + double precision, intent(in) :: u_0(sze) + double precision :: v_0(sze) + double precision :: u_dot_v + e_0 = 0.d0 + v_0 = 0.d0 + call matrix_vector_product(u_0,v_0,hmatrix,sze,sze) + e_0 = u_dot_v(v_0,u_0,sze) +end diff --git a/src/Determinants/utils.irp.f b/src/Determinants/utils.irp.f index 22faee83..cc191970 100644 --- a/src/Determinants/utils.irp.f +++ b/src/Determinants/utils.irp.f @@ -1,15 +1,21 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] + use bitmasks implicit none BEGIN_DOC ! H matrix on the basis of the slater determinants defined by psi_det END_DOC - integer :: i,j + integer :: i,j,k double precision :: hij + integer :: degree(N_det),idx(0:N_det) call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) - !$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij) & + !$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij,degree,idx,k) & !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets) do i =1,N_det - do j =i,N_det +! call get_excitation_degree_vector(psi_det,psi_det(1,1,i),degree,N_int,N_det,idx) +! do k =1, idx(0) +! j = idx(k) +! if(j.lt.i)cycle + do j = i, N_det call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) H_matrix_all_dets(i,j) = hij H_matrix_all_dets(j,i) = hij @@ -18,3 +24,33 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] !$OMP END PARALLEL DO END_PROVIDER + +subroutine provide_big_matrix_stored_with_current_dets(sze,dets_in,big_matrix_stored) + use bitmasks + integer, intent(in) :: sze + integer(bit_kind), intent(in) :: dets_in(N_int,2,sze) + double precision, intent(out) :: big_matrix_stored(sze,sze) + integer :: i,j,k + double precision :: hij + integer :: degree(N_det),idx(0:N_det) + call i_H_j(dets_in(1,1,1),dets_in(1,1,1),N_int,hij) + print*, 'providing big_matrix_stored' + print*, n_det_max_stored + !$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij,degree,idx,k) & + !$OMP SHARED (sze, dets_in, N_int,big_matrix_stored) + do i =1,sze +! call get_excitation_degree_vector(dets_in,dets_in(1,1,i),degree,N_int,sze,idx) +! do k =1, idx(0) +! j = idx(k) + do j = i, sze + if(j.lt.i)cycle + call i_H_j(dets_in(1,1,i),dets_in(1,1,j),N_int,hij) + big_matrix_stored(i,j) = hij + big_matrix_stored(j,i) = hij + enddo + enddo + !$OMP END PARALLEL DO + print*, 'big_matrix_stored provided !!' + + +end diff --git a/src/Integrals_Bielec/EZFIO.cfg b/src/Integrals_Bielec/EZFIO.cfg index 01b87fc1..0d5c5832 100644 --- a/src/Integrals_Bielec/EZFIO.cfg +++ b/src/Integrals_Bielec/EZFIO.cfg @@ -12,6 +12,22 @@ interface: ezfio,provider,ocaml default: False ezfio_name: no_vvvv_integrals + +[write_ao_map_after_transfo] +type: logical +doc: If True, you dump all the ao integrals after having transformed the mo integrals +interface: ezfio,provider,ocaml +default: False +ezfio_name: write_ao_map_after_transfo + +[clear_ao_map_after_mo_transfo] +type: logical +doc: If True, you clear all the ao integrals after having done the transformation +interface: ezfio,provider,ocaml +default: False +ezfio_name: clear_ao_map_after_mo_transfo + + [no_ivvv_integrals] type: logical doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual index and 1 belonging to the core inactive active orbitals @@ -19,6 +35,13 @@ interface: ezfio,provider,ocaml default: False ezfio_name: no_ivvv_integrals +[no_vvv_integrals] +type: logical +doc: Can be switched on only if no_vvvv_integrals is True, then do not computes the integrals having 3 virtual orbitals +interface: ezfio,provider,ocaml +default: False +ezfio_name: no_vvv_integrals + [disk_access_mo_integrals] type: Disk_access doc: Read/Write MO integrals from/to disk [ Write | Read | None ] diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index eb443701..54bcc1c4 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -349,6 +349,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integral = ao_bielec_integral(1,1,1,1) real :: map_mb + print*, 'read_ao_integrals',read_ao_integrals + print*, 'disk_access_ao_integrals',disk_access_ao_integrals if (read_ao_integrals) then integer :: load_ao_integrals print*,'Reading the AO integrals' diff --git a/src/Integrals_Bielec/integrals_3_index.irp.f b/src/Integrals_Bielec/integrals_3_index.irp.f new file mode 100644 index 00000000..b9ee29b9 --- /dev/null +++ b/src/Integrals_Bielec/integrals_3_index.irp.f @@ -0,0 +1,22 @@ + BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num_align,mo_tot_num, mo_tot_num)] +&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num_align,mo_tot_num, mo_tot_num)] + implicit none + integer :: i,j,k,l + double precision :: get_mo_bielec_integral_schwartz + double precision :: integral + + do k = 1, mo_tot_num + do i = 1, mo_tot_num + do j = 1, mo_tot_num + l = j + integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map) + big_array_coulomb_integrals(j,i,k) = integral + l = j + integral = get_mo_bielec_integral_schwartz(i,j,l,k,mo_integrals_map) + big_array_exchange_integrals(j,i,k) = integral + enddo + enddo + enddo + + +END_PROVIDER diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index dc35f278..305abee3 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -414,6 +414,73 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map) deallocate(pairs,hash,iorder,tmp_val) end +subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,l,i,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +subroutine get_mo_bielec_integrals_exch_ii(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1) + ! for k,l fixed. + END_DOC + integer, intent(in) :: k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_bielec_integrals_in_map + + integer :: kk + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(k,i,i,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(map, hash, out_val, sze) + else + call map_get_many(map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + + integer*8 function get_mo_map_size() implicit none BEGIN_DOC diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 85fdde10..bf23ad1f 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -20,6 +20,7 @@ end BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] +use map_module implicit none integer(bit_kind) :: mask_ijkl(N_int,4) integer(bit_kind) :: mask_ijk(N_int,3) @@ -40,7 +41,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] if(no_vvvv_integrals)then integer :: i,j,k,l -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 4 ! print*, '' @@ -52,7 +53,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) enddo call add_integrals_to_map(mask_ijkl) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 2 (virt) ^2 ! = J_iv print*, '' @@ -76,17 +77,19 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] mask_ijkl(i,4) = virt_bitmask(i,1) enddo call add_integrals_to_map(mask_ijkl) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! + if(.not.no_vvv_integrals)then print*, '' - print*, '' + print*, ' and ' do i = 1,N_int - mask_ijk(i,1) = virt_bitmask(i,1) - mask_ijk(i,2) = virt_bitmask(i,1) - mask_ijk(i,3) = virt_bitmask(i,1) + mask_ijk(i,1) = virt_bitmask(i,1) + mask_ijk(i,2) = virt_bitmask(i,1) + mask_ijk(i,3) = virt_bitmask(i,1) enddo call add_integrals_to_map_three_indices(mask_ijk) + endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 3 (virt) ^1 ! print*, '' @@ -101,9 +104,9 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 1 (virt) ^3 ! - print*, '' - print*, '' if(.not.no_ivvv_integrals)then + print*, '' + print*, '' do i = 1,N_int mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) mask_ijkl(i,2) = virt_bitmask(i,1) @@ -116,6 +119,21 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] else call add_integrals_to_map(full_ijkl_bitmask_4) endif + + if(write_ao_map_after_transfo)then + call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') + disk_access_ao_integrals = "Read" + touch disk_access_ao_integrals + call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read") + endif + if(clear_ao_map_after_mo_transfo)then + call clear_ao_map + integer (map_size_kind) :: get_ao_map_size + print*, '^^^^^^^^^^^^^^^^^^^^^' + print *, 'get_ao_map_size',get_ao_map_size + print*, '^^^^^^^^^^^^^^^^^^^^^' + FREE ao_bielec_integrals_in_map + endif END_PROVIDER subroutine set_integrals_jj_into_map @@ -391,64 +409,41 @@ subroutine add_integrals_to_map(mask_ijkl) endif j1 = ishft((l*l-l),-1) do j0 = 1, n_j -! print*, 'l :: j0',l j = list_ijkl(j0,2) -! print*, 'j :: 2',j if (j > l) then -! print*, 'j>l' -! print*, j,l exit endif j1 += 1 do k0 = 1, n_k k = list_ijkl(k0,3) -! print*, 'l :: k0',l -! print*, 'k :: 3',j i1 = ishft((k*k-k),-1) if (i1<=j1) then continue else -! print*, 'k>l' -! print*, k,l exit endif bielec_tmp_1 = 0.d0 do i0 = 1, n_i i = list_ijkl(i0,1) -! print*, 'l :: i0',l -! print*, 'i :: 1',i if (i>k) then -! print*, 'i>k' -! print*, i,k exit endif bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) ! i1+=1 enddo -! do i = 1, min(k,j1-i1+list_ijkl(1,1)) -! do i = 1, min(k,j1-i1+list_ijkl(1,1)-1) do i0 = 1, n_i i = list_ijkl(i0,1) if(i> min(k,j1-i1+list_ijkl(1,1)-1))then -! if (i>k) then !min(k,j1-i1) exit endif -! print*, i,j,k,l -! print*, k,j1,i1,j1-i1 if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then cycle endif -! print*, i,j,k,l n_integrals += 1 buffer_value(n_integrals) = bielec_tmp_1(i) !DEC$ FORCEINLINE call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) -! if(i==12.and.k==12 .and.j==12.and.l==12)then -! print*, i,j,k,l,buffer_i(n_integrals) -! accu_bis += buffer_value(n_integrals) -! print*, buffer_value(n_integrals),accu_bis -! endif if (n_integrals == size_buffer) then call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& real(mo_integrals_threshold,integral_kind)) @@ -631,7 +626,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk) bielec_tmp_2 = 0.d0 do j1 = 1,ao_num call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) - ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) enddo do j1 = 1,ao_num kmax = 0 @@ -732,9 +726,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk) j = list_ijkl(j0,2) do i0 = 1, n_i i = list_ijkl(i0,1) -! if(m==2)then -! if(i==j .and. j == k)cycle -! endif if (i>k) then exit endif diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 13138499..f9645aa4 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -561,3 +561,18 @@ end +subroutine matrix_vector_product(u0,u1,matrix,sze,lda) + implicit none + BEGIN_DOC +! performs u1 += u0 * matrix + END_DOC + integer, intent(in) :: sze,lda + double precision, intent(in) :: u0(sze) + double precision, intent(inout) :: u1(sze) + double precision, intent(in) :: matrix(lda,sze) + integer :: i,j + integer :: incx,incy + incx = 1 + incy = 1 + call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy) +end