From 257d7004aa77aedd4bff1e75982bd98a45053df0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Jul 2018 01:13:04 +0200 Subject: [PATCH] OpenMP in save_MRCC_wf --- plugins/CAS_SD_ZMQ/cassd_zmq.irp.f | 76 +++++++++++++++++++- plugins/Full_CI_ZMQ/fci_zmq.irp.f | 4 +- plugins/Generators_CAS/generators.irp.f | 2 +- plugins/dress_zmq/alpha_factory.irp.f | 10 +-- plugins/mrcepa0/dressing.irp.f | 4 +- plugins/mrcepa0/save_mrcc_wavefunction.irp.f | 8 +++ src/Determinants/slater_rules.irp.f | 1 - 7 files changed, 94 insertions(+), 11 deletions(-) diff --git a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f index 44f9afbc..0cef5d6f 100644 --- a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f +++ b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f @@ -5,16 +5,28 @@ program cassd_zmq integer :: degree integer :: n_det_before, to_select double precision :: threshold_davidson_in + double precision :: error(N_states) allocate (pt2(N_states)) double precision :: hf_energy_ref logical :: has + integer :: N_states_p + character*(512) :: fmt + character*(8) :: pt2_string + pt2 = -huge(1.d0) + error = 0.d0 threshold_davidson_in = threshold_davidson threshold_davidson = threshold_davidson_in * 100.d0 SOFT_TOUCH threshold_davidson + if (do_pt2) then + pt2_string = ' ' + else + pt2_string = '(approx)' + endif + call diagonalize_CI call save_wavefunction @@ -45,7 +57,6 @@ program cassd_zmq double precision :: E_CI_before(N_states) - print*,'Beginning the selection ...' if (.True.) then ! Avoid pre-calculation of CI_energy E_CI_before(1:N_states) = CI_energy(1:N_states) endif @@ -60,6 +71,8 @@ program cassd_zmq (maxval(abs(pt2(1:N_states))) > pt2_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) + write(*,'(A)') '--------------------------------------------------------------------------------' + correlation_energy_ratio = (CI_energy(1) - hf_energy_ref) / & (E_CI_before(1) + pt2(1) - hf_energy_ref) @@ -98,6 +111,67 @@ program cassd_zmq to_select = min(to_select, N_det_max-n_det_before) call ZMQ_selection(to_select, pt2) + N_states_p = min(N_det,N_states) + + print *, '' + print '(A,I12)', 'Summary at N_det = ', N_det + print '(A)', '-----------------------------------' + print *, '' + call write_double(6,correlation_energy_ratio, 'Correlation ratio') + print *, '' + + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))' + write(*,fmt) ('State',k, k=1,N_states_p) + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))' + write(*,fmt) '# E ', E_CI_before(1:N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', E_CI_before(1:N_states_p)-E_CI_before(1) + write(*,fmt) '# Excit. (eV)', (E_CI_before(1:N_states_p)-E_CI_before(1))*27.211396641308d0 + endif + write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))' + write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p) + write(*,'(A)') '#' + write(*,fmt) '# E+PT2 ', (E_CI_before(k)+pt2(k),error(k), k=1,N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1)), & + dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1))*27.211396641308d0, & + dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p) + endif + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + print *, '' + + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print*, 'correlation_ratio = ', correlation_energy_ratio + + do k=1, N_states_p + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', E_CI_before(k) + print *, 'E+PT2'//pt2_string//' = ', E_CI_before(k)+pt2(k), ' +/- ', error(k) + enddo + + print *, '-----' + if(N_states.gt.1)then + print *, 'Variational Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (E_CI_before(i) - E_CI_before(1)), & + (E_CI_before(i) - E_CI_before(1)) * 27.211396641308d0 + enddo + print *, '-----' + print*, 'Variational + perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))), & + (E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))) * 27.211396641308d0 + enddo + endif + PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index 2d3c1f29..933056d4 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -22,8 +22,8 @@ program fci_zmq threshold_davidson = threshold_davidson_in * 100.d0 SOFT_TOUCH threshold_davidson - call diagonalize_CI - call save_wavefunction +! call diagonalize_CI +! call save_wavefunction call ezfio_has_hartree_fock_energy(has) if (has) then diff --git a/plugins/Generators_CAS/generators.irp.f b/plugins/Generators_CAS/generators.irp.f index 4be8c061..4d99f3a1 100644 --- a/plugins/Generators_CAS/generators.irp.f +++ b/plugins/Generators_CAS/generators.irp.f @@ -35,10 +35,10 @@ END_PROVIDER integer, external :: number_of_holes,number_of_particles integer, allocatable :: nongen(:) integer :: inongen - inongen = 0 allocate(nongen(N_det)) + inongen = 0 m=0 do i=1,N_det good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 ) diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f index 261966be..f590f5d1 100644 --- a/plugins/dress_zmq/alpha_factory.irp.f +++ b/plugins/dress_zmq/alpha_factory.irp.f @@ -66,7 +66,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index integer(bit_kind), allocatable:: preinteresting_det(:,:,:) integer ,allocatable :: abuf(:), labuf(:) - allocate(abuf(N_det*6), labuf(N_det)) + allocate(abuf(0:N_det*6), labuf(0:N_det)) allocate(preinteresting_det(N_int,2,N_det)) PROVIDE fragment_count @@ -387,7 +387,7 @@ subroutine alpha_callback_mask(delta_ij_loc, i_gen, sp, mask, bannedOrb, banned, integer(bit_kind), allocatable :: det_minilist(:,:,:) - allocate(abuf(siz), labuf(N_det), putten(N_det), det_minilist(N_int, 2, N_det)) + allocate(abuf(0:siz), labuf(0:N_det), putten(N_det), det_minilist(N_int, 2, N_det)) do i=1,siz abuf(i) = psi_from_sorted_gen(rabuf(i)) @@ -638,7 +638,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, indexes, ab integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) logical, intent(inout) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2) integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num) - integer, intent(inout) :: abuf(*) + integer, intent(inout) :: abuf(0:*) integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt, s integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) integer :: phasemask(2,N_int*bit_kind_size) @@ -704,7 +704,7 @@ subroutine get_d2(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp) implicit none integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) - integer, intent(inout) :: abuf(*) + integer, intent(inout) :: abuf(0:*) integer, intent(in) :: i_gen logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num) @@ -832,7 +832,7 @@ subroutine get_d1(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp) implicit none integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) - integer, intent(inout) :: abuf(*) + integer, intent(inout) :: abuf(0:*) integer,intent(in) :: i_gen logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) integer(bit_kind) :: det(N_int, 2) diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 319d0603..f7ee5ace 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -236,6 +236,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b double precision, intent(inout) :: contrib(N_states) double precision :: sdress, hdress + PROVIDE n_act_orb elec_num if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat @@ -1025,7 +1026,7 @@ subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_m if (good) then if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then N_tq += 1 - do k=1,N_int + do k=1,Nint tq(k,1,N_tq) = det_buffer(k,1,i) tq(k,2,N_tq) = det_buffer(k,2,i) enddo @@ -1131,6 +1132,7 @@ subroutine get_cc_coef(tq,c_alpha) integer :: i_state, k_sd, l_sd, i_I logical :: ok + PROVIDE n_act_orb elec_num if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat endif diff --git a/plugins/mrcepa0/save_mrcc_wavefunction.irp.f b/plugins/mrcepa0/save_mrcc_wavefunction.irp.f index 49a38557..9afa2b71 100644 --- a/plugins/mrcepa0/save_mrcc_wavefunction.irp.f +++ b/plugins/mrcepa0/save_mrcc_wavefunction.irp.f @@ -3,6 +3,7 @@ program save_mrcc_wf threshold_generators = 1.d0 threshold_selectors = 1.d0 + PROVIDE N_int psi_det TOUCH threshold_generators threshold_selectors mrmode=5 @@ -21,12 +22,19 @@ subroutine run1 double precision :: c_alpha(N_states) call set_generators_bitmasks_as_holes_and_particles + call get_cc_coef(psi_det(1,1,1), c_alpha) + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP PRIVATE(k,c_alpha) SCHEDULE(static,64) do k=1,N_det if (maxval(abs(psi_coef(k,1:N_states))) == 0.d0) then + if (iand(k,1023) == 0) then + print *, k, '/', N_det + endif call get_cc_coef(psi_det(1,1,k), c_alpha) psi_coef(k,1:N_states) = c_alpha(1:N_states) endif enddo + !$OMP END PARALLEL DO SOFT_TOUCH psi_coef end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index a3120be9..6e972114 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -568,7 +568,6 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) case (0) - print *,irp_here,": ZERO" double precision, external :: diag_S_mat_elem s2 = diag_S_mat_elem(key_i,Nint) hij = diag_H_mat_elem(key_i,Nint)