diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index bd2d37a3..9d33523b 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -32,19 +32,23 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, if(degree.ne.2)then return endif - - call bitstring_to_list_ab(key_i, occ, Ne, Nint) + integer :: degree_i,degree_j + call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int) + call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int) call get_double_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) if(s1.ne.s2)then ! opposite spin two-body -! key_j, key_i htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) if(three_body_h_tc)then if(.not.double_normal_ord)then - call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) - elseif(double_normal_ord.and.+Ne(1).gt.2)then + if(degree_i>degree_j)then + call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) + else + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + endif + elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? endif endif @@ -56,8 +60,12 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) if(three_body_h_tc)then if(.not.double_normal_ord)then - call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) - elseif(double_normal_ord.and.+Ne(1).gt.2)then + if(degree_i>degree_j)then + call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) + else + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + endif + elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? endif diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 66ca2e6a..99352162 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,13 +11,16 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid -! call test_slater_tc_opt - call timing_hij +! call test_slater_tc_opt + call timing_tot +! call timing_diag +! call timing_single +! call timing_double end subroutine test_slater_tc_opt implicit none - integer :: i,j + integer :: i,j,degree double precision :: hmono, htwoe, htot, hthree double precision :: hnewmono, hnewtwoe, hnewthree, hnewtot double precision :: accu_d ,i_count, accu @@ -25,84 +28,193 @@ subroutine test_slater_tc_opt accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det -! do i = 1,1 - call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) - call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) - call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) -! print*,hthree,hnewthree -! print*,htot,hnewtot,dabs(hnewtot-htot) - accu_d += dabs(htot-hnewtot) - if(dabs(htot-hnewtot).gt.1.d-8)then - print*,i - print*,htot,hnewtot,dabs(htot-hnewtot) - endif -! do j = 319,319 do j = 1,N_det - if(i==j)cycle - integer :: degree - call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) -! if(degree .ne. 1)cycle call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot) -! if(dabs(hthree).gt.1.d-15)then if(dabs(htot).gt.1.d-15)then -! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then i_count += 1.D0 accu += dabs(htot-hnewtot) -! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then - if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i,degree - call debug_det(psi_det(1,1,i),N_int) - call debug_det(psi_det(1,1,j),N_int) - print*,htot,hnewtot,dabs(htot-hnewtot) -! print*,hthree,hnewthree,dabs(hthree-hnewthree) - stop - endif -! print*,htot,hnewtot,dabs(htot-hnewtot) + if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + print*,j,i,degree + call debug_det(psi_det(1,1,i),N_int) + call debug_det(psi_det(1,1,j),N_int) + print*,htot,hnewtot,dabs(htot-hnewtot) + print*,hthree,hnewthree,dabs(hthree-hnewthree) + stop + endif endif enddo enddo - print*,'accu_d = ',accu_d/N_det print*,'accu = ',accu/i_count end - -subroutine timing_hij +subroutine timing_tot implicit none integer :: i,j double precision :: wall0, wall1 double precision, allocatable :: mat_old(:,:),mat_new(:,:) - double precision :: hmono, htwoe, hthree, htot -! allocate(mat_old(N_det,N_det)) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det +! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + i_count += 1.d0 + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for old hij for total = ',wall1 - wall0 + + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det +! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + i_count += 1.d0 + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij for total = ',wall1 - wall0 + call i_H_j(psi_det(1,1,1), psi_det(1,1,2),N_int,htot) + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_H_j(psi_det(1,1,j), psi_det(1,1,i),N_int,htot) + i_count += 1.d0 + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij STANDARD = ',wall1 - wall0 + +end + +subroutine timing_diag + implicit none + integer :: i,j + double precision :: wall0, wall1 + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) call wall_time(wall0) + i_count = 0.d0 do i = 1, N_det - do j = 1, N_det + do j = i,i + i_count += 1.d0 call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) -! mat_old(j,i) = htot enddo enddo call wall_time(wall1) - print*,'time for old hij = ',wall1 - wall0 + print*,'i_count = ',i_count + print*,'time for old hij for diagonal= ',wall1 - wall0 -! allocate(mat_new(N_det,N_det)) call wall_time(wall0) + i_count = 0.d0 do i = 1, N_det - do j = 1, N_det + do j = i,i + i_count += 1.d0 call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) -! mat_new(j,i) = htot enddo enddo call wall_time(wall1) - print*,'time for new hij = ',wall1 - wall0 - double precision :: accu + print*,'i_count = ',i_count + print*,'time for new hij for diagonal= ',wall1 - wall0 + +end + +subroutine timing_single + implicit none + integer :: i,j + double precision :: wall0, wall1,accu + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + i_count = 0.d0 accu = 0.d0 do i = 1, N_det do j = 1, N_det -! accu += dabs(mat_new(j,i) - mat_old(j,i)) + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.1)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 enddo enddo - print*,'accu = ',accu + print*,'i_count = ',i_count + print*,'time for old hij for singles = ',accu + + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.1)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + print*,'i_count = ',i_count + print*,'time for new hij for singles = ',accu end + +subroutine timing_double + implicit none + integer :: i,j + double precision :: wall0, wall1,accu + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.2)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + print*,'i_count = ',i_count + print*,'time for old hij for doubles = ',accu + + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.2)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij for doubles = ',accu + +end +