diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 413a9d9a..060722a1 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -1,18 +1,20 @@ subroutine run_cipsi implicit none + use selection_types BEGIN_DOC ! Selected Full Configuration Interaction with deterministic selection and ! stochastic PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:) + type(pt2_type) :: pt2_data + double precision, allocatable :: rpt2(:), zeros(:) integer :: to_select logical, external :: qp_stop double precision :: threshold_generators_save double precision :: rss double precision, external :: memory_of_double - PROVIDE H_apply_buffer_allocated + PROVIDE H_apply_buffer_allocated N_iter = 1 threshold_generators = 1.d0 @@ -21,7 +23,8 @@ subroutine run_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states)) + allocate (zeros(N_states), rpt2(N_states)) + allocate (pt2_data % pt2(N_states), pt2_data % norm2(N_states), pt2_data % variance(N_states)) double precision :: hf_energy_ref logical :: has @@ -30,10 +33,10 @@ subroutine run_cipsi relative_error=PT2_relative_error zeros = 0.d0 - pt2 = -huge(1.e0) rpt2 = -huge(1.e0) - norm = 0.d0 - variance = huge(1.e0) + pt2_data % pt2 = -huge(1.e0) + pt2_data % norm2 = 0.d0 + pt2_data % variance = huge(1.e0) if (s2_eig) then call make_s2_eigenfunction @@ -67,8 +70,8 @@ subroutine run_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(variance(1:N_states))) > variance_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' @@ -77,22 +80,21 @@ subroutine run_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) if (do_pt2) then - pt2 = 0.d0 - variance = 0.d0 - norm = 0.d0 + pt2_data % pt2 = 0.d0 + pt2_data % variance = 0.d0 + pt2_data % norm2 = 0.d0 threshold_generators_save = threshold_generators threshold_generators = 1.d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm, 0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, 0) ! Stochastic PT2 threshold_generators = threshold_generators_save SOFT_TOUCH threshold_generators else - call ZMQ_selection(to_select, pt2, variance, norm) + call ZMQ_selection(to_select, pt2_data) endif do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & @@ -100,7 +102,11 @@ subroutine run_cipsi correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') - call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -121,7 +127,7 @@ subroutine run_cipsi call diagonalize_CI call save_wavefunction call save_energy(psi_energy_with_nucl_rep, zeros) - if (qp_stop()) exit + if (qp_stop()) exit enddo if (.not.qp_stop()) then @@ -132,13 +138,12 @@ subroutine run_cipsi endif if (do_pt2) then - pt2(:) = 0.d0 - variance(:) = 0.d0 - norm(:) = 0.d0 + pt2_data % pt2(:) = 0.d0 + pt2_data % variance(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 threshold_generators = 1d0 SOFT_TOUCH threshold_generators - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & - norm,0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 SOFT_TOUCH threshold_generators endif print *, 'N_det = ', N_det @@ -148,11 +153,15 @@ subroutine run_cipsi do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo call save_energy(psi_energy_with_nucl_rep, rpt2) - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep(1:N_states), & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 52831c16..ca4b3ba7 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -107,7 +107,7 @@ end function -subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) +subroutine ZMQ_pt2(E, pt2_data, relative_error, error, N_in) use f77_zmq use selection_types @@ -117,10 +117,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) integer, intent(in) :: N_in ! integer, intent(inout) :: N_in double precision, intent(in) :: relative_error, E(N_states) - double precision, intent(out) :: pt2(N_states),error(N_states) - double precision, intent(out) :: variance(N_states),norm2(N_states) - - + double precision, intent(out) :: error(N_states) + type(pt2_type), intent(inout) :: pt2_data +! integer :: i, N double precision :: state_average_weight_save(N_states), w(N_states,4) @@ -138,10 +137,10 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) endif if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then - pt2=0.d0 - variance=0.d0 - norm2=0.d0 - call ZMQ_selection(N_in, pt2, variance, norm2) + pt2_data % pt2=0.d0 + pt2_data % variance=0.d0 + pt2_data % norm2=0.d0 + call ZMQ_selection(N_in, pt2_data % pt2, pt2_data % variance, pt2_data % norm2) error(:) = 0.d0 else @@ -304,10 +303,11 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) if (i==0) then call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, w(1,1), w(1,2), w(1,3), w(1,4), b, N) - pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) + pt2_data % pt2(pt2_stoch_istate) = w(pt2_stoch_istate,1) error(pt2_stoch_istate) = w(pt2_stoch_istate,2) - variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) - norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) + pt2_data % variance(pt2_stoch_istate) = w(pt2_stoch_istate,3) + pt2_data % norm2(pt2_stoch_istate) = w(pt2_stoch_istate,4) + !TODO SEGV else call pt2_slave_inproc(i) @@ -335,10 +335,10 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm2, N_in) TOUCH state_average_weight endif do k=N_det+1,N_states - pt2(k) = 0.d0 + pt2_data % pt2(k) = 0.d0 enddo - call update_pt2_and_variance_weights(pt2, variance, norm2, N_states) + call update_pt2_and_variance_weights(pt2_data, N_states) end subroutine diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 0ca08a0f..288a155c 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -19,15 +19,17 @@ BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] variance_match_weight(:) = 1.d0 END_PROVIDER -subroutine update_pt2_and_variance_weights(pt2, variance, norm2, N_st) +subroutine update_pt2_and_variance_weights(pt2_data, N_st) implicit none + use selection_types BEGIN_DOC ! Updates the PT2- and Variance- matching weights. END_DOC integer, intent(in) :: N_st - double precision, intent(in) :: pt2(N_st) - double precision, intent(in) :: variance(N_st) - double precision, intent(in) :: norm2(N_st) + type(pt2_type), intent(in) :: pt2_data + double precision :: pt2(N_st) + double precision :: variance(N_st) + double precision :: norm2(N_st) double precision :: avg, rpt2(N_st), element, dt, x integer :: k @@ -35,6 +37,10 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm2, N_st) integer, parameter :: i_itermax = 1 double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) + pt2(:) = pt2_data % pt2(:) + variance(:) = pt2_data % variance(:) + norm2(:) = pt2_data % norm2(:) + if (i_iter == 0) then allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) memo_pt2(:,:) = 1.d0 diff --git a/src/cipsi/selection_types.f90 b/src/cipsi/selection_types.f90 index 29e48524..fd76215a 100644 --- a/src/cipsi/selection_types.f90 +++ b/src/cipsi/selection_types.f90 @@ -5,5 +5,11 @@ module selection_types double precision, pointer :: val(:) double precision :: mini endtype + + type pt2_type + double precision, allocatable :: pt2(:) + double precision, allocatable :: variance(:) + double precision, allocatable :: norm2(:) + endtype end module diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index b1534eae..449fd61a 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -1,12 +1,15 @@ subroutine run_stochastic_cipsi + use selection_types implicit none BEGIN_DOC ! Selected Full Configuration Interaction with Stochastic selection and PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm2(:), rpt2(:), zeros(:) + double precision, allocatable :: rpt2(:), zeros(:) integer :: to_select - logical, external :: qp_stop + type(pt2_type) :: pt2_data + logical, external :: qp_stop + double precision :: rss double precision, external :: memory_of_double @@ -19,7 +22,11 @@ subroutine run_stochastic_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm2(N_states), variance(N_states)) + allocate (zeros(N_states), rpt2(N_states)) + allocate( pt2_data % pt2(N_states) ) + allocate( pt2_data % variance(N_states) ) + allocate( pt2_data % norm2(N_states) ) + double precision :: hf_energy_ref logical :: has @@ -28,10 +35,10 @@ subroutine run_stochastic_cipsi relative_error=PT2_relative_error zeros = 0.d0 - pt2 = -huge(1.e0) + pt2_data % pt2 = -huge(1.e0) rpt2 = -huge(1.e0) - norm2 = 0.d0 - variance = huge(1.e0) + pt2_data % norm2 = 0.d0 + pt2_data % variance = huge(1.e0) if (s2_eig) then call make_s2_eigenfunction @@ -65,8 +72,8 @@ subroutine run_stochastic_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(variance(1:N_states))) > variance_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' @@ -75,14 +82,15 @@ subroutine run_stochastic_cipsi to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) - pt2 = 0.d0 - variance = 0.d0 - norm2 = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm2, to_select) ! Stochastic PT2 and selection + + pt2_data % pt2 = 0.d0 + pt2_data % variance = 0.d0 + pt2_data % norm2 = 0.d0 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error, & + to_select) ! Stochastic PT2 and selection do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & @@ -90,7 +98,11 @@ subroutine run_stochastic_cipsi correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') - call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, rpt2) @@ -121,18 +133,21 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, zeros) endif - pt2(:) = 0.d0 - variance(:) = 0.d0 - norm2(:) = 0.d0 - call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & - norm2,0) ! Stochastic PT2 + pt2_data % pt2(:) = 0.d0 + pt2_data % variance(:) = 0.d0 + pt2_data % norm2(:) = 0.d0 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, relative_error, error, 0) ! Stochastic PT2 do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo call save_energy(psi_energy_with_nucl_rep, rpt2) - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 22919117..19e8b33b 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -1,4 +1,4 @@ -subroutine ZMQ_selection(N_in, pt2, variance, norm2) +subroutine ZMQ_selection(N_in, pt2_data) use f77_zmq use selection_types @@ -9,9 +9,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) type(selection_buffer) :: b integer :: i, N integer, external :: omp_get_thread_num - double precision, intent(out) :: pt2(N_states) - double precision, intent(out) :: variance(N_states) - double precision, intent(out) :: norm2(N_states) + type(pt2_type), intent(inout) :: pt2_data ! PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators @@ -112,19 +110,20 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) enddo endif - !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2, variance, norm2) PRIVATE(i) NUM_THREADS(nproc_target+1) + !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2_data) PRIVATE(i) NUM_THREADS(nproc_target+1) i = omp_get_thread_num() if (i==0) then - call selection_collector(zmq_socket_pull, b, N, pt2, variance, norm2) + call selection_collector(zmq_socket_pull, b, N, & + pt2_data % pt2, pt2_data % variance, pt2_data % norm2) else call selection_slave_inproc(i) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'selection') do i=N_det+1,N_states - pt2(i) = 0.d0 - variance(i) = 0.d0 - norm2(i) = 0.d0 + pt2_data % pt2(i) = 0.d0 + pt2_data % variance(i) = 0.d0 + pt2_data % norm2(i) = 0.d0 enddo if (N_in > 0) then if (s2_eig) then @@ -134,12 +133,12 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm2) endif call delete_selection_buffer(b) do k=1,N_states - pt2(k) = pt2(k) * f(k) - variance(k) = variance(k) * f(k) - norm2(k) = norm2(k) * f(k) + pt2_data % pt2(k) = pt2_data % pt2(k) * f(k) + pt2_data % variance(k) = pt2_data % variance(k) * f(k) + pt2_data % norm2(k) = pt2_data % norm2(k) * f(k) enddo - call update_pt2_and_variance_weights(pt2, variance, norm2, N_states) + call update_pt2_and_variance_weights(pt2_data, N_states) end subroutine diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index 83994781..8c076ade 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -28,35 +28,41 @@ end subroutine run implicit none + use selection_types integer :: i,j,k logical, external :: detEq - double precision :: pt2(N_states) + type(pt2_type) :: pt2_data integer :: degree integer :: n_det_before, to_select double precision :: threshold_davidson_in double precision :: E_CI_before(N_states), relative_error, error(N_states), variance(N_states), norm2(N_states), rpt2(N_states) - pt2(:) = 0.d0 + allocate( pt2_data % pt2(N_states) ) + allocate( pt2_data % variance(N_states) ) + allocate( pt2_data % norm2(N_states) ) E_CI_before(:) = psi_energy(:) + nuclear_repulsion relative_error=PT2_relative_error if (do_pt2) then - call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & - norm2,0) ! Stochastic PT2 + call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,relative_error,error,0) ! Stochastic PT2 else - call ZMQ_selection(0, pt2, variance, norm2) + call ZMQ_selection(0, pt2_data) endif do k=1,N_states - rpt2(k) = pt2(k)/(1.d0 + norm2(k)) + rpt2(k) = pt2_data % pt2(k)/(1.d0 + pt2_data % norm2(k)) enddo - call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep(1:N_states), & + pt2_data % pt2, error, & + pt2_data % variance, & + pt2_data % norm2, & + N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(E_CI_before,pt2) + call save_energy(E_CI_before,pt2_data % pt2) end