From 2f246780eb0a08ee87ab5d98fe9c0e2f17685594 Mon Sep 17 00:00:00 2001 From: ydamour Date: Tue, 13 Jun 2023 14:05:13 +0200 Subject: [PATCH 1/3] fix bug in get_excitation_general --- src/utils_cc/org/phase.org | 2 ++ src/utils_cc/phase.irp.f | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/utils_cc/org/phase.org b/src/utils_cc/org/phase.org index 5f67859c..2156a251 100644 --- a/src/utils_cc/org/phase.org +++ b/src/utils_cc/org/phase.org @@ -137,6 +137,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n1(j) + if (k > n_anni(j)) exit if (l1(i,j) /= list_anni(k,j)) cycle pos_anni(k,j) = i k = k + 1 @@ -147,6 +148,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n2(j) + if (k > n_crea(j)) exit if (l2(i,j) /= list_crea(k,j)) cycle pos_crea(k,j) = i k = k + 1 diff --git a/src/utils_cc/phase.irp.f b/src/utils_cc/phase.irp.f index 01b41f49..e0703fb8 100644 --- a/src/utils_cc/phase.irp.f +++ b/src/utils_cc/phase.irp.f @@ -96,6 +96,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n1(j) + if (k > n_anni(j)) exit if (l1(i,j) /= list_anni(k,j)) cycle pos_anni(k,j) = i k = k + 1 @@ -106,6 +107,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N do j = 1, 2 k = 1 do i = 1, n2(j) + if (k > n_crea(j)) exit if (l2(i,j) /= list_crea(k,j)) cycle pos_crea(k,j) = i k = k + 1 From a56644a808e0aea4c16d72c61b717ac4b2e0cabc Mon Sep 17 00:00:00 2001 From: ydamour Date: Tue, 13 Jun 2023 14:24:39 +0200 Subject: [PATCH 2/3] remove old stuffs --- src/mo_optimization/my_providers.irp.f | 141 ------------------------- 1 file changed, 141 deletions(-) delete mode 100644 src/mo_optimization/my_providers.irp.f diff --git a/src/mo_optimization/my_providers.irp.f b/src/mo_optimization/my_providers.irp.f deleted file mode 100644 index 7469ffd5..00000000 --- a/src/mo_optimization/my_providers.irp.f +++ /dev/null @@ -1,141 +0,0 @@ -! Dimensions of MOs - - -BEGIN_PROVIDER [ integer, n_mo_dim ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of MOs we can build, - ! with i>j - END_DOC - - n_mo_dim = mo_num*(mo_num-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_core ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of core MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_core = dim_list_core_orb*(dim_list_core_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_act ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of active MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_act = dim_list_act_orb*(dim_list_act_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_inact ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of inactive MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_inact = dim_list_inact_orb*(dim_list_inact_orb-1)/2 - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_mo_dim_virt ] - implicit none - BEGIN_DOC - ! Number of different pairs (i,j) of virtual MOs we can build, - ! with i>j - END_DOC - - n_mo_dim_virt = dim_list_virt_orb*(dim_list_virt_orb-1)/2 - -END_PROVIDER - -! Energies/criterions - -BEGIN_PROVIDER [ double precision, my_st_av_energy ] - implicit none - BEGIN_DOC - ! State average CI energy - END_DOC - - !call update_st_av_ci_energy(my_st_av_energy) - call state_average_energy(my_st_av_energy) - -END_PROVIDER - -! With all the MOs - -BEGIN_PROVIDER [ double precision, my_gradient_opt, (n_mo_dim) ] -&BEGIN_PROVIDER [ double precision, my_CC1_opt ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, for all the MOs. - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision :: norm_grad - - PROVIDE mo_two_e_integrals_in_map - - call gradient_opt(n_mo_dim, my_gradient_opt, my_CC1_opt, norm_grad) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, my_hessian_opt, (n_mo_dim, n_mo_dim) ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, for all the MOs. - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision, allocatable :: h_f(:,:,:,:) - - PROVIDE mo_two_e_integrals_in_map - - allocate(h_f(mo_num, mo_num, mo_num, mo_num)) - - call hessian_list_opt(n_mo_dim, my_hessian_opt, h_f) - -END_PROVIDER - -! With the list of active MOs -! Can be generalized to any mo_class by changing the list/dimension - -BEGIN_PROVIDER [ double precision, my_gradient_list_opt, (n_mo_dim_act) ] -&BEGIN_PROVIDER [ double precision, my_CC2_opt ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, only for the active MOs ! - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision :: norm_grad - - PROVIDE mo_two_e_integrals_in_map !one_e_dm_mo two_e_dm_mo mo_one_e_integrals - - call gradient_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_list_opt, my_CC2_opt, norm_grad) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, my_hessian_list_opt, (n_mo_dim_act, n_mo_dim_act) ] - implicit none - BEGIN_DOC - ! - Gradient of the energy with respect to the MO rotations, only for the active MOs ! - ! - Maximal element of the gradient in absolute value - END_DOC - - double precision, allocatable :: h_f(:,:,:,:) - - PROVIDE mo_two_e_integrals_in_map - - allocate(h_f(dim_list_act_orb, dim_list_act_orb, dim_list_act_orb, dim_list_act_orb)) - - call hessian_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_list_opt, h_f) - -END_PROVIDER From 88f168724e65038b4d9ce9d4a566252de62f1fb5 Mon Sep 17 00:00:00 2001 From: ydamour Date: Thu, 15 Jun 2023 14:46:17 +0200 Subject: [PATCH 3/3] fix binary search (T) --- src/ccsd/ccsd_t_space_orb_stoch.irp.f | 82 +++++++++++++++------------ 1 file changed, 47 insertions(+), 35 deletions(-) diff --git a/src/ccsd/ccsd_t_space_orb_stoch.irp.f b/src/ccsd/ccsd_t_space_orb_stoch.irp.f index b669025e..31fe67ce 100644 --- a/src/ccsd/ccsd_t_space_orb_stoch.irp.f +++ b/src/ccsd/ccsd_t_space_orb_stoch.irp.f @@ -104,17 +104,17 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ integer*8, allocatable :: sampled(:) ! integer(omp_lock_kind), allocatable :: lock(:) integer*2 , allocatable :: abc(:,:) - integer*8 :: Nabc, i8 + integer*8 :: Nabc, i8,kiter integer*8, allocatable :: iorder(:) double precision :: eocc double precision :: norm - integer :: kiter, isample + integer :: isample ! Prepare table of triplets (a,b,c) Nabc = (int(nV,8) * int(nV+1,8) * int(nV+2,8))/6_8 - nV - allocate (memo(Nabc), sampled(Nabc), Pabc(Nabc), waccu(Nabc)) + allocate (memo(Nabc), sampled(Nabc), Pabc(Nabc), waccu(0:Nabc)) allocate (abc(4,Nabc), iorder(Nabc)) !, lock(Nabc)) ! eocc = 3.d0/dble(nO) * sum(f_o(1:nO)) @@ -124,21 +124,21 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ do c = b+1, nV Nabc = Nabc + 1_8 Pabc(Nabc) = -1.d0/(f_v(a) + f_v(b) + f_v(c)) - abc(1,Nabc) = a - abc(2,Nabc) = b - abc(3,Nabc) = c + abc(1,Nabc) = int(a,2) + abc(2,Nabc) = int(b,2) + abc(3,Nabc) = int(c,2) enddo Nabc = Nabc + 1_8 - abc(1,Nabc) = a - abc(2,Nabc) = b - abc(3,Nabc) = a + abc(1,Nabc) = int(a,2) + abc(2,Nabc) = int(b,2) + abc(3,Nabc) = int(a,2) Pabc(Nabc) = -1.d0/(2.d0*f_v(a) + f_v(b)) Nabc = Nabc + 1_8 - abc(1,Nabc) = b - abc(2,Nabc) = a - abc(3,Nabc) = b + abc(1,Nabc) = int(b,2) + abc(2,Nabc) = int(a,2) + abc(3,Nabc) = int(b,2) Pabc(Nabc) = -1.d0/(f_v(a) + 2.d0*f_v(b)) enddo enddo @@ -169,6 +169,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ waccu(i8) = waccu(i8+1) - Pabc(i8+1) enddo waccu(:) = waccu(:) + 1.d0 + waccu(0) = 0.d0 logical :: converged, do_comp double precision :: eta, variance, error, sample @@ -222,8 +223,12 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ do kiter=1,Nabc !$OMP MASTER - do while ((imin <= Nabc).and.(sampled(imin)>-1_8)) - imin = imin+1 + do while (imin <= Nabc) + if (sampled(imin)>-1_8) then + imin = imin+1 + else + exit + endif enddo ! Deterministic part @@ -301,6 +306,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ endif enddo + isample = min(isample,nbuckets) do ieta=bounds(1,isample), Nabc w = dble(max(sampled(ieta),0_8)) tmp = w * memo(ieta) * Pabc(ieta) @@ -331,33 +337,39 @@ end -integer*8 function binary_search(arr, key, size) +integer*8 function binary_search(arr, key, sze) implicit none BEGIN_DOC -! Searches the key in array arr(1:size) between l_in and r_in, and returns its index +! Searches the key in array arr(1:sze) between l_in and r_in, and returns its index END_DOC - integer*8 :: size, i, j, mid, l_in, r_in - double precision, dimension(size) :: arr(1:size) + integer*8 :: sze, i, j, mid + double precision :: arr(0:sze) double precision :: key - i = 1_8 - j = size + if ( key < arr(1) ) then + binary_search = 0_8 + return + end if - do while (j >= i) - mid = i + (j - i) / 2 - if (arr(mid) >= key) then - if (mid > 1 .and. arr(mid - 1) < key) then - binary_search = mid - return - end if - j = mid - 1 - else if (arr(mid) < key) then - i = mid + 1 - else - binary_search = mid + 1 - return - end if + if ( key >= arr(sze) ) then + binary_search = sze + return + end if + + i = 0_8 + j = sze + 1_8 + + do while (.True.) + mid = (i + j) / 2_8 + if ( key >= arr(mid) ) then + i = mid + else + j = mid + end if + if (j-i <= 1_8) then + binary_search = i + return + endif end do - binary_search = i end function binary_search