From c604b34a0b8483f13d6745210aeedfe508e893db Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Feb 2021 18:40:29 +0100 Subject: [PATCH 1/9] Renormalized selection for CSF --- src/cipsi/selection.irp.f | 53 +++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5afb514e..799fce95 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -668,10 +668,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d w = 0d0 -! integer(bit_kind) :: occ(N_int,2), n -! call configuration_of_det(det,occ,N_int) -! call configuration_to_dets_size(occ,n,elec_alpha_num,N_int) - e_pert = 0.d0 coef = 0.d0 logical :: do_diag @@ -699,7 +695,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision :: eigvalues(N_states+1) double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2) - integer :: iwork(3+5*(N_states+1)), info, k ,n + integer :: iwork(3+5*(N_states+1)), info, k if (do_diag) then double precision :: pt2_matrix(N_states+1,N_states+1) @@ -765,36 +761,43 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case(5) ! Variance selection -! w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate) - w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) -! do jstate=1,N_states -! if (istate == jstate) cycle -! w = w + dabs(alpha_h_psi*mat(jstate,p1,p2)) * s_weight(istate,jstate) -! enddo + if (h0_type == 'CFG') then + w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) & + / c0_weight(istate) + else + w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) + endif case(6) -! w = w - coef(istate) * coef(istate) * s_weight(istate,istate) - w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate)) -! do jstate=1,N_states -! if (istate == jstate) cycle -! w = w + dabs(coef(istate)*coef(jstate)) * s_weight(istate,jstate) -! enddo + if (h0_type == 'CFG') then + w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate)) & + / c0_weight(istate) + else + w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate)) + endif case default ! Energy selection -! w = w + e_pert(istate) * s_weight(istate,istate) - w = min(w, e_pert(istate) * s_weight(istate,istate)) -! do jstate=1,N_states -! if (istate == jstate) cycle -! w = w + dabs(X(istate)*X(jstate)) * s_weight(istate,jstate) -! enddo + if (h0_type == 'CFG') then + w = min(w, e_pert(istate) * s_weight(istate,istate)) / c0_weight(istate) + else + w = min(w, e_pert(istate) * s_weight(istate,istate)) + endif end select end do -! w = dble(n) * w - + integer(bit_kind) :: occ(N_int,2), n + if (h0_type == 'CFG') then + do k=1,N_int + occ(k,1) = ieor(det(k,1),det(k,2)) + occ(k,2) = iand(det(k,1),det(k,2)) + enddo + call configuration_to_dets_size(occ,n,elec_alpha_num,N_int) + n = max(n,1) + w *= dsqrt(dble(n)) + endif if(w <= buf%mini) then call add_to_selection_buffer(buf, det, w) From 9f26e533518530e9d7008f15913b60b8a41acd2a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 12 Mar 2021 11:41:58 +0100 Subject: [PATCH 2/9] Fixes #148 --- src/cipsi/pt2_stoch_routines.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 489ffaaf..a0107513 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -187,7 +187,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then stop 'Unable to put pt2_stoch_istate on ZMQ server' endif - if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then stop 'Unable to put threshold_generators on ZMQ server' endif From 2a264766134460ca3849f91470f64ecdd9699c45 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 12 Mar 2021 11:42:15 +0100 Subject: [PATCH 3/9] Change in weights --- src/cipsi/selection.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 799fce95..0b7e99c1 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -796,7 +796,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d enddo call configuration_to_dets_size(occ,n,elec_alpha_num,N_int) n = max(n,1) - w *= dsqrt(dble(n)) + w *= dble(n) endif if(w <= buf%mini) then From d217a6b9f1d897aef2375e3cfb9ef7fbc780733b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 00:26:31 +0100 Subject: [PATCH 4/9] Optimized get_all_spin_singles --- src/determinants/spindeterminants.irp.f | 46 +++++++++---------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 8ba405e0..765ee394 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -1193,16 +1193,10 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo if ( degree == 4 ) then @@ -1247,23 +1241,19 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 2).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo - if ( degree == 2 ) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 + if ( degree /= 2 ) then + cycle endif + singles(n_singles) = idx(i) + n_singles = n_singles+1 + enddo n_singles = n_singles-1 @@ -1296,23 +1286,19 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 + if ( degree /= 4 ) then + cycle endif + doubles(n_doubles) = idx(i) + n_doubles = n_doubles+1 + enddo n_doubles = n_doubles-1 From d79724b7e77357fa075d8da83678f1b511733d7d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 01:04:23 +0100 Subject: [PATCH 5/9] Optimized phasemask --- src/cipsi/selection.irp.f | 77 ++++++++++++++++++++++++++++----------- 1 file changed, 56 insertions(+), 21 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5afb514e..9b92dc9d 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -8,27 +8,56 @@ subroutine get_mask_phase(det1, pm, Nint) integer(bit_kind), intent(out) :: pm(Nint,2) integer(bit_kind) :: tmp1, tmp2 integer :: i - pm(1:Nint,1:2) = det1(1:Nint,1:2) tmp1 = 0_8 tmp2 = 0_8 - do i=1,Nint - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do + select case (Nint) + +BEGIN_TEMPLATE + case ($Nint) + do i=1,$Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +SUBST [ Nint ] +1;; +2;; +3;; +4;; +END_TEMPLATE + case default + do i=1,Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do + end select end subroutine @@ -445,11 +474,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d endif do i=1,fullinteresting(0) - fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i)) + do k=1,N_int + fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) + fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) + enddo enddo do i=1,interesting(0) - minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i)) + do k=1,N_int + minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) + minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) + enddo enddo do s2=s1,2 From 882d8ee3509a0fe522630b882ad8fa9af9ff7168 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 01:04:23 +0100 Subject: [PATCH 6/9] Improved print_energy --- src/cipsi/selection.irp.f | 77 ++++++++++++++++++++++++++---------- src/tools/print_energy.irp.f | 21 +--------- 2 files changed, 57 insertions(+), 41 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5afb514e..9b92dc9d 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -8,27 +8,56 @@ subroutine get_mask_phase(det1, pm, Nint) integer(bit_kind), intent(out) :: pm(Nint,2) integer(bit_kind) :: tmp1, tmp2 integer :: i - pm(1:Nint,1:2) = det1(1:Nint,1:2) tmp1 = 0_8 tmp2 = 0_8 - do i=1,Nint - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do + select case (Nint) + +BEGIN_TEMPLATE + case ($Nint) + do i=1,$Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +SUBST [ Nint ] +1;; +2;; +3;; +4;; +END_TEMPLATE + case default + do i=1,Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do + end select end subroutine @@ -445,11 +474,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d endif do i=1,fullinteresting(0) - fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i)) + do k=1,N_int + fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) + fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) + enddo enddo do i=1,interesting(0) - minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i)) + do k=1,N_int + minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) + minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) + enddo enddo do s2=s1,2 diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f index f78dffc8..4fe1572c 100644 --- a/src/tools/print_energy.irp.f +++ b/src/tools/print_energy.irp.f @@ -14,24 +14,5 @@ end subroutine run implicit none - integer :: i,j - double precision :: i_H_psi_array(N_states) - double precision :: E(N_states) - double precision :: norm(N_states) - - E(1:N_states) = nuclear_repulsion - norm(1:N_states) = 0.d0 - do i=1,N_det - call i_H_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, & - size(psi_coef,1), N_states, i_H_psi_array) - do j=1,N_states - norm(j) += psi_coef(i,j)*psi_coef(i,j) - E(j) += i_H_psi_array(j) * psi_coef(i,j) - enddo - enddo - - print *, 'Energy:' - do i=1,N_states - print *, E(i)/norm(i) - enddo + print *, psi_energy + nuclear_repulsion end From 697fbddde64144e4fc5b0eed4465be7d35fb8339 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 31 Mar 2021 01:50:53 +0200 Subject: [PATCH 7/9] Optimized get_all_spin_doubles_1 --- src/davidson/davidson_parallel.irp.f | 6 ++++-- src/determinants/spindeterminants.irp.f | 28 +++++++++++-------------- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 32f89979..10b41f4c 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -438,7 +438,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) ipos=1 do imin=1,N_det,tasksize imax = min(N_det,imin-1+tasksize) - if (imin==1) then + if (imin<=N_det/2) then istep = 2 else istep = 1 @@ -505,7 +505,9 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) print *, irp_here, ': Failed in zmq_set_running' endif - call omp_set_max_active_levels(4) + + call omp_set_max_active_levels(5) + !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() if (ithread == 0 ) then diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 765ee394..dea4a566 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -1075,19 +1075,17 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, integer :: i include 'utils/constants.include.F' integer :: degree - + integer :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /) + integer :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /) n_singles = 1 n_doubles = 1 do i=1,size_buffer degree = popcnt( xor( spindet, buffer(i) ) ) - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 - else if ( degree == 2 ) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 - endif + doubles(n_doubles) = idx(i) + singles(n_singles) = idx(i) + n_doubles = n_doubles+add_double(degree) + n_singles = n_singles+add_single(degree) enddo n_singles = n_singles-1 n_doubles = n_doubles-1 @@ -1113,15 +1111,14 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_ integer :: i integer(bit_kind) :: v integer :: degree + integer :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /) include 'utils/constants.include.F' n_singles = 1 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - if (degree == 2) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 - endif + singles(n_singles) = idx(i) + n_singles = n_singles+add_single(degree) enddo n_singles = n_singles-1 @@ -1145,14 +1142,13 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_ integer :: i include 'utils/constants.include.F' integer :: degree + integer :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /) n_doubles = 1 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 - endif + doubles(n_doubles) = idx(i) + n_doubles = n_doubles+add_double(degree) enddo n_doubles = n_doubles-1 From 8102d43aff65c0c23ff8d16e5c407516cd46dcc8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 31 Mar 2021 13:33:06 +0200 Subject: [PATCH 8/9] Removed internal reads in zmq --- .../integrals_in_map_slave.irp.f | 2 +- src/cipsi/run_selection_slave.irp.f | 2 +- src/utils/c_funcs.c | 34 +++++ src/utils/c_functions.f90 | 127 +++++++++++++++++- src/zmq/utils.irp.f | 80 ++++++----- 5 files changed, 195 insertions(+), 50 deletions(-) diff --git a/src/ao_two_e_ints/integrals_in_map_slave.irp.f b/src/ao_two_e_ints/integrals_in_map_slave.irp.f index a91bdecb..122fa2ac 100644 --- a/src/ao_two_e_ints/integrals_in_map_slave.irp.f +++ b/src/ao_two_e_ints/integrals_in_map_slave.irp.f @@ -116,7 +116,7 @@ subroutine ao_two_e_integrals_in_map_slave(thread,iproc) exit endif if (task_id == 0) exit - read(task,*) j, l + call sscanf_dd(task, j, l) integer, external :: task_done_to_taskserver call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) == -1) then diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index c2ba2379..91bd3a38 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -52,7 +52,7 @@ subroutine run_selection_slave(thread,iproc,energy) ctask = ctask - 1 else integer :: i_generator, N, subset, bsize - read(task,*) subset, i_generator, N + call sscanf_ddd(task, subset, i_generator, N) if(buf%N == 0) then ! Only first time call create_selection_buffer(N, N*2, buf) diff --git a/src/utils/c_funcs.c b/src/utils/c_funcs.c index 5b7035fb..16b53256 100644 --- a/src/utils/c_funcs.c +++ b/src/utils/c_funcs.c @@ -1,6 +1,40 @@ #include +#include +#include void usleep_c(int s) { usleep((useconds_t) s); } + +void sscanf_ssds_c(const char* str, char* s1, char* s2, int* i, char* s3) +{ + sscanf(str, "%s %s %d %s", s1, s2, i, s3); + s1[strlen(s1)] = ' '; + s2[strlen(s2)] = ' '; + s3[strlen(s3)] = ' '; +} + +void sscanf_dd_c(const char* str, int* i1, int* i2) +{ + sscanf(str, "%d %d", i1, i2); +} + +void sscanf_ddd_c(const char* str, int* i1, int* i2, int* i3) +{ + sscanf(str, "%d %d %d", i1, i2, i3); +} + +void sscanf_ss_c(const char* str, char* s1, char* s2) +{ + sscanf(str, "%s %s", s1, s2); + s1[strlen(s1)] = ' '; + s2[strlen(s2)] = ' '; +} + +void sscanf_sd_c(const char* str, char* s1, int* i) +{ + sscanf(str, "%s %d", s1, i); + s1[strlen(s1)] = ' '; +} + diff --git a/src/utils/c_functions.f90 b/src/utils/c_functions.f90 index 425aafd6..65d4ad62 100644 --- a/src/utils/c_functions.f90 +++ b/src/utils/c_functions.f90 @@ -2,20 +2,133 @@ module c_functions use iso_c_binding interface - subroutine usleep_c(us) bind (C,name="usleep_c") - use iso_c_binding - integer(c_int), value :: us - end subroutine usleep_c + subroutine usleep_c(us) bind (C,name="usleep_c") + use iso_c_binding + integer(c_int), value :: us + end subroutine usleep_c end interface -end module + interface + integer(c_int) function atoi_c(a) bind (C,name="atoi") + use iso_c_binding + character(kind=c_char), intent(in) :: a(*) + end function atoi_c + end interface -subroutine usleep(us) + interface + subroutine sscanf_ss_c(str,s1, s2) bind (C) + use iso_c_binding + character(kind=c_char), intent(in ) :: str(*) + character(kind=c_char), intent(out) :: s1(*),s2(*) + end subroutine sscanf_ss_c + end interface + + interface + subroutine sscanf_ssds_c(str, s1, s2, i, s3) bind (C) + use iso_c_binding + character(kind=c_char), intent(in ) :: str(*) + character(kind=c_char), intent(out) :: s1(*),s2(*),s3(*) + integer(kind=c_int) , intent(out) :: i + end subroutine sscanf_ssds_c + end interface + + interface + subroutine sscanf_dd_c(str, i1, i2) bind (C) + use iso_c_binding + character(kind=c_char), intent(in ) :: str(*) + integer(kind=c_int) , intent(out) :: i1, i2 + end subroutine sscanf_dd_c + end interface + + interface + subroutine sscanf_ddd_c(str, i1, i2, i3) bind (C) + use iso_c_binding + character(kind=c_char), intent(in ) :: str(*) + integer(kind=c_int) , intent(out) :: i1, i2, i3 + end subroutine sscanf_ddd_c + end interface + + interface + subroutine sscanf_sd_c(str,s1, i) bind (C) + use iso_c_binding + character(kind=c_char), intent(in ) :: str(*) + character(kind=c_char), intent(out) :: s1(*) + integer(kind=c_int) , intent(out) :: i + end subroutine sscanf_sd_c + end interface + +contains + + integer function atoi(a) + implicit none + character(len=*), intent(in) :: a + atoi = atoi_c(trim(a)//c_null_char) + end function atoi + +end module c_functions + +subroutine sscanf_ss(str, s1,s2) use c_functions use iso_c_binding implicit none + character(*), intent(in) :: str + character(*), intent(out) :: s1,s2 + s1 = ' ' + s2 = ' ' + call sscanf_ss_c(trim(str)//c_null_char, s1, s2) +end subroutine sscanf_ss + +subroutine sscanf_sd(str, s1,i) + use c_functions + use iso_c_binding + implicit none + character(*), intent(in) :: str + character(*), intent(out) :: s1 + integer, intent(out) :: i + s1 = ' ' + call sscanf_sd_c(trim(str)//c_null_char, s1, i) +end subroutine sscanf_sd + +subroutine sscanf_ssds(str, s1,s2,i,s3) + use c_functions + use iso_c_binding + implicit none + character(*), intent(in) :: str + character(*), intent(out) :: s1,s2,s3 + integer, intent(out) :: i + s1 = ' ' + s2 = ' ' + s3 = ' ' + call sscanf_ssds_c(trim(str)//c_null_char, s1, s2, i, s3) +end subroutine sscanf_ssds + +subroutine sscanf_dd(str, i1,i2) + use c_functions + use iso_c_binding + implicit none + character(*), intent(in) :: str + integer, intent(out) :: i1, i2 + call sscanf_dd_c(trim(str)//c_null_char, i1, i2) +end subroutine sscanf_dd + +subroutine sscanf_ddd(str, i1,i2,i3) + use c_functions + use iso_c_binding + implicit none + character(*), intent(in) :: str + integer, intent(out) :: i1, i2, i3 + call sscanf_ddd_c(trim(str)//c_null_char, i1, i2, i3) +end subroutine sscanf_ddd + + +subroutine usleep(us) + use iso_c_binding + use c_functions + implicit none integer, intent(in) :: us integer(c_int) :: u u = us call usleep_c(u) -end +end subroutine usleep + + diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index 058f6bca..93dbd16a 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -16,6 +16,7 @@ END_PROVIDER BEGIN_PROVIDER [ character*(128), qp_run_address ] &BEGIN_PROVIDER [ integer, zmq_port_start ] use f77_zmq + use c_functions implicit none BEGIN_DOC ! Address of the qp_run socket @@ -32,14 +33,15 @@ END_PROVIDER do i=len(buffer),1,-1 if ( buffer(i:i) == ':') then qp_run_address = trim(buffer(1:i-1)) - read(buffer(i+1:), *, err=10,end=10) zmq_port_start + zmq_port_start = atoi(buffer(i+1:)) exit endif enddo - return - 10 continue - print *, irp_here, ': Error in read' - stop -1 + + if (zmq_port_start == 0) then + print *, irp_here, ': zmq_port_start is 0' + stop -1 + endif END_PROVIDER BEGIN_PROVIDER [ character*(128), zmq_socket_pull_tcp_address ] @@ -84,6 +86,7 @@ end subroutine switch_qp_run_to_master use f77_zmq + use c_functions implicit none BEGIN_DOC ! Address of the master qp_run socket @@ -102,16 +105,17 @@ subroutine switch_qp_run_to_master do i=len(buffer),1,-1 if ( buffer(i:i) == ':') then qp_run_address = trim(buffer(1:i-1)) - read(buffer(i+1:), *, end=10, err=10) zmq_port_start + zmq_port_start = atoi(buffer(i+1:)) exit endif enddo call reset_zmq_addresses return - 10 continue - print *, irp_here, ': Error in read' - stop -1 + if (zmq_port_start == 0) then + print *, irp_here, ': zmq_port_start is 0' + stop -1 + endif end @@ -650,12 +654,17 @@ integer function connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0) message = trim(message(1:rc)) if(message(1:5) == "error") then - go to 10 + connect_to_taskserver = -1 + return end if - read(message,*, end=10, err=10) reply, state, worker_id, address + + call sscanf_ssds(message, reply, state, worker_id, address) + if (trim(reply) /= 'connect_reply') then - go to 10 + connect_to_taskserver = -1 + return endif + if (trim(state) /= zmq_state) then integer, external :: disconnect_from_taskserver_state if (disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_id, state) == -1) then @@ -663,13 +672,8 @@ integer function connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) continue endif connect_to_taskserver = -1 - return endif - return - 10 continue -! print *, irp_here//': '//trim(message) - connect_to_taskserver = -1 end integer function disconnect_from_taskserver(zmq_to_qp_run_socket, worker_id) @@ -698,7 +702,7 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i character*(512) :: message, reply character*(128) :: state_tmp - disconnect_from_taskserver_state = 0 + disconnect_from_taskserver_state = -1 write(message,*) 'disconnect '//trim(state), worker_id @@ -718,21 +722,15 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i rc = min(510,rc) message = trim(message(1:rc)) - read(message,*, end=10, err=10) reply, state_tmp - if ((trim(reply) == 'disconnect_reply').and.(trim(state_tmp) == trim(state))) then - return - endif - if (trim(message) == 'error Wrong state') then - disconnect_from_taskserver_state = -1 - return - else if (trim(message) == 'error No job is running') then - disconnect_from_taskserver_state = -1 + call sscanf_ss(message, reply, state_tmp) + + if (trim(state_tmp) /= trim(state)) then return endif - return - 10 continue - disconnect_from_taskserver_state = -1 + if ((trim(reply) == 'disconnect_reply')) then + disconnect_from_taskserver_state = 0 + endif end integer function add_task_to_taskserver(zmq_to_qp_run_socket,task) @@ -898,7 +896,7 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id character*(1024) :: message character*(64) :: reply - integer :: rc, sze + integer :: rc, sze, i get_task_from_taskserver = 0 @@ -912,12 +910,15 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id endif task_id = 0 - message = repeat(' ',1024) + message = ' ' rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 1024, 0) - rc = min(64,rc) - read(message(1:rc),*, end=10, err=10) reply - if (trim(reply) == 'get_task_reply') then - read(message(1:rc),*, end=10, err=10) reply, task_id + i = 1 + do while (message(i:i) /= ' ') + i = i+1 + enddo + reply = message(1:i-1) + if (reply == 'get_task_reply') then + call sscanf_sd(message, reply, task_id) rc = 15 do while (rc < 1024 .and. message(rc:rc) == ' ') rc += 1 @@ -937,15 +938,12 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id get_task_from_taskserver = -1 return endif - return - - 10 continue - get_task_from_taskserver = -1 end integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id,task,n_tasks) + use c_functions use f77_zmq implicit none BEGIN_DOC @@ -1000,7 +998,7 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i return endif rc = min(1024,rc) - read(message(1:rc),*, end=10, err=10) task_id(i) + task_id(i) = atoi(message(1:rc)) if (task_id(i) == 0) then task(i) = 'terminate' n_tasks = i From 83ae7ff58156af011c45a421395b344d39e1d6e5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 1 Apr 2021 21:47:17 +0200 Subject: [PATCH 9/9] Fixed 2RDM broken in 533413277d --- src/two_body_rdm/two_e_dm_mo.irp.f | 1 + 1 file changed, 1 insertion(+) diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index 47f5983f..107161c8 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -34,6 +34,7 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num,1)] enddo enddo enddo + two_e_dm_mo(:,:,:,:,:) = two_e_dm_mo(:,:,:,:,:) * 2.d0 END_PROVIDER