From 98059e87f5a630dc26ef4cbe5102ce049fa6d76f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 20 Sep 2016 17:29:02 +0200 Subject: [PATCH] Extracting davidson from slater rules --- config/ifort.cfg | 2 +- plugins/All_singles/NEEDED_CHILDREN_MODULES | 2 +- plugins/CAS_SD/NEEDED_CHILDREN_MODULES | 2 +- plugins/CID/NEEDED_CHILDREN_MODULES | 2 +- plugins/CIS/NEEDED_CHILDREN_MODULES | 2 +- plugins/CISD/NEEDED_CHILDREN_MODULES | 2 +- plugins/Casino/NEEDED_CHILDREN_MODULES | 2 +- plugins/DDCI_selected/NEEDED_CHILDREN_MODULES | 2 +- plugins/FCIdump/NEEDED_CHILDREN_MODULES | 2 +- plugins/FOBOCI/NEEDED_CHILDREN_MODULES | 2 +- plugins/Full_CI/H_apply.irp.f | 22 -- plugins/Full_CI/NEEDED_CHILDREN_MODULES | 2 +- plugins/Perturbation/NEEDED_CHILDREN_MODULES | 2 +- .../Perturbation/delta_rho_perturbation.irp.f | 77 ------- plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES | 2 +- plugins/QmcChem/NEEDED_CHILDREN_MODULES | 2 +- src/Davidson/davidson.irp.f | 3 + .../diagonalization.irp.f} | 65 +----- .../diagonalize_CI.irp.f | 16 -- .../diagonalize_CI_SC2.irp.f | 0 .../diagonalize_CI_mono.irp.f | 0 src/Davidson/parameters.irp.f | 62 ++++++ src/Davidson/u0Hu0.irp.f | 190 ++++++++++++++++++ src/Determinants/determinants.irp.f | 19 +- src/Determinants/s2.irp.f | 2 +- src/Determinants/slater_rules.irp.f | 189 ----------------- 26 files changed, 291 insertions(+), 382 deletions(-) delete mode 100644 plugins/Perturbation/delta_rho_perturbation.irp.f create mode 100644 src/Davidson/davidson.irp.f rename src/{Determinants/davidson.irp.f => Davidson/diagonalization.irp.f} (89%) rename src/{Determinants => Davidson}/diagonalize_CI.irp.f (96%) rename src/{Determinants => Davidson}/diagonalize_CI_SC2.irp.f (100%) rename src/{Determinants => Davidson}/diagonalize_CI_mono.irp.f (100%) create mode 100644 src/Davidson/parameters.irp.f create mode 100644 src/Davidson/u0Hu0.irp.f diff --git a/config/ifort.cfg b/config/ifort.cfg index b04506d4..255438cc 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags # [OPT] FC : -traceback -FCFLAGS : -xHost -O2 -ip -ftz -g +FCFLAGS : -xHost -O0 -ip -ftz -g # Profiling flags ################# diff --git a/plugins/All_singles/NEEDED_CHILDREN_MODULES b/plugins/All_singles/NEEDED_CHILDREN_MODULES index bb97ddb9..ee0ff040 100644 --- a/plugins/All_singles/NEEDED_CHILDREN_MODULES +++ b/plugins/All_singles/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Generators_restart Perturbation Properties Selectors_no_sorted Utils +Generators_restart Perturbation Properties Selectors_no_sorted Utils Davidson diff --git a/plugins/CAS_SD/NEEDED_CHILDREN_MODULES b/plugins/CAS_SD/NEEDED_CHILDREN_MODULES index f7264a0f..0b7ce8a9 100644 --- a/plugins/CAS_SD/NEEDED_CHILDREN_MODULES +++ b/plugins/CAS_SD/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_CAS +Perturbation Selectors_full Generators_CAS Davidson diff --git a/plugins/CID/NEEDED_CHILDREN_MODULES b/plugins/CID/NEEDED_CHILDREN_MODULES index afc8cfd4..1632a44d 100644 --- a/plugins/CID/NEEDED_CHILDREN_MODULES +++ b/plugins/CID/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod +Selectors_full SingleRefMethod Davidson diff --git a/plugins/CIS/NEEDED_CHILDREN_MODULES b/plugins/CIS/NEEDED_CHILDREN_MODULES index afc8cfd4..1632a44d 100644 --- a/plugins/CIS/NEEDED_CHILDREN_MODULES +++ b/plugins/CIS/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod +Selectors_full SingleRefMethod Davidson diff --git a/plugins/CISD/NEEDED_CHILDREN_MODULES b/plugins/CISD/NEEDED_CHILDREN_MODULES index afc8cfd4..1632a44d 100644 --- a/plugins/CISD/NEEDED_CHILDREN_MODULES +++ b/plugins/CISD/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod +Selectors_full SingleRefMethod Davidson diff --git a/plugins/Casino/NEEDED_CHILDREN_MODULES b/plugins/Casino/NEEDED_CHILDREN_MODULES index aae89501..34de8ddb 100644 --- a/plugins/Casino/NEEDED_CHILDREN_MODULES +++ b/plugins/Casino/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants +Determinants Davidson diff --git a/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES b/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES index f7264a0f..0b7ce8a9 100644 --- a/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES +++ b/plugins/DDCI_selected/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_CAS +Perturbation Selectors_full Generators_CAS Davidson diff --git a/plugins/FCIdump/NEEDED_CHILDREN_MODULES b/plugins/FCIdump/NEEDED_CHILDREN_MODULES index aae89501..34de8ddb 100644 --- a/plugins/FCIdump/NEEDED_CHILDREN_MODULES +++ b/plugins/FCIdump/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants +Determinants Davidson diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES index f6c0c1c4..e40934be 100644 --- a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES +++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_no_sorted Hartree_Fock +Perturbation Selectors_no_sorted Hartree_Fock Davidson diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 596c947a..9e22b965 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -30,28 +30,6 @@ s.unset_openmp() print s -s = H_apply("select_mono_delta_rho") -s.unset_double_excitations() -s.set_selection_pt2("delta_rho_one_point") -s.unset_openmp() -print s - -s = H_apply("pt2_mono_delta_rho") -s.unset_double_excitations() -s.set_perturbation("delta_rho_one_point") -s.unset_openmp() -print s - -s = H_apply("select_mono_di_delta_rho") -s.set_selection_pt2("delta_rho_one_point") -s.unset_openmp() -print s - -s = H_apply("pt2_mono_di_delta_rho") -s.set_perturbation("delta_rho_one_point") -s.unset_openmp() -print s - END_SHELL diff --git a/plugins/Full_CI/NEEDED_CHILDREN_MODULES b/plugins/Full_CI/NEEDED_CHILDREN_MODULES index 04ce9e78..ad5f053f 100644 --- a/plugins/Full_CI/NEEDED_CHILDREN_MODULES +++ b/plugins/Full_CI/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full +Perturbation Selectors_full Generators_full Davidson diff --git a/plugins/Perturbation/NEEDED_CHILDREN_MODULES b/plugins/Perturbation/NEEDED_CHILDREN_MODULES index e29a6721..eba3650e 100644 --- a/plugins/Perturbation/NEEDED_CHILDREN_MODULES +++ b/plugins/Perturbation/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Properties Hartree_Fock +Properties Hartree_Fock Davidson diff --git a/plugins/Perturbation/delta_rho_perturbation.irp.f b/plugins/Perturbation/delta_rho_perturbation.irp.f deleted file mode 100644 index c95972a6..00000000 --- a/plugins/Perturbation/delta_rho_perturbation.irp.f +++ /dev/null @@ -1,77 +0,0 @@ -subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,n_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st) - double precision :: i_O1_psi_array(N_st) - double precision :: i_H_psi_array(N_st) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant - ! - ! for the various n_st states, at various level of theory. - ! - ! c_pert(i) = /( - ) - ! - ! e_2_pert(i) = c_pert(i) * - ! - ! H_pert_diag(i) = c_pert(i)^2 * - ! - ! To get the contribution of the first order : - ! - ! = sum(over i) e_2_pert(i) - ! - ! To get the contribution of the diagonal elements of the second order : - ! - ! [ + + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2] - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem,diag_o1_mat_elem_alpha_beta - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase,delta_e,h,oii,diag_o1_mat_elem - integer :: h1,h2,p1,p2,s1,s2 - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - -! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int) -! if(degree.gt.degree_max_generators+1)then -! H_pert_diag = 0.d0 -! e_2_pert = 0.d0 -! c_pert = 0.d0 -! return -! endif - call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array) - - !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(det_pert,Nint) - oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int) - - - do i =1,N_st - if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then - c_pert(i) = -1.d0 - e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) - e_2_pert(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) ) + c_pert(i) * c_pert(i) * oii - H_pert_diag(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) ) - else - c_pert(i) = -1.d0 - e_2_pert(i) = -dabs(i_H_psi_array(i)) - H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i) - endif - enddo - - -end - diff --git a/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES index 7e790003..107c1643 100644 --- a/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES +++ b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Psiref_Utils +Psiref_Utils Davidson diff --git a/plugins/QmcChem/NEEDED_CHILDREN_MODULES b/plugins/QmcChem/NEEDED_CHILDREN_MODULES index aae89501..34de8ddb 100644 --- a/plugins/QmcChem/NEEDED_CHILDREN_MODULES +++ b/plugins/QmcChem/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants +Determinants Davidson diff --git a/src/Davidson/davidson.irp.f b/src/Davidson/davidson.irp.f new file mode 100644 index 00000000..abe3b504 --- /dev/null +++ b/src/Davidson/davidson.irp.f @@ -0,0 +1,3 @@ +program davidson + stop 1 +end diff --git a/src/Determinants/davidson.irp.f b/src/Davidson/diagonalization.irp.f similarity index 89% rename from src/Determinants/davidson.irp.f rename to src/Davidson/diagonalization.irp.f index 2d7b32bf..8ece1c5c 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Davidson/diagonalization.irp.f @@ -1,20 +1,3 @@ -BEGIN_PROVIDER [ integer, davidson_iter_max ] - implicit none - BEGIN_DOC - ! Max number of Davidson iterations - END_DOC - davidson_iter_max = 100 -END_PROVIDER - -BEGIN_PROVIDER [ integer, davidson_sze_max ] - implicit none - BEGIN_DOC - ! Max number of Davidson sizes - END_DOC - ASSERT (davidson_sze_max <= davidson_iter_max) - davidson_sze_max = max(8,2*N_states_diag) -END_PROVIDER - subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) use bitmasks implicit none @@ -69,6 +52,9 @@ end logical function det_inf(key1, key2, Nint) use bitmasks implicit none + BEGIN_DOC +! Ordering function for determinants + END_DOC integer,intent(in) :: Nint integer(bit_kind),intent(in) :: key1(Nint, 2), key2(Nint, 2) integer :: i,j @@ -91,7 +77,6 @@ end function subroutine tamiser(key, idx, no, n, Nint, N_key) use bitmasks implicit none - BEGIN_DOC ! Uncodumented : TODO END_DOC @@ -619,47 +604,3 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ) end -BEGIN_PROVIDER [ character(64), davidson_criterion ] - implicit none - BEGIN_DOC - ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] - END_DOC - davidson_criterion = 'residual' -END_PROVIDER - -subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) - implicit none - BEGIN_DOC -! True if the Davidson algorithm is converged - END_DOC - integer, intent(in) :: N_st, iterations - logical, intent(out) :: converged - double precision, intent(in) :: energy(N_st), residual(N_st) - double precision, intent(in) :: wall, cpu - double precision :: E(N_st), time - double precision, allocatable, save :: energy_old(:) - - if (.not.allocated(energy_old)) then - allocate(energy_old(N_st)) - energy_old = 0.d0 - endif - - E = energy - energy_old - energy_old = energy - if (davidson_criterion == 'energy') then - converged = dabs(maxval(E(1:N_st))) < threshold_davidson - else if (davidson_criterion == 'residual') then - converged = dabs(maxval(residual(1:N_st))) < threshold_davidson - else if (davidson_criterion == 'both') then - converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) & - < threshold_davidson - else if (davidson_criterion == 'wall_time') then - call wall_time(time) - converged = time - wall > threshold_davidson - else if (davidson_criterion == 'cpu_time') then - call cpu_time(time) - converged = time - cpu > threshold_davidson - else if (davidson_criterion == 'iterations') then - converged = iterations >= int(threshold_davidson) - endif -end diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f similarity index 96% rename from src/Determinants/diagonalize_CI.irp.f rename to src/Davidson/diagonalize_CI.irp.f index 58118461..07001902 100644 --- a/src/Determinants/diagonalize_CI.irp.f +++ b/src/Davidson/diagonalize_CI.irp.f @@ -1,19 +1,3 @@ -BEGIN_PROVIDER [ character*(64), diag_algorithm ] - implicit none - BEGIN_DOC - ! Diagonalization algorithm (Davidson or Lapack) - END_DOC - if (N_det > N_det_max_jacobi) then - diag_algorithm = "Davidson" - else - diag_algorithm = "Lapack" - endif - - if (N_det < N_states_diag) then - diag_algorithm = "Lapack" - endif - -END_PROVIDER BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] implicit none diff --git a/src/Determinants/diagonalize_CI_SC2.irp.f b/src/Davidson/diagonalize_CI_SC2.irp.f similarity index 100% rename from src/Determinants/diagonalize_CI_SC2.irp.f rename to src/Davidson/diagonalize_CI_SC2.irp.f diff --git a/src/Determinants/diagonalize_CI_mono.irp.f b/src/Davidson/diagonalize_CI_mono.irp.f similarity index 100% rename from src/Determinants/diagonalize_CI_mono.irp.f rename to src/Davidson/diagonalize_CI_mono.irp.f diff --git a/src/Davidson/parameters.irp.f b/src/Davidson/parameters.irp.f new file mode 100644 index 00000000..66c500c3 --- /dev/null +++ b/src/Davidson/parameters.irp.f @@ -0,0 +1,62 @@ +BEGIN_PROVIDER [ integer, davidson_iter_max ] + implicit none + BEGIN_DOC + ! Max number of Davidson iterations + END_DOC + davidson_iter_max = 100 +END_PROVIDER + +BEGIN_PROVIDER [ integer, davidson_sze_max ] + implicit none + BEGIN_DOC + ! Max number of Davidson sizes + END_DOC + ASSERT (davidson_sze_max <= davidson_iter_max) + davidson_sze_max = max(8,2*N_states_diag) +END_PROVIDER + + +BEGIN_PROVIDER [ character(64), davidson_criterion ] + implicit none + BEGIN_DOC + ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] + END_DOC + davidson_criterion = 'residual' +END_PROVIDER + +subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) + implicit none + BEGIN_DOC +! True if the Davidson algorithm is converged + END_DOC + integer, intent(in) :: N_st, iterations + logical, intent(out) :: converged + double precision, intent(in) :: energy(N_st), residual(N_st) + double precision, intent(in) :: wall, cpu + double precision :: E(N_st), time + double precision, allocatable, save :: energy_old(:) + + if (.not.allocated(energy_old)) then + allocate(energy_old(N_st)) + energy_old = 0.d0 + endif + + E = energy - energy_old + energy_old = energy + if (davidson_criterion == 'energy') then + converged = dabs(maxval(E(1:N_st))) < threshold_davidson + else if (davidson_criterion == 'residual') then + converged = dabs(maxval(residual(1:N_st))) < threshold_davidson + else if (davidson_criterion == 'both') then + converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) & + < threshold_davidson + else if (davidson_criterion == 'wall_time') then + call wall_time(time) + converged = time - wall > threshold_davidson + else if (davidson_criterion == 'cpu_time') then + call cpu_time(time) + converged = time - cpu > threshold_davidson + else if (davidson_criterion == 'iterations') then + converged = iterations >= int(threshold_davidson) + endif +end diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f new file mode 100644 index 00000000..0fd7e52e --- /dev/null +++ b/src/Davidson/u0Hu0.irp.f @@ -0,0 +1,190 @@ +subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes e_0 = / + ! + ! n : number of determinants + ! + END_DOC + integer, intent(in) :: n,Nint + double precision, intent(out) :: e_0 + double precision, intent(in) :: u_0(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + call u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,1,n) +end + +subroutine u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8) + use bitmasks + implicit none + BEGIN_DOC + ! Computes e_0 = / + ! + ! n : number of determinants + ! + END_DOC + integer, intent(in) :: n,Nint, N_st, sze_8 + double precision, intent(out) :: e_0(N_st) + double precision, intent(in) :: u_0(sze_8,N_st) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + + double precision, allocatable :: H_jj(:), v_0(:,:) + double precision :: u_dot_u,u_dot_v,diag_H_mat_elem + integer :: i,j + allocate (H_jj(n), v_0(sze_8,N_st)) + do i = 1, n + H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint) + enddo + + call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) + do i=1,N_st + e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) + enddo +end + + +subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + END_DOC + integer, intent(in) :: n,Nint + double precision, intent(out) :: v_0(n) + double precision, intent(in) :: u_0(n) + double precision, intent(in) :: H_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,1,n) +end + +subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + END_DOC + integer, intent(in) :: N_st,n,Nint, sze_8 + double precision, intent(out) :: v_0(sze_8,N_st) + double precision, intent(in) :: u_0(sze_8,N_st) + double precision, intent(in) :: H_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + double precision :: hij + double precision, allocatable :: vt(:,:) + integer :: i,j,k,l, jj,ii + integer :: i0, j0 + + integer, allocatable :: shortcut(:,:), sort_idx(:,:) + integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) + integer(bit_kind) :: sorted_i(Nint) + + integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate + + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (n>0) + PROVIDE ref_bitmask_energy + + allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) + v_0 = 0.d0 + + call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) + call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8) + allocate(vt(sze_8,N_st)) + Vt = 0.d0 + + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,1) + do sh2=sh,shortcut(0,1) + exa = 0 + do ni=1,Nint + exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) + end do + if(exa > 2) then + cycle + end if + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + if(sh==sh2) then + endi = i-1 + else + endi = shortcut(sh2+1,1)-1 + end if + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh2,1),endi + org_j = sort_idx(j,1) + ext = exa + do ni=1,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + end do + if(ext <= 4) then + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + do istate=1,N_st + vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) + vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) + enddo + endif + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,2) + do i=shortcut(sh,2),shortcut(sh+1,2)-1 + org_i = sort_idx(i,2) + do j=shortcut(sh,2),i-1 + org_j = sort_idx(j,2) + ext = 0 + do ni=1,Nint + ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) + end do + if(ext == 4) then + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + do istate=1,N_st + vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) + vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) + enddo + end if + end do + end do + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do istate=1,N_st + do i=n,1,-1 + v_0(i,istate) = v_0(i,istate) + vt(i,istate) + enddo + enddo + !$OMP END CRITICAL + + deallocate(vt) + !$OMP END PARALLEL + + do istate=1,N_st + do i=1,n + v_0(i,istate) += H_jj(i) * u_0(i,istate) + enddo + enddo + deallocate (shortcut, sort_idx, sorted, version) +end + + diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index d0827f3b..0e62d524 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -1,5 +1,22 @@ use bitmasks +BEGIN_PROVIDER [ character*(64), diag_algorithm ] + implicit none + BEGIN_DOC + ! Diagonalization algorithm (Davidson or Lapack) + END_DOC + if (N_det > N_det_max_jacobi) then + diag_algorithm = "Davidson" + else + diag_algorithm = "Lapack" + endif + + if (N_det < N_states_diag) then + diag_algorithm = "Lapack" + endif +END_PROVIDER + + BEGIN_PROVIDER [ integer, N_det ] implicit none BEGIN_DOC @@ -872,4 +889,4 @@ subroutine apply_hole(det, s1, h1, res, ok, Nint) res(ii, s1) = ibclr(res(ii, s1), pos) ok = .true. -end subroutine \ No newline at end of file +end subroutine diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index c66bc714..3d3eda0b 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -162,7 +162,7 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8) ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) - PROVIDE ref_bitmask_energy davidson_criterion + PROVIDE ref_bitmask_energy allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) v_0 = 0.d0 diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 5e27253a..0d4a1567 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1634,195 +1634,6 @@ subroutine get_occ_from_key(key,occ,Nint) end -subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint) - use bitmasks - implicit none - BEGIN_DOC - ! Computes e_0 = / - ! - ! n : number of determinants - ! - END_DOC - integer, intent(in) :: n,Nint - double precision, intent(out) :: e_0 - double precision, intent(in) :: u_0(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - call u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,1,n) -end - -subroutine u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8) - use bitmasks - implicit none - BEGIN_DOC - ! Computes e_0 = / - ! - ! n : number of determinants - ! - END_DOC - integer, intent(in) :: n,Nint, N_st, sze_8 - double precision, intent(out) :: e_0(N_st) - double precision, intent(in) :: u_0(sze_8,N_st) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - - double precision, allocatable :: H_jj(:), v_0(:,:) - double precision :: u_dot_u,u_dot_v,diag_H_mat_elem - integer :: i,j - allocate (H_jj(n), v_0(sze_8,N_st)) - do i = 1, n - H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint) - enddo - - call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) - do i=1,N_st - e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) - enddo -end - - -subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) - use bitmasks - implicit none - BEGIN_DOC - ! Computes v_0 = H|u_0> - ! - ! n : number of determinants - ! - ! H_jj : array of - END_DOC - integer, intent(in) :: n,Nint - double precision, intent(out) :: v_0(n) - double precision, intent(in) :: u_0(n) - double precision, intent(in) :: H_jj(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,1,n) -end - -subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) - use bitmasks - implicit none - BEGIN_DOC - ! Computes v_0 = H|u_0> - ! - ! n : number of determinants - ! - ! H_jj : array of - END_DOC - integer, intent(in) :: N_st,n,Nint, sze_8 - double precision, intent(out) :: v_0(sze_8,N_st) - double precision, intent(in) :: u_0(sze_8,N_st) - double precision, intent(in) :: H_jj(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - double precision :: hij - double precision, allocatable :: vt(:,:) - integer :: i,j,k,l, jj,ii - integer :: i0, j0 - - integer, allocatable :: shortcut(:,:), sort_idx(:,:) - integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) - integer(bit_kind) :: sorted_i(Nint) - - integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate - - - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - ASSERT (n>0) - PROVIDE ref_bitmask_energy davidson_criterion - - allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - v_0 = 0.d0 - - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) - call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8) - allocate(vt(sze_8,N_st)) - Vt = 0.d0 - - !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0,1) - do sh2=sh,shortcut(0,1) - exa = 0 - do ni=1,Nint - exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) - end do - if(exa > 2) then - cycle - end if - - do i=shortcut(sh,1),shortcut(sh+1,1)-1 - org_i = sort_idx(i,1) - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1,1)-1 - end if - do ni=1,Nint - sorted_i(ni) = sorted(ni,i,1) - enddo - - do j=shortcut(sh2,1),endi - org_j = sort_idx(j,1) - ext = exa - do ni=1,Nint - ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) - end do - if(ext <= 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) - enddo - endif - enddo - enddo - enddo - enddo - !$OMP END DO NOWAIT - - !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0,2) - do i=shortcut(sh,2),shortcut(sh+1,2)-1 - org_i = sort_idx(i,2) - do j=shortcut(sh,2),i-1 - org_j = sort_idx(j,2) - ext = 0 - do ni=1,Nint - ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) - end do - if(ext == 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) - enddo - end if - end do - end do - enddo - !$OMP END DO NOWAIT - - !$OMP CRITICAL - do istate=1,N_st - do i=n,1,-1 - v_0(i,istate) = v_0(i,istate) + vt(i,istate) - enddo - enddo - !$OMP END CRITICAL - - deallocate(vt) - !$OMP END PARALLEL - - do istate=1,N_st - do i=1,n - v_0(i,istate) += H_jj(i) * u_0(i,istate) - enddo - enddo - deallocate (shortcut, sort_idx, sorted, version) -end - subroutine get_double_excitation_phase(det1,det2,exc,phase,Nint) use bitmasks