From 8a7d6444abe7d837b33863f0e59ce85b6a95e23b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 18 Oct 2021 08:56:25 +0200 Subject: [PATCH 1/7] removed some stuffs to clean the CIPSI module --- src/cipsi/EZFIO.cfg | 6 - src/cipsi/NEED | 1 - src/cipsi/pert_rdm_providers.irp.f | 183 -------------- src/cipsi/pt2_stoch_routines.irp.f | 2 +- src/cipsi/selection.irp.f | 30 +-- src/cipsi/update_2rdm.irp.f | 223 ------------------ src/cipsi/zmq_selection.irp.f | 2 +- src/davidson_dressed/diagonalize_ci.irp.f | 2 +- src/dft_one_e/mu_erf_dft.irp.f | 1 + src/iterations/print_summary.irp.f | 2 +- src/tools/NEED | 1 + .../print_e_components.irp.f | 0 12 files changed, 21 insertions(+), 432 deletions(-) delete mode 100644 src/cipsi/pert_rdm_providers.irp.f delete mode 100644 src/cipsi/update_2rdm.irp.f rename src/{davidson => two_body_rdm}/print_e_components.irp.f (100%) diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 19b45ac1..7fcf19eb 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -1,9 +1,3 @@ -[pert_2rdm] -type: logical -doc: If true, computes the one- and two-body rdms with perturbation theory -interface: ezfio,provider,ocaml -default: False - [save_wf_after_selection] type: logical doc: If true, saves the wave function after the selection, before the diagonalization diff --git a/src/cipsi/NEED b/src/cipsi/NEED index bfbc559a..85d01f79 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -2,5 +2,4 @@ perturbation zmq mpi iterations -two_body_rdm csf diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f deleted file mode 100644 index eca8decc..00000000 --- a/src/cipsi/pert_rdm_providers.irp.f +++ /dev/null @@ -1,183 +0,0 @@ - -use bitmasks -use omp_lib - -BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock] - use f77_zmq - implicit none - call omp_init_lock(pert_2rdm_lock) -END_PROVIDER - -BEGIN_PROVIDER [integer, n_orb_pert_rdm] - implicit none - n_orb_pert_rdm = n_act_orb -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)] - implicit none - list_orb_reverse_pert_rdm = list_act_reverse - -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)] - implicit none - list_orb_pert_rdm = list_act - -END_PROVIDER - -BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] - implicit none - pert_2rdm_provider = 0.d0 - -END_PROVIDER - -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection) - use bitmasks - use selection_types - implicit none - - integer, intent(in) :: n_det_connection - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - integer, intent(in) :: i_generator, sp, h1, h2 - double precision, intent(in) :: mat(N_states, mo_num, mo_num) - logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) - double precision, intent(in) :: fock_diag_tmp(mo_num) - double precision, intent(in) :: E0(N_states) - type(pt2_type), intent(inout) :: pt2_data - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate, jstate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) - double precision, external :: diag_H_mat_elem_fock - double precision :: E_shift - - logical, external :: detEq - double precision, allocatable :: values(:) - integer, allocatable :: keys(:,:) - integer :: nkeys - integer :: sze_buff - sze_buff = 5 * mo_num ** 2 - allocate(keys(4,sze_buff),values(sze_buff)) - nkeys = 0 - if(sp == 3) then - s1 = 1 - s2 = 2 - else - s1 = sp - s2 = sp - end if - call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) - E_shift = 0.d0 - - if (h0_type == 'CFG') then - j = det_to_configuration(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j) - endif - - do p1=1,mo_num - if(bannedOrb(p1, s1)) cycle - ib = 1 - if(sp /= 3) ib = p1+1 - - do p2=ib,mo_num - -! ----- -! /!\ Generating only single excited determinants doesn't work because a -! determinant generated by a single excitation may be doubly excited wrt -! to a determinant of the future. In that case, the determinant will be -! detected as already generated when generating in the future with a -! double excitation. -! -! if (.not.do_singles) then -! if ((h1 == p1) .or. (h2 == p2)) then -! cycle -! endif -! endif -! -! if (.not.do_doubles) then -! if ((h1 /= p1).and.(h2 /= p2)) then -! cycle -! endif -! endif -! ----- - - if(bannedOrb(p2, s2)) cycle - if(banned(p1,p2)) cycle - - - if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle - call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) - - if (do_only_cas) then - integer, external :: number_of_holes, number_of_particles - if (number_of_particles(det)>0) then - cycle - endif - if (number_of_holes(det)>0) then - cycle - endif - endif - - if (do_ddci) then - logical, external :: is_a_two_holes_two_particles - if (is_a_two_holes_two_particles(det)) then - cycle - endif - endif - - if (do_only_1h1p) then - logical, external :: is_a_1h1p - if (.not.is_a_1h1p(det)) cycle - endif - - - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - - sum_e_pert = 0d0 - integer :: degree - call get_excitation_degree(det,HF_bitmask,degree,N_int) - if(degree == 2)cycle - do istate=1,N_states - delta_E = E0(istate) - Hii + E_shift - alpha_h_psi = mat(istate, p1, p2) - val = alpha_h_psi + alpha_h_psi - tmp = dsqrt(delta_E * delta_E + val * val) - if (delta_E < 0.d0) then - tmp = -tmp - endif - e_pert = 0.5d0 * (tmp - delta_E) - coef(istate) = e_pert / alpha_h_psi - print*,e_pert,coef,alpha_h_psi - pt2_data % pt2(istate) += e_pert - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - enddo - - do istate=1,N_states - alpha_h_psi = mat(istate, p1, p2) - e_pert = coef(istate) * alpha_h_psi - do jstate=1,N_states - pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate) - enddo - - if (weight_selection /= 5) then - ! Energy selection - sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) - - else - ! Variance selection - sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) - endif - end do - call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - - if(sum_e_pert <= buf%mini) then - call add_to_selection_buffer(buf, det, sum_e_pert) - end if - end do - end do - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -end - - diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index b366a268..3594aaf2 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then PROVIDE psi_configuration_hii det_to_configuration diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index eda9642c..f1ec6ff6 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -464,14 +464,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) - if(pert_2rdm)then - allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) - do i=1,fullinteresting(0) - do j = 1, N_states - coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) - enddo - enddo - endif +! if(pert_2rdm)then +! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) +! do i=1,fullinteresting(0) +! do j = 1, N_states +! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) +! enddo +! enddo +! endif do i=1,fullinteresting(0) do k=1,N_int @@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - if(.not.pert_2rdm)then +! if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) - else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) - endif +! else +! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) +! endif end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) - if(pert_2rdm)then - deallocate(coef_fullminilist_rev) - endif +! if(pert_2rdm)then +! deallocate(coef_fullminilist_rev) +! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f deleted file mode 100644 index 260c48fd..00000000 --- a/src/cipsi/update_2rdm.irp.f +++ /dev/null @@ -1,223 +0,0 @@ -use bitmasks - -subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: n_det_connection,sze_buff - double precision, intent(in) :: coef(N_states) - integer(bit_kind), intent(in) :: det(N_int,2) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer, intent(inout) :: keys(4,sze_buff),nkeys - double precision, intent(inout) :: values(sze_buff) - integer :: i,j - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase, contrib - do i = 1, n_det_connection - call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) - if(degree.gt.2)cycle - contrib = 0.d0 - do j = 1, N_states - contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j) - enddo - ! case of single excitations - if(degree == 1)then - if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - nkeys = 0 - endif - call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff) - else - !! case of double excitations - ! if (nkeys + 4 .ge. sze_buff)then - ! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - ! nkeys = 0 - ! endif - ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - endif - enddo -!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -!nkeys = 0 - -end - -subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer(bit_kind), intent(in) :: det1(N_int,2) - integer(bit_kind), intent(in) :: det2(N_int,2) - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2),ispin,other_spin - integer :: h1,h2,p1,p2,i - call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int) - - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - p1 = exc(1,2,1) - ispin = 1 - other_spin = 2 - else - ! Mono beta - h1 = exc(1,1,2) - p1 = exc(1,2,2) - ispin = 2 - other_spin = 1 - endif - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - !update the alpha/beta part - do i = 1, n_occ_ab(other_spin) - h2 = occ(i,other_spin) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - !update the same spin part -!do i = 1, n_occ_ab(ispin) -! h2 = occ(i,ispin) -! if(list_orb_reverse_pert_rdm(h2).lt.0)return -! h2 = list_orb_reverse_pert_rdm(h2) - -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 -! -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 -!enddo - -end - -subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - integer :: h1,h2,p1,p2 - - if (exc(0,1,1) == 1) then - ! Double alpha/beta - h1 = exc(1,1,1) - h2 = exc(1,1,2) - p1 = exc(1,2,1) - p2 = exc(1,2,2) - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - - else - if (exc(0,1,1) == 2) then - ! Double alpha/alpha - h1 = exc(1,1,1) - h2 = exc(2,1,1) - p1 = exc(1,2,1) - p2 = exc(2,2,1) - else if (exc(0,1,2) == 2) then - ! Double beta - h1 = exc(1,1,2) - h2 = exc(2,1,2) - p1 = exc(1,2,2) - p2 = exc(2,2,2) - endif - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - -end - - diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 58630709..1bfe87c0 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') diff --git a/src/davidson_dressed/diagonalize_ci.irp.f b/src/davidson_dressed/diagonalize_ci.irp.f index 7619532a..709ee0e6 100644 --- a/src/davidson_dressed/diagonalize_ci.irp.f +++ b/src/davidson_dressed/diagonalize_ci.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] enddo do j=1,min(N_det,N_states) write(st,'(I4)') j - call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) + call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st)) call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) enddo diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index a1d7a5c1..21f6dc21 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -58,6 +58,7 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft_average] r(3) = final_grid_points(3,i) call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) rho = dm_a + dm_b + if(mu_of_r_dft(i).gt.1.d+3)cycle mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i) enddo mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 641ee209..8e6285e2 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s enddo endif - call print_energy_components() +! call print_energy_components() end subroutine diff --git a/src/tools/NEED b/src/tools/NEED index c07c9109..0f4e17b0 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -2,3 +2,4 @@ fci mo_two_e_erf_ints aux_quantities hartree_fock +two_body_rdm diff --git a/src/davidson/print_e_components.irp.f b/src/two_body_rdm/print_e_components.irp.f similarity index 100% rename from src/davidson/print_e_components.irp.f rename to src/two_body_rdm/print_e_components.irp.f From 699f655b8943b5aa7486be845294d8dec8d86305 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 25 Oct 2021 10:35:22 +0200 Subject: [PATCH 2/7] added int grad --- src/dft_utils_in_r/ints_grad.irp.f | 39 ++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/dft_utils_in_r/ints_grad.irp.f diff --git a/src/dft_utils_in_r/ints_grad.irp.f b/src/dft_utils_in_r/ints_grad.irp.f new file mode 100644 index 00000000..239fe554 --- /dev/null +++ b/src/dft_utils_in_r/ints_grad.irp.f @@ -0,0 +1,39 @@ + BEGIN_PROVIDER [ double precision, mo_grad_ints, (mo_num, mo_num,3)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + mo_grad_ints = 0.d0 + do m = 1, 3 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints(i,j,m) += mos_grad_in_r_array(j,ipoint,m) * mos_in_r_array(i,ipoint) * weight + enddo + enddo + enddo + enddo + + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, mo_grad_ints_transp, (3,mo_num, mo_num)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + do m = 1, 3 + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints_transp(m,i,j) = mo_grad_ints(i,j,m) + enddo + enddo + enddo + + +END_PROVIDER From 1f0c48023d20540a3101f5c46b9317479a980eb6 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 26 Oct 2021 19:11:04 +0200 Subject: [PATCH 3/7] dav_double_dress_ext_rout.irp.f --- .../dav_double_dress_ext_rout.irp.f | 520 ++++++++++++++++++ src/dav_general_mat/dav_ext_rout.irp.f | 2 +- src/utils/sort.irp.f | 209 ++----- 3 files changed, 557 insertions(+), 174 deletions(-) create mode 100644 src/dav_general_mat/dav_double_dress_ext_rout.irp.f diff --git a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f new file mode 100644 index 00000000..9bcfc5ab --- /dev/null +++ b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f @@ -0,0 +1,520 @@ +subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc) + use mmap_module + BEGIN_DOC + ! Generic Davidson diagonalization with TWO DRESSING VECTORS + ! + ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian + ! + ! Dressing_vec : COLUMN / LINE DRESSING VECTOR + ! + ! idx_dress : position of the basis function used to use the Dressing_vec (usually the largest coeff) + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! u_in : guess coefficients on the various states. Overwritten on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + implicit none + integer, intent(in) :: sze, N_st, N_st_diag, idx_dress + double precision, intent(in) :: H_jj(sze),Dress_jj(sze),Dressing_vec(sze,N_st) + double precision, intent(inout) :: u_in(sze,N_st_diag) + double precision, intent(out) :: energies(N_st_diag) + logical, intent(out) :: converged + external hcalc + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + do k=1,N_st + do i=1,sze + H_jj_tmp(i) += u_in(i,k) * Dressing_vec(i,k) + enddo + enddo + + integer :: iter + integer :: i,j,k,l,m + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: s_tmp(:,:) + double precision, allocatable :: residual_norm(:),inv_c_idx_dress_vec(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st),inv_c_idx_dress + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + logical :: state_ok(N_st_diag*davidson_sze_max) + integer :: nproc_target + integer :: order(N_st_diag) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag*davidson_sze_max) + + + allocate(inv_c_idx_dress_vec(N_st)) + inv_c_idx_dress = 1.d0/u_in(idx_dress,1) + do i = 1, N_st + inv_c_idx_dress_vec(i) = 1.d0/u_in(idx_dress,i) + enddo + include 'constants.include.F' + + integer :: N_st_diag_in + N_st_diag_in = N_st_diag + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag_in*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag_in*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag_in))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag_in*itermax, N_st_diag_in*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag_in*itermax) &! U + + 1.0d0*dble(sze*m)*(N_st_diag_in*itermax) &! W + + 3.0d0*(N_st_diag_in*itermax)**2 &! h,y,s_tmp + + 1.d0*(N_st_diag_in*itermax) &! lambda + + 1.d0*(N_st_diag_in) &! residual_norm + ! In H_u_0_nstates_zmq + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on collector + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on slave + + 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag_in,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions ') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag_in*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag_in*itermax), & + + ! Small + h(N_st_diag_in*itermax,N_st_diag_in*itermax), & + y(N_st_diag_in*itermax,N_st_diag_in*itermax), & + s_tmp(N_st_diag_in*itermax,N_st_diag_in*itermax), & + residual_norm(N_st_diag_in), & + lambda(N_st_diag_in*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + s_tmp = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag_in >= N_st) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + ! Davidson iterations + ! =================== + + converged = .False. + + do k=N_st+1,N_st_diag_in + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + enddo + u_in(k,k) = u_in(k,k) + 10.d0 + enddo + do k=1,N_st_diag_in + call normalize(u_in(1,k),sze) + enddo + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag_in*(iter-1) + shift2 = N_st_diag_in*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag_in,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + + if (N_st == 1) then + + l = idx_dress + double precision :: f + f = inv_c_idx_dress + do istate=1,N_st_diag_in + do i=1,sze + W(i,shift+istate) += Dressing_vec(i,1) *f * U(l,shift+istate) + W(l,shift+istate) += Dressing_vec(i,1) *f * U(i,shift+istate) + enddo + enddo + + else + print*,'dav_double_dressed routine not yet implemented for N_st > 1' +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! psi_coef, size(psi_coef,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! Dressing_vec, size(Dressing_vec,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! Dressing_vec, size(Dressing_vec,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + ! Diagonalize h + ! --------------- + + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif + + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag_in + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag_in + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + do k=1,N_st + if (residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag_in + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + enddo + + + call nullify_small_elements(sze,N_st_diag_in,U,size(U,1),threshold_davidson_pt2) + do k=1,N_st_diag_in + do i=1,sze + u_in(i,k) = U(i,k) + enddo + enddo + + do k=1,N_st_diag_in + energies(k) = lambda(k) + enddo + write_buffer = '======' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) + + deallocate(W) + + deallocate ( & + residual_norm, & + U, overlap, & + h, y, s_tmp, & + lambda & + ) + FREE nthreads_davidson +end + + +subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze) + implicit none + BEGIN_DOC + ! Routine that computes the diagonal part of the dressing + ! + ! v(i) += u(i) * dress_diag(i) + ! + ! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized + ! + ! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st),dress_diag(sze) + double precision, intent(inout) :: v(sze,N_st) + integer :: i,istate + do istate = 1, N_st + do i = 1, sze + v(i,istate) += dress_diag(i) * u(i,istate) + enddo + enddo +end + + + + + + + + + + + diff --git a/src/dav_general_mat/dav_ext_rout.irp.f b/src/dav_general_mat/dav_ext_rout.irp.f index ff2167e4..0f2d7680 100644 --- a/src/dav_general_mat/dav_ext_rout.irp.f +++ b/src/dav_general_mat/dav_ext_rout.irp.f @@ -3,7 +3,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co use mmap_module implicit none BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix + ! Generic Davidson diagonalization ! ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson ! diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index a63eb4a3..2a655eed 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -38,7 +38,15 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer, external :: omp_get_num_threads - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + if (omp_get_num_threads() == 1) then + !$OMP PARALLEL DEFAULT(SHARED) + !$OMP SINGLE + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + !$OMP END SINGLE + !$OMP END PARALLEL + else + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + endif end recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) @@ -49,7 +57,7 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j - + if(isize<2)return c = x( shiftr(first+last,1) ) @@ -81,11 +89,16 @@ BEGIN_TEMPLATE endif else if (first < i-1) then + !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level) call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) + !$OMP END TASK endif if (j+1 < last) then + !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level) call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) + !$OMP END TASK endif + !$OMP TASKWAIT endif end @@ -249,60 +262,7 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE - -!---------------------- INTEL -IRP_IF INTEL - BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - character, allocatable :: tmp(:) - if (isize < 2) return - call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) - allocate(tmp(n)) - call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) - deallocate(tmp) - iorder(1:isize) = iorder(1:isize)+1 - call $Xset_order(x,iorder,isize) - end - - subroutine $Xsort_noidx(x,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer :: n - character, allocatable :: tmp(:) - if (isize < 2) return - call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) - allocate(tmp(n)) - call ippsSortRadixAscend_$ityp_I(x, isize, tmp) - deallocate(tmp) - end - -SUBST [ X, type, ityp, n, ippsz ] - ; real ; 32f ; 4 ; 13 ;; - i ; integer ; 32s ; 4 ; 11 ;; - i2 ; integer*2 ; 16s ; 2 ; 7 ;; -END_TEMPLATE - -BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -329,12 +289,12 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type ] - d ; double precision ;; +SUBST [ X, type, Y ] + ; real ; i ;; + d ; double precision ; i8 ;; END_TEMPLATE BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -346,112 +306,8 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n - if (isize < 2) then - return - endif - call sorted_$Xnumber(x,isize,n) - if (isize == n) then - return - endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else - call $Xradix_sort(x,iorder,isize,-1) - endif - end subroutine $Xsort - -SUBST [ X, type ] - i8 ; integer*8 ;; -END_TEMPLATE - -!---------------------- END INTEL -IRP_ELSE -!---------------------- NON-INTEL -BEGIN_TEMPLATE - - subroutine $Xsort_noidx(x,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer, allocatable :: iorder(:) - integer :: i - allocate(iorder(isize)) - do i=1,isize - iorder(i)=i - enddo - call $Xsort(x,iorder,isize) - deallocate(iorder) - end subroutine $Xsort_noidx - -SUBST [ X, type ] - ; real ;; - d ; double precision ;; - i ; integer ;; - i8 ; integer*8 ;; - i2 ; integer*2 ;; -END_TEMPLATE - -BEGIN_TEMPLATE - - subroutine $Xsort(x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - if (isize < 2) then - return - endif -! call sorted_$Xnumber(x,isize,n) -! if (isize == n) then -! return -! endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else -! call heap_$Xsort(x,iorder,isize) - call quick_$Xsort(x,iorder,isize) - endif - end subroutine $Xsort - -SUBST [ X, type ] - ; real ;; - d ; double precision ;; -END_TEMPLATE - -BEGIN_TEMPLATE - - subroutine $Xsort(x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - if (isize < 2) then - return - endif - call sorted_$Xnumber(x,isize,n) - if (isize == n) then - return - endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else - call $Xradix_sort(x,iorder,isize,-1) - endif +! call $Xradix_sort(x,iorder,isize,-1) + call quick_$Xsort(x,iorder,isize) end subroutine $Xsort SUBST [ X, type ] @@ -460,11 +316,6 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE -IRP_ENDIF -!---------------------- END NON-INTEL - - - BEGIN_TEMPLATE subroutine $Xset_order(x,iorder,isize) implicit none @@ -562,12 +413,10 @@ SUBST [ X, type ] i2; integer*2 ;; END_TEMPLATE - BEGIN_TEMPLATE -recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) + recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) implicit none - BEGIN_DOC ! Sort integer array x(isize) using the radix sort algorithm. ! iorder in input should be (1,2,3,...,isize), and in output @@ -703,14 +552,24 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) endif +! !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000) +! !$OMP SINGLE if (i3>1_$int_type) then +! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000) call $Xradix_sort$big(x,iorder,i3,iradix_new-1) +! !$OMP END TASK endif if (isize-i3>1_$int_type) then +! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) +! !$OMP END TASK endif +! !$OMP TASKWAIT +! !$OMP END SINGLE +! !$OMP END PARALLEL + return endif @@ -765,11 +624,16 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) if (i1>1_$int_type) then + !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) + !$OMP END TASK endif if (i0>1) then + !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000) call $Xradix_sort$big(x,iorder,i0,iradix-1) + !$OMP END TASK endif + !$OMP TASKWAIT end @@ -782,4 +646,3 @@ SUBST [ X, type, integer_size, is_big, big, int_type ] END_TEMPLATE - From fccd7e2d1a22c3e1599311cd431bed602314cc71 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 26 Oct 2021 20:49:48 +0200 Subject: [PATCH 4/7] reput the sort.irp.f --- src/utils/sort.irp.f | 209 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 173 insertions(+), 36 deletions(-) diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index 2a655eed..a63eb4a3 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -38,15 +38,7 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer, external :: omp_get_num_threads - if (omp_get_num_threads() == 1) then - !$OMP PARALLEL DEFAULT(SHARED) - !$OMP SINGLE - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - !$OMP END SINGLE - !$OMP END PARALLEL - else - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - endif + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) end recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) @@ -57,7 +49,7 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j - + if(isize<2)return c = x( shiftr(first+last,1) ) @@ -89,16 +81,11 @@ BEGIN_TEMPLATE endif else if (first < i-1) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level) call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) - !$OMP END TASK endif if (j+1 < last) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level) call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) - !$OMP END TASK endif - !$OMP TASKWAIT endif end @@ -262,7 +249,60 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE + +!---------------------- INTEL +IRP_IF INTEL + BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) + deallocate(tmp) + iorder(1:isize) = iorder(1:isize)+1 + call $Xset_order(x,iorder,isize) + end + + subroutine $Xsort_noidx(x,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixAscend_$ityp_I(x, isize, tmp) + deallocate(tmp) + end + +SUBST [ X, type, ityp, n, ippsz ] + ; real ; 32f ; 4 ; 13 ;; + i ; integer ; 32s ; 4 ; 11 ;; + i2 ; integer*2 ; 16s ; 2 ; 7 ;; +END_TEMPLATE + +BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -289,12 +329,12 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type, Y ] - ; real ; i ;; - d ; double precision ; i8 ;; +SUBST [ X, type ] + d ; double precision ;; END_TEMPLATE BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -306,8 +346,112 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif + end subroutine $Xsort + +SUBST [ X, type ] + i8 ; integer*8 ;; +END_TEMPLATE + +!---------------------- END INTEL +IRP_ELSE +!---------------------- NON-INTEL +BEGIN_TEMPLATE + + subroutine $Xsort_noidx(x,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call $Xsort(x,iorder,isize) + deallocate(iorder) + end subroutine $Xsort_noidx + +SUBST [ X, type ] + ; real ;; + d ; double precision ;; + i ; integer ;; + i8 ; integer*8 ;; + i2 ; integer*2 ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif +! call sorted_$Xnumber(x,isize,n) +! if (isize == n) then +! return +! endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else +! call heap_$Xsort(x,iorder,isize) + call quick_$Xsort(x,iorder,isize) + endif + end subroutine $Xsort + +SUBST [ X, type ] + ; real ;; + d ; double precision ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif end subroutine $Xsort SUBST [ X, type ] @@ -316,6 +460,11 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE +IRP_ENDIF +!---------------------- END NON-INTEL + + + BEGIN_TEMPLATE subroutine $Xset_order(x,iorder,isize) implicit none @@ -413,10 +562,12 @@ SUBST [ X, type ] i2; integer*2 ;; END_TEMPLATE + BEGIN_TEMPLATE - recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) +recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) implicit none + BEGIN_DOC ! Sort integer array x(isize) using the radix sort algorithm. ! iorder in input should be (1,2,3,...,isize), and in output @@ -552,24 +703,14 @@ BEGIN_TEMPLATE endif -! !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000) -! !$OMP SINGLE if (i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000) call $Xradix_sort$big(x,iorder,i3,iradix_new-1) -! !$OMP END TASK endif if (isize-i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) -! !$OMP END TASK endif -! !$OMP TASKWAIT -! !$OMP END SINGLE -! !$OMP END PARALLEL - return endif @@ -624,16 +765,11 @@ BEGIN_TEMPLATE if (i1>1_$int_type) then - !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) - !$OMP END TASK endif if (i0>1) then - !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000) call $Xradix_sort$big(x,iorder,i0,iradix-1) - !$OMP END TASK endif - !$OMP TASKWAIT end @@ -646,3 +782,4 @@ SUBST [ X, type, integer_size, is_big, big, int_type ] END_TEMPLATE + From b19569905057f36aedbd478620d57cf4f8b4542b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 28 Oct 2021 15:19:45 +0200 Subject: [PATCH 5/7] added H_matrix_diag_all_dets --- src/determinants/utils.irp.f | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 3aec16f9..bd2ee293 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -20,6 +20,25 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] !$OMP END PARALLEL DO END_PROVIDER +BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ] + use bitmasks + implicit none + BEGIN_DOC + ! |H| matrix on the basis of the Slater determinants defined by psi_det + END_DOC + integer :: i + double precision :: hij + integer :: degree(N_det) + call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) & + !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets) + do i =1,N_det + call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij) + H_matrix_diag_all_dets(i) = hij + enddo + !$OMP END PARALLEL DO +END_PROVIDER + BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ] use bitmasks From 18bb5ee9172841bd4c8a33ca144f65be2b469a22 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 29 Oct 2021 11:44:37 +0200 Subject: [PATCH 6/7] added davidson for diagonal dress --- .../dav_diag_dressed_ext_rout.irp.f | 481 ++++++++++++++++++ .../dav_double_dress_ext_rout.irp.f | 2 +- .../dav_dressed_ext_rout.irp.f | 4 +- src/dav_general_mat/dav_ext_rout.irp.f | 8 +- src/determinants/utils.irp.f | 2 + 5 files changed, 491 insertions(+), 6 deletions(-) create mode 100644 src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f diff --git a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f new file mode 100644 index 00000000..2f3d7f80 --- /dev/null +++ b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f @@ -0,0 +1,481 @@ + +subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) + use mmap_module + implicit none + BEGIN_DOC + ! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR + ! + ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! u_in : guess coefficients on the various states. Overwritten on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + integer, intent(in) :: sze, N_st, N_st_diag_in + double precision, intent(in) :: H_jj(sze),Dress_jj(sze) + double precision, intent(inout) :: u_in(sze,N_st_diag_in) + double precision, intent(out) :: energies(N_st) + external hcalc + + integer :: iter, N_st_diag + integer :: i,j,k,l,m + logical, intent(inout) :: converged + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + integer :: nproc_target + integer :: order(N_st_diag_in) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:)!, S_d(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag_in*davidson_sze_max) + + include 'constants.include.F' + + N_st_diag = N_st_diag_in + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag*itermax, N_st_diag*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + provide threshold_davidson !nthreads_davidson + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = sze + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag*itermax) &! U + + 1.d0*dble(sze*m)*(N_st_diag*itermax) &! W + + 2.0d0*(N_st_diag*itermax)**2 &! h,y + + 2.d0*(N_st_diag*itermax) &! s2,lambda + + 1.d0*(N_st_diag) &! residual_norm + ! In H_S2_u_0_nstates_zmq + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on collector + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on slave + + 0.5d0*maxab &! idx0 in H_S2_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag*itermax), & + ! Small + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + residual_norm(N_st_diag), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + ! Initialize from N_st to N_st_diat with gaussian random numbers + ! to be sure to have overlap with any eigenvectors + do k=N_st+1,N_st_diag + u_in(k,k) = 10.d0 + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) + enddo + enddo + ! Normalize all states + do k=1,N_st_diag + call normalize(u_in(1,k),sze) + enddo + + ! Copy from the guess input "u_in" to the working vectors "U" + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + + ! Gram-Schmidt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) + call ortho_qr(U,size(U,1),sze,shift2) + ! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag) + ! where sze is the size of the vector, N_st_diag is the number of states + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + + ! Diagonalize h y = lambda y + ! --------------- + + call lapack_diag(lambda,y,h,size(h,1),shift2) + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + do k=1,N_st + l = order(k) + if (k /= l) then + lambda(k) = overlap(l,1) + endif + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + + do k=1,N_st + if (residual_norm(k) > 1.e8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + do j=1,N_st_diag + k=1 + do while ((k 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - ! Gram-Schmidt to orthogonalize all new guess with the previous vectors call ortho_qr(U,size(U,1),sze,shift2) call ortho_qr(U,size(U,1),sze,shift2) @@ -345,6 +344,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co enddo + ! Re-contract U and update W + ! -------------------------------- + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) do k=1,N_st_diag @@ -360,8 +362,8 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co U(i,k) = u_in(i,k) enddo enddo - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) do j=1,N_st_diag k=1 do while ((k Date: Fri, 5 Nov 2021 15:17:39 +0100 Subject: [PATCH 7/7] removed stupid assert --- .../dav_double_dress_ext_rout.irp.f | 2 -- src/mo_two_e_ints/map_integrals.irp.f | 30 +++++++++---------- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f index 4f4572fa..884fd672 100644 --- a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f +++ b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f @@ -207,8 +207,6 @@ subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies ASSERT (N_st > 0) ASSERT (N_st_diag_in >= N_st) ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) ! Davidson iterations ! =================== diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 9f73d518..272916e3 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -302,21 +302,21 @@ end integer(key_kind) :: idx double precision :: tmp - icount = 1 ! Avoid division by zero - do j=1,mo_num - do i=1,j-1 - call two_e_integrals_index(i,j,j,i,idx) - !DIR$ FORCEINLINE - call map_get(mo_integrals_map,idx,tmp) - banned_excitation(i,j) = dabs(tmp) < 1.d-14 - banned_excitation(j,i) = banned_excitation(i,j) - if (banned_excitation(i,j)) icount = icount+2 - enddo - enddo - use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% - if (use_banned_excitation) then - print *, 'Using sparsity of exchange integrals' - endif +!icount = 1 ! Avoid division by zero +!do j=1,mo_num +! do i=1,j-1 +! call two_e_integrals_index(i,j,j,i,idx) +! !DIR$ FORCEINLINE +! call map_get(mo_integrals_map,idx,tmp) +! banned_excitation(i,j) = dabs(tmp) < 1.d-14 +! banned_excitation(j,i) = banned_excitation(i,j) +! if (banned_excitation(i,j)) icount = icount+2 +! enddo +!enddo +!use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% +!if (use_banned_excitation) then +! print *, 'Using sparsity of exchange integrals' +!endif END_PROVIDER