From 8953f052b477e77b5f598aaa2ac95905e09ad49a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 4 Jan 2021 22:15:46 +0100 Subject: [PATCH 01/22] Accelerate H Psi when c_i=0 --- config/ifort_avx.cfg | 2 +- src/davidson/u0_h_u0.irp.f | 27 +++++++++++++++++++++++++-- 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/config/ifort_avx.cfg b/config/ifort_avx.cfg index f410e3a6..d689050e 100644 --- a/config/ifort_avx.cfg +++ b/config/ifort_avx.cfg @@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags # [OPT] FC : -traceback -FCFLAGS : -mavx -O2 -ip -ftz -g +FCFLAGS : -xAVX -O2 -ip -ftz -g # Profiling flags ################# diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index 302b8423..e56fdec4 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -3,7 +3,7 @@ implicit none BEGIN_DOC ! psi_energy(i) = $\langle \Psi_i | H | \Psi_i \rangle$ -! +! ! psi_s2(i) = $\langle \Psi_i | S^2 | \Psi_i \rangle$ END_DOC call u_0_H_u_0(psi_energy,psi_s2,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size) @@ -250,7 +250,7 @@ compute_singles=.True. !$OMP singles_beta_csc,singles_beta_csc_idx) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & !$OMP lcol, lrow, l_a, l_b, utl, kk, & - !$OMP buffer, doubles, n_doubles, & + !$OMP buffer, doubles, n_doubles, umax, & !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, ratio, & !$OMP n_singles_b, k8, last_found,left,right,right_max) @@ -399,8 +399,11 @@ compute_singles=.True. ! Loop over alpha singles ! ----------------------- + double precision :: umax + !DIR$ LOOP COUNT avg(1000) do k = 1,n_singles_a,block_size + umax = 0.d0 ! Prefetch u_t(:,l_a) do kk=0,block_size-1 if (k+kk > n_singles_a) exit @@ -409,8 +412,10 @@ compute_singles=.True. do l=1,N_st utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) enddo enddo + if (umax < 1.d-20) cycle do kk=0,block_size-1 if (k+kk > n_singles_a) exit @@ -490,6 +495,7 @@ compute_singles=.True. tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) !DIR$ LOOP COUNT avg(1000) do i=1,n_singles_a,block_size + umax = 0.d0 ! Prefetch u_t(:,l_a) do kk=0,block_size-1 if (i+kk > n_singles_a) exit @@ -498,8 +504,10 @@ compute_singles=.True. do l=1,N_st utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) enddo enddo + if (umax < 1.d-20) cycle do kk=0,block_size-1 if (i+kk > n_singles_a) exit @@ -524,6 +532,7 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(50000) do i=1,n_doubles,block_size + umax = 0.d0 ! Prefetch u_t(:,l_a) do kk=0,block_size-1 if (i+kk > n_doubles) exit @@ -532,8 +541,10 @@ compute_singles=.True. do l=1,N_st utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) enddo enddo + if (umax < 1.d-20) cycle do kk=0,block_size-1 if (i+kk > n_doubles) exit @@ -599,6 +610,7 @@ compute_singles=.True. tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) !DIR$ LOOP COUNT avg(1000) do i=1,n_singles_b,block_size + umax = 0.d0 do kk=0,block_size-1 if (i+kk > n_singles_b) exit l_b = singles_b(i+kk) @@ -609,8 +621,10 @@ compute_singles=.True. do l=1,N_st utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) enddo enddo + if (umax < 1.d-20) cycle do kk=0,block_size-1 if (i+kk > n_singles_b) exit @@ -634,6 +648,7 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(50000) do i=1,n_doubles,block_size + umax = 0.d0 do kk=0,block_size-1 if (i+kk > n_doubles) exit l_b = doubles(i+kk) @@ -643,8 +658,10 @@ compute_singles=.True. do l=1,N_st utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) enddo enddo + if (umax < 1.d-20) cycle do kk=0,block_size-1 if (i+kk > n_doubles) exit @@ -671,6 +688,12 @@ compute_singles=.True. ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- + umax = 0.d0 + do l=1,N_st + umax = max(umax, dabs(u_t(l,k_a))) + enddo + if (umax < 1.d-20) cycle + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) From bac477cf39c485de513ef6bd7121314ffe3de24e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 14 Jan 2021 23:28:37 +0100 Subject: [PATCH 02/22] Added debug_cfg --- src/bitmask/bitmasks_routines.irp.f | 54 +++++++++++++++++++++++ src/determinants/create_excitations.irp.f | 3 +- 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 5c4bf347..c34d54dc 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -125,6 +125,41 @@ subroutine bitstring_to_str( output, string, Nint ) output(ibuf:ibuf) = '|' end +subroutine configuration_to_str( output, string, Nint ) + use bitmasks + implicit none + BEGIN_DOC +! Transform the bit string of a configuration to a string for printing + END_DOC + character*(*), intent(out) :: output + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + + integer :: i, j, ibuf + integer(bit_kind) :: itemp + + ibuf = 1 + output = '' + output(ibuf:ibuf) = '|' + ibuf = ibuf+1 + do i=1,Nint + itemp = 1_bit_kind + do j=1,bit_kind_size + if (iand(itemp,string(i,2)) == itemp) then + output(ibuf:ibuf) = '2' + else if (iand(itemp,string(i,1)) == itemp) then + output(ibuf:ibuf) = '1' + else + output(ibuf:ibuf) = '0' + endif + ibuf = ibuf+1 + itemp = shiftl(itemp,1) + enddo + enddo + output(ibuf:ibuf) = '|' +end + + subroutine bitstring_to_hexa( output, string, Nint ) use bitmasks @@ -166,6 +201,25 @@ subroutine debug_det(string,Nint) end +subroutine debug_cfg(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant in '+-' notation and + ! hexadecimal representation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + character*(2048) :: output(2) + call bitstring_to_hexa( output(1), string(1,1), Nint ) + call bitstring_to_hexa( output(2), string(1,2), Nint ) + print *, trim(output(1)) , '|', trim(output(2)) + + call configuration_to_str( output(1), string, Nint ) + print *, trim(output(1)) + +end + subroutine print_det(string,Nint) use bitmasks implicit none diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index cec87901..6f3ec521 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -99,7 +99,8 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin) other_spin(1) = 2 other_spin(2) = 1 if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then - ! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip" + ! There is a spin "ispin" in the orbital i_flip AND + ! There is no electron of opposit spin in the same orbital "i_flip" is_spin_flip_possible = .True. return else From 23f38509042e9daec27064d19396b40272f100ff Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 14 Jan 2021 23:29:56 +0100 Subject: [PATCH 03/22] Added do_single_excitation_cfg, single-exc wrt cfg --- src/determinants/configurations.irp.f | 1 + src/determinants/create_excitations.irp.f | 76 ++++++++++++++++++++++- 2 files changed, 75 insertions(+), 2 deletions(-) diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 156e747c..192376b9 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -584,3 +584,4 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d i += sze enddo END_PROVIDER + diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index 6f3ec521..4500a873 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -24,7 +24,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok) ! check whether position j is occupied if (iand(key_in(k,ispin),mask) /= 0_bit_kind) then key_in(k,ispin) = ibclr(key_in(k,ispin),j) - else + else i_ok= -1 return end if @@ -35,7 +35,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok) mask = ibset(0_bit_kind,j) if (iand(key_in(k,ispin),mask) == 0_bit_kind) then key_in(k,ispin) = ibset(key_in(k,ispin),j) - else + else i_ok= -1 return end if @@ -108,3 +108,75 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin) endif end +subroutine do_single_excitation_cfg(key_in,i_hole,i_particle,i_ok) + use bitmasks + implicit none + BEGIN_DOC + ! Applies the signle excitation operator to a configuration + ! If the excitation is possible, i_ok is True + END_DOC + integer, intent(in) :: i_hole,i_particle + integer(bit_kind), intent(inout) :: key_in(N_int,2) + integer, intent(out) :: i_ok + integer :: k,j,i + integer(bit_kind) :: mask + + ASSERT (i_hole > 0) + ASSERT (i_particle <= mo_num) + + i_ok = 1 + + ! hole + k = shiftr(i_hole-1,bit_kind_shift)+1 + j = i_hole-shiftl(k-1,bit_kind_shift)-1 + mask = ibset(0_bit_kind,j) + + ! Check if the position j is singly occupied + ! 1 -> 0 (SOMO) + ! 0 0 (DOMO) + if (iand(key_in(k,1),mask) /= 0_bit_kind) then + key_in(k,1) = ibclr(key_in(k,1),j) + + ! Check if the position j is doubly occupied + ! 0 -> 1 (SOMO) + ! 1 0 (DOMO) + else if (iand(key_in(k,2),mask) /= 0_bit_kind) then + key_in(k,1) = ibset(key_in(k,1),j) + key_in(k,2) = ibclr(key_in(k,2),j) + + ! The position j is unoccupied: Not OK + ! 0 -> 0 (SOMO) + ! 0 0 (DOMO) + else + i_ok = -1 + return + endif + + + ! particle + k = shiftr(i_particle-1,bit_kind_shift)+1 + j = i_particle-shiftl(k-1,bit_kind_shift)-1 + mask = ibset(0_bit_kind,j) + + ! Check if the position j is singly occupied + ! 1 -> 0 (SOMO) + ! 0 1 (DOMO) + if (iand(key_in(k,1),mask) /= 0_bit_kind) then + key_in(k,1) = ibclr(key_in(k,1),j) + key_in(k,2) = ibset(key_in(k,2),j) + + ! Check if the position j is doubly occupied : Not OK + ! 0 -> 1 (SOMO) + ! 1 0 (DOMO) + else if (iand(key_in(k,2),mask) /= 0_bit_kind) then + i_ok = -1 + return + + ! Position at j is unoccupied + ! 0 -> 0 (SOMO) + ! 0 0 (DOMO) + else + key_in(k,1) = ibset(key_in(k,1),j) + endif + +end From 6a10d02c193e5944374c89577e6577979a31ae62 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 15 Jan 2021 00:07:59 +0100 Subject: [PATCH 04/22] Singly excited cfg --- src/determinants/create_excitations.irp.f | 71 ++++++++++++++++++----- 1 file changed, 55 insertions(+), 16 deletions(-) diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index 4500a873..17173106 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -108,23 +108,25 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin) endif end -subroutine do_single_excitation_cfg(key_in,i_hole,i_particle,i_ok) +subroutine do_single_excitation_cfg(key_in,key_out,i_hole,i_particle,ok) use bitmasks implicit none BEGIN_DOC ! Applies the signle excitation operator to a configuration - ! If the excitation is possible, i_ok is True + ! If the excitation is possible, ok is True END_DOC integer, intent(in) :: i_hole,i_particle integer(bit_kind), intent(inout) :: key_in(N_int,2) - integer, intent(out) :: i_ok + logical , intent(out) :: ok integer :: k,j,i integer(bit_kind) :: mask + integer(bit_kind) :: key_out(N_int,2) ASSERT (i_hole > 0) ASSERT (i_particle <= mo_num) - i_ok = 1 + ok = .True. + key_out(:,:) = key_in(:,:) ! hole k = shiftr(i_hole-1,bit_kind_shift)+1 @@ -134,21 +136,21 @@ subroutine do_single_excitation_cfg(key_in,i_hole,i_particle,i_ok) ! Check if the position j is singly occupied ! 1 -> 0 (SOMO) ! 0 0 (DOMO) - if (iand(key_in(k,1),mask) /= 0_bit_kind) then - key_in(k,1) = ibclr(key_in(k,1),j) + if (iand(key_out(k,1),mask) /= 0_bit_kind) then + key_out(k,1) = ibclr(key_out(k,1),j) ! Check if the position j is doubly occupied ! 0 -> 1 (SOMO) ! 1 0 (DOMO) - else if (iand(key_in(k,2),mask) /= 0_bit_kind) then - key_in(k,1) = ibset(key_in(k,1),j) - key_in(k,2) = ibclr(key_in(k,2),j) + else if (iand(key_out(k,2),mask) /= 0_bit_kind) then + key_out(k,1) = ibset(key_out(k,1),j) + key_out(k,2) = ibclr(key_out(k,2),j) ! The position j is unoccupied: Not OK ! 0 -> 0 (SOMO) ! 0 0 (DOMO) else - i_ok = -1 + ok =.False. return endif @@ -161,22 +163,59 @@ subroutine do_single_excitation_cfg(key_in,i_hole,i_particle,i_ok) ! Check if the position j is singly occupied ! 1 -> 0 (SOMO) ! 0 1 (DOMO) - if (iand(key_in(k,1),mask) /= 0_bit_kind) then - key_in(k,1) = ibclr(key_in(k,1),j) - key_in(k,2) = ibset(key_in(k,2),j) + if (iand(key_out(k,1),mask) /= 0_bit_kind) then + key_out(k,1) = ibclr(key_out(k,1),j) + key_out(k,2) = ibset(key_out(k,2),j) ! Check if the position j is doubly occupied : Not OK ! 0 -> 1 (SOMO) ! 1 0 (DOMO) - else if (iand(key_in(k,2),mask) /= 0_bit_kind) then - i_ok = -1 + else if (iand(key_out(k,2),mask) /= 0_bit_kind) then + ok = .False. return ! Position at j is unoccupied ! 0 -> 0 (SOMO) ! 0 0 (DOMO) else - key_in(k,1) = ibset(key_in(k,1),j) + key_out(k,1) = ibset(key_out(k,1),j) endif end + +subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint) + implicit none + use bitmasks + BEGIN_DOC + ! Generate all single excitation wrt a configuration + ! + ! n_singles : on input, max number of singles : + ! elec_alpha_num * (mo_num - elec_beta_num) + ! on output, number of generated singles + END_DOC + integer, intent(in) :: Nint + integer, intent(inout) :: n_singles + integer(bit_kind), intent(in) :: cfg(Nint,2) + integer(bit_kind), intent(out) :: singles(Nint,2,*) + + integer :: i,k, n_singles_ma, i_hole, i_particle + integer(bit_kind) :: single(Nint,2) + logical :: i_ok + + n_singles = 0 + !TODO + !Make list of Somo and Domo for holes + !Make list of Unocc and Somo for particles + do i_hole = 1, mo_num + do i_particle = 1, mo_num + call do_single_excitation_cfg(cfg,single,i_hole,i_particle,i_ok) + if (i_ok) then + n_singles = n_singles + 1 + do k=1,Nint + singles(k,1,n_singles) = single(k,1) + singles(k,2,n_singles) = single(k,2) + enddo + endif + enddo + enddo +end From 8059bf27458185f9d150b65c5ad8ed57e500cdef Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Jan 2021 18:28:00 +0100 Subject: [PATCH 05/22] Added save_wf_after_selection keyword --- src/cipsi/EZFIO.cfg | 6 ++++++ src/cipsi/cipsi.irp.f | 5 ++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 7d3a5fde..3fdf74d2 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -4,6 +4,12 @@ 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 +interface: ezfio,provider,ocaml +default: False + [seniority_max] type: integer doc: Maximum number of allowed open shells. Using -1 selects all determinants diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 112ebb56..6e715531 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -114,7 +114,10 @@ subroutine run_cipsi ! Add selected determinants call copy_H_apply_buffer_to_wf() -! call save_wavefunction + + if (save_wf_after_selection) then + call save_wavefunction + endif PROVIDE psi_coef PROVIDE psi_det From ec9b50373ab72676556f00de6fb82e010891cec8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Jan 2021 18:29:27 +0100 Subject: [PATCH 06/22] RELEASE_NOTES --- RELEASE_NOTES.org | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index 6e46d868..a913149e 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -53,6 +53,7 @@ - Added ~print_energy~ - Added ~print_hamiltonian~ - Added input for two body RDM + - Added keyword ~save_wf_after_selection~ *** Code @@ -79,3 +80,4 @@ + From de43df3173f8ab6f76ee8bdcf171f57948dc7e25 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Jan 2021 12:35:15 +0100 Subject: [PATCH 07/22] Save wf after selection --- src/cipsi/stochastic_cipsi.irp.f | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index bd0c6a80..77c6ed48 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -104,7 +104,9 @@ subroutine run_stochastic_cipsi ! Add selected determinants call copy_H_apply_buffer_to_wf() -! call save_wavefunction + if (save_wf_after_selection) then + call save_wavefunction + endif PROVIDE psi_coef PROVIDE psi_det From 362dd42ddcba23dea1830275baf41610f4905814 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Jan 2021 13:17:35 +0100 Subject: [PATCH 08/22] Increased task size in Davidson --- src/davidson/davidson_parallel.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 59318cd6..3fdddad4 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -428,7 +428,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) integer :: istep, imin, imax, ishift, ipos integer, external :: add_task_to_taskserver - integer, parameter :: tasksize=10000 + integer, parameter :: tasksize=100000 character*(100000) :: task istep=1 ishift=0 From 340f7076c10ee236b311c66d1363944605052feb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jan 2021 10:00:53 +0100 Subject: [PATCH 09/22] Removed += in selection --- src/cipsi/selection.irp.f | 67 +++++++-------- src/davidson/u0_h_u0.irp.f | 168 ++++++++++++++++++++++++++----------- 2 files changed, 149 insertions(+), 86 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 0edd3f63..8b21798f 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -737,12 +737,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi = mat(istate, p1, p2) - do jstate=1,N_states - pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) - enddo - - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - pt2_data % pt2(istate) += e_pert(istate) + pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate) + pt2_data % variance(istate) = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi + pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate) !!!DEBUG ! delta_E = E0(istate) - Hii + E_shift @@ -1573,7 +1570,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * phase end if - mat(:, p1, p2) += coefs(:) * hij + mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij end do end do else ! AA BB @@ -1590,7 +1587,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, else hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do end if @@ -1649,18 +1646,18 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) end do do putj=hfix+1, mo_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(1:N_states,putj) += hij * coefs(1:N_states) + tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) end do if(ma == 1) then - mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num) + mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num) else - mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num) + mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num) end if end if @@ -1674,22 +1671,22 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, putj = p1 if(.not. banned(putj,puti,bant)) then hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) - tmp_row(:,puti) += hij * coefs(:) + tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) end if putj = p2 if(.not. banned(putj,puti,bant)) then hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) - tmp_row2(:,puti) += hij * coefs(:) + tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) end if end do if(mi == 1) then - mat(:,:,p1) += tmp_row(:,:) - mat(:,:,p2) += tmp_row2(:,:) + mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:) + mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:) else - mat(:,p1,:) += tmp_row(:,:) - mat(:,p2,:) += tmp_row2(:,:) + mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:) + mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:) end if else if(p(0,ma) == 3) then @@ -1702,16 +1699,16 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(:,putj) += hij * coefs(:) + tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) end do do putj=hfix+1,mo_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(:,putj) += hij * coefs(:) + tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) end do - mat(:, :puti-1, puti) += tmp_row(:,:puti-1) - mat(:, puti, puti:) += tmp_row(:,puti:) + mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1) + mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:) end do else hfix = h(1,mi) @@ -1725,19 +1722,19 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, putj = p2 if(.not. banned(puti,putj,1)) then hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) - tmp_row(:,puti) += hij * coefs(:) + tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) end if putj = p1 if(.not. banned(puti,putj,1)) then hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) - tmp_row2(:,puti) += hij * coefs(:) + tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) end if end do - mat(:,:p2-1,p2) += tmp_row(:,:p2-1) - mat(:,p2,p2:) += tmp_row(:,p2:) - mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) - mat(:,p1,p1:) += tmp_row2(:,p1:) + mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1) + mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:) + mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1) + mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:) end if end if deallocate(lbanned) @@ -1760,7 +1757,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) call i_h_j(gen, det, N_int, hij) - mat(:, p1, p2) += coefs(:) * hij + mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij end do end do end @@ -1813,9 +1810,9 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) if(ma == 1) then - mat(:, putj, puti) += coefs(:) * hij + mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij else - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end if end do else @@ -1831,7 +1828,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p1 = p(turn2(i), 1) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do end if @@ -1851,7 +1848,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p1 = p(i1, ma) p2 = p(i2, ma) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do else if(tip == 3) then @@ -1865,7 +1862,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p2 = p(i, ma) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) - mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij + mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij end do else ! tip == 4 puti = p(1, sp) @@ -1876,7 +1873,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, h1 = h(1, mi) h2 = h(2, mi) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int) - mat(:, puti, putj) += coefs(:) * hij + mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end if end if end if diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index e56fdec4..99bc3816 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -211,6 +211,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, double precision :: rss, mem, ratio double precision, allocatable :: utl(:,:) integer, parameter :: block_size=128 + logical :: u_is_sparse ! call resident_memory(rss) ! mem = dble(singles_beta_csc_size) / 1024.d0**3 @@ -222,6 +223,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, ! endif compute_singles=.True. + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 allocate(idx0(maxab)) @@ -249,7 +251,7 @@ compute_singles=.True. !$OMP singles_alpha_csc,singles_alpha_csc_idx, & !$OMP singles_beta_csc,singles_beta_csc_idx) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & - !$OMP lcol, lrow, l_a, l_b, utl, kk, & + !$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, & !$OMP buffer, doubles, n_doubles, umax, & !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, ratio, & @@ -266,6 +268,19 @@ compute_singles=.True. kcol_prev=-1 + ! Check if u has multiple zeros + kk=0 + do k=1,N_det + umax = 0.d0 + do l=1,N_st + umax = max(umax, dabs(u_t(l,k))) + enddo + if (umax < 1.d-20) then + kk = kk+1 + endif + enddo + u_is_sparse = N_det / kk < 10 + ASSERT (iend <= N_det) ASSERT (istart > 0) ASSERT (istep > 0) @@ -405,16 +420,25 @@ compute_singles=.True. do k = 1,n_singles_a,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (k+kk > n_singles_a) exit - l_a = singles_a(k+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (k+kk > n_singles_a) exit + l_a = singles_a(k+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (k+kk > n_singles_a) exit + l_a = singles_a(k+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -497,16 +521,25 @@ compute_singles=.True. do i=1,n_singles_a,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (i+kk > n_singles_a) exit - l_a = singles_a(i+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_singles_a) exit + l_a = singles_a(i+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_singles_a) exit + l_a = singles_a(i+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -534,16 +567,25 @@ compute_singles=.True. do i=1,n_doubles,block_size umax = 0.d0 ! Prefetch u_t(:,l_a) - do kk=0,block_size-1 - if (i+kk > n_doubles) exit - l_a = doubles(i+kk) - ASSERT (l_a <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_a = doubles(i+kk) + ASSERT (l_a <= N_det) - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_a = doubles(i+kk) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -611,19 +653,29 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(1000) do i=1,n_singles_b,block_size umax = 0.d0 - do kk=0,block_size-1 - if (i+kk > n_singles_b) exit - l_b = singles_b(i+kk) - ASSERT (l_b <= N_det) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_singles_b) exit + l_b = singles_b(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) - l_a = psi_bilinear_matrix_transp_order(l_b) - ASSERT (l_a <= N_det) - - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_singles_b) exit + l_b = singles_b(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -649,18 +701,28 @@ compute_singles=.True. !DIR$ LOOP COUNT avg(50000) do i=1,n_doubles,block_size umax = 0.d0 - do kk=0,block_size-1 - if (i+kk > n_doubles) exit - l_b = doubles(i+kk) - ASSERT (l_b <= N_det) - l_a = psi_bilinear_matrix_transp_order(l_b) - ASSERT (l_a <= N_det) - - do l=1,N_st - utl(l,kk+1) = u_t(l,l_a) - umax = max(umax, dabs(utl(l,kk+1))) + if (u_is_sparse) then + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_b = doubles(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + do l=1,N_st + utl(l,kk+1) = u_t(l,l_a) + umax = max(umax, dabs(utl(l,kk+1))) + enddo enddo - enddo + else + do kk=0,block_size-1 + if (i+kk > n_doubles) exit + l_b = doubles(i+kk) + l_a = psi_bilinear_matrix_transp_order(l_b) + ASSERT (l_b <= N_det) + ASSERT (l_a <= N_det) + enddo + umax = 1.d0 + endif if (umax < 1.d-20) cycle do kk=0,block_size-1 @@ -688,10 +750,14 @@ compute_singles=.True. ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - umax = 0.d0 - do l=1,N_st - umax = max(umax, dabs(u_t(l,k_a))) - enddo + if (u_is_sparse) then + umax = 0.d0 + do l=1,N_st + umax = max(umax, dabs(u_t(l,k_a))) + enddo + else + umax = 1.d0 + endif if (umax < 1.d-20) cycle krow = psi_bilinear_matrix_rows(k_a) From 0c3c5fd92649fd6075e4fe6f14fb89eb22ba3cc1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jan 2021 10:52:33 +0100 Subject: [PATCH 10/22] Davidson optimization --- src/davidson/u0_h_u0.irp.f | 7 ++++++- src/mo_two_e_ints/map_integrals.irp.f | 26 ++++++++++++++++++++------ 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index 99bc3816..f4c5638d 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -269,7 +269,7 @@ compute_singles=.True. kcol_prev=-1 ! Check if u has multiple zeros - kk=0 + kk=1 ! Avoid division by zero do k=1,N_det umax = 0.d0 do l=1,N_st @@ -436,6 +436,7 @@ compute_singles=.True. if (k+kk > n_singles_a) exit l_a = singles_a(k+kk) ASSERT (l_a <= N_det) + utl(:,kk+1) = u_t(:,l_a) enddo umax = 1.d0 endif @@ -537,6 +538,7 @@ compute_singles=.True. if (i+kk > n_singles_a) exit l_a = singles_a(i+kk) ASSERT (l_a <= N_det) + utl(:,kk+1) = u_t(:,l_a) enddo umax = 1.d0 endif @@ -583,6 +585,7 @@ compute_singles=.True. if (i+kk > n_doubles) exit l_a = doubles(i+kk) ASSERT (l_a <= N_det) + utl(:,kk+1) = u_t(:,l_a) enddo umax = 1.d0 endif @@ -673,6 +676,7 @@ compute_singles=.True. l_a = psi_bilinear_matrix_transp_order(l_b) ASSERT (l_b <= N_det) ASSERT (l_a <= N_det) + utl(:,kk+1) = u_t(:,l_a) enddo umax = 1.d0 endif @@ -720,6 +724,7 @@ compute_singles=.True. l_a = psi_bilinear_matrix_transp_order(l_b) ASSERT (l_b <= N_det) ASSERT (l_a <= N_det) + utl(:,kk+1) = u_t(:,l_a) enddo umax = 1.d0 endif diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 354dd994..3d301725 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -99,9 +99,15 @@ double precision function get_two_e_integral(i,j,k,l,map) type(map_type), intent(inout) :: map real(integral_kind) :: tmp PROVIDE mo_two_e_integrals_in_map mo_integrals_cache - if (banned_excitation(i,k) .or. banned_excitation(j,l)) then - get_two_e_integral = 0.d0 - return + if (use_banned_excitation) then + if (banned_excitation(i,k)) then + get_two_e_integral = 0.d0 + return + endif + if (banned_excitation(j,l)) then + get_two_e_integral = 0.d0 + return + endif endif ii = l-mo_integrals_cache_min ii = ior(ii, k-mo_integrals_cache_min) @@ -282,17 +288,19 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) end -BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] + BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ logical, use_banned_excitation ] implicit none use map_module BEGIN_DOC ! If true, the excitation is banned in the selection. Useful with local MOs. END_DOC banned_excitation = .False. - integer :: i,j + integer :: i,j, icount integer(key_kind) :: idx double precision :: tmp -! double precision :: buffer(mo_num) + + 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) @@ -300,8 +308,14 @@ BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] 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+1 enddo enddo + use_banned_excitation = (mo_num*mo_num) / icount <= 10 + if (use_banned_excitation) then + print *, 'Using sparsity of exchange integrals' + endif + END_PROVIDER From 7e4d08f51ffa052b79220d5f1b0b716a8b54fd20 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jan 2021 11:03:19 +0100 Subject: [PATCH 11/22] Update Davidson parameters --- src/davidson/davidson_parallel.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 3fdddad4..9b2dc045 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -441,7 +441,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) if (imin==1) then istep = 2 else - istep = 1 + istep = 2 endif do ishift=0,istep-1 write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|' From 9ea755c757307223562b33f6b683ca53cd867f71 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jan 2021 22:28:13 +0100 Subject: [PATCH 12/22] Binary search in det to configuration --- src/determinants/configurations.irp.f | 39 +++++++++------------------ 1 file changed, 13 insertions(+), 26 deletions(-) diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 7e4e6de8..ecd5bfea 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -319,7 +319,7 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] ! Returns the index of the configuration for each determinant END_DOC integer :: i,j,k,r,l - integer*8 :: key + integer*8 :: key, key2 integer(bit_kind) :: occ(N_int,2) logical :: found integer*8, allocatable :: bit_tmp(:) @@ -340,36 +340,23 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] key = configuration_search_key(occ,N_int) - ! TODO: Binary search - l = 1 - r = N_configuration -! do while(r-l > 32) -! j = shiftr(r+l,1) -! if (bit_tmp(j) < key) then -! l = j -! else -! r = j -! endif -! enddo - do j=l,r - found = .True. - do k=1,N_int - if ( (occ(k,1) /= psi_configuration(k,1,j)) & - .or. (occ(k,2) /= psi_configuration(k,2,j)) ) then - found = .False. - exit - endif - enddo - if (found) then + l = 0 + r = N_configuration+1 + j = shiftr(r-l,1) + do while (j>=1) + j = j+l + key2 = configuration_search_key(psi_configuration(1,1,j),N_int) + if (key2 == key) then det_to_configuration(i) = j exit + else if (key2 > key) then + r = j + else + l = j endif + j = shiftr(r-l,1) enddo - if (.not.found) then - print *, '3 bug in ', irp_here - stop -1 - endif enddo !$OMP END PARALLEL DO deallocate(bit_tmp) From a441497098774b2df7838cd3718482cf664053b8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Jan 2021 23:28:52 +0100 Subject: [PATCH 13/22] d-16 threshold --- src/cipsi/selection.irp.f | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 8b21798f..5f331fef 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -166,6 +166,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer :: l_a, nmax, idx integer, allocatable :: indices(:), exc_degree(:), iorder(:) + double precision, parameter :: norm_thr = 1.d-16 allocate (indices(N_det), & exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) @@ -183,7 +184,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d i = psi_bilinear_matrix_rows(l_a) if (nt + exc_degree(i) <= 4) then idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) - if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then + if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 endif @@ -207,7 +208,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d idx = psi_det_sorted_order( & psi_bilinear_matrix_order( & psi_bilinear_matrix_transp_order(l_a))) - if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then + if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 endif From 6553b53d1dbefa51ba027e248c735dd6ef60d98e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 24 Jan 2021 11:22:03 +0100 Subject: [PATCH 14/22] Restore Davidson task size --- src/davidson/davidson_parallel.irp.f | 4 ++-- src/determinants/configurations.irp.f | 1 + src/two_body_rdm/example.irp.f | 4 ++++ 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 9b2dc045..59318cd6 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -428,7 +428,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) integer :: istep, imin, imax, ishift, ipos integer, external :: add_task_to_taskserver - integer, parameter :: tasksize=100000 + integer, parameter :: tasksize=10000 character*(100000) :: task istep=1 ishift=0 @@ -441,7 +441,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) if (imin==1) then istep = 2 else - istep = 2 + istep = 1 endif do ishift=0,istep-1 write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|' diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index ecd5bfea..42e16880 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -42,6 +42,7 @@ subroutine configuration_to_dets_size(o,sze,n_alpha,Nint) amax -= popcnt( o(k,2) ) enddo if (binom_int(bmax, amax) > huge(1)) then + print *, bmax, amax print *, irp_here, ': Too many determinants to generate' stop 1 endif diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 4400613c..c85e560a 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -284,3 +284,7 @@ subroutine routine_full_mos print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3 end + +program test + call routine_active_only +end From 7fa1637b91abd00d730eefa5fb3ad7cd1d3280b7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 24 Jan 2021 23:09:37 +0100 Subject: [PATCH 15/22] Changed banned --- src/davidson/davidson_parallel.irp.f | 2 +- src/davidson/u0_h_u0.irp.f | 5 ++++- src/mo_two_e_ints/map_integrals.irp.f | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 59318cd6..32f89979 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -428,7 +428,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) integer :: istep, imin, imax, ishift, ipos integer, external :: add_task_to_taskserver - integer, parameter :: tasksize=10000 + integer, parameter :: tasksize=20000 character*(100000) :: task istep=1 ishift=0 diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index f4c5638d..bc6acdb9 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -270,16 +270,19 @@ compute_singles=.True. ! Check if u has multiple zeros kk=1 ! Avoid division by zero + !$OMP DO do k=1,N_det umax = 0.d0 do l=1,N_st umax = max(umax, dabs(u_t(l,k))) enddo if (umax < 1.d-20) then + !$OMP ATOMIC kk = kk+1 endif enddo - u_is_sparse = N_det / kk < 10 + !$OMP END DO + u_is_sparse = N_det / kk < 20 ! 5% ASSERT (iend <= N_det) ASSERT (istart > 0) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 3d301725..8756ee47 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -308,10 +308,10 @@ end 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+1 + if (banned_excitation(i,j)) icount = icount+2 enddo enddo - use_banned_excitation = (mo_num*mo_num) / icount <= 10 + use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% if (use_banned_excitation) then print *, 'Using sparsity of exchange integrals' endif From 0ffc173b0b9c3d9add24c6f683c09d62f359f5f2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 10:21:25 +0100 Subject: [PATCH 16/22] Fixed 2rdm compilation --- docs/source/users_guide/quickstart.rst | 6 ++++++ src/two_body_rdm/example.irp.f | 3 --- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/source/users_guide/quickstart.rst b/docs/source/users_guide/quickstart.rst index f0620c5a..e55adcb7 100644 --- a/docs/source/users_guide/quickstart.rst +++ b/docs/source/users_guide/quickstart.rst @@ -128,6 +128,12 @@ and the atomic basis set: ao_two_e_erf_ints density_for_dft electrons mo_two_e_ints scf_utils ao_two_e_ints determinants ezfio nuclei work +If you need to run using an already existing EZFIO database, use + +.. code:: bash + + qp set_file hcn + Run a Hartree-Fock calculation ------------------------------ diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index c85e560a..e9cbd1a2 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -285,6 +285,3 @@ subroutine routine_full_mos end -program test - call routine_active_only -end From 6e1c626f7e4ab5213f692b666acedd3aa2d78339 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 10:24:49 +0100 Subject: [PATCH 17/22] Fixed travis scripts --- travis/compilation.sh | 2 +- travis/configuration.sh | 2 +- travis/testing.sh | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/travis/compilation.sh b/travis/compilation.sh index 1aa26dda..071b4872 100755 --- a/travis/compilation.sh +++ b/travis/compilation.sh @@ -8,7 +8,7 @@ tar -zxf $HOME/cache/config.tgz # Configure QP2 cd qp2 source ./quantum_package.rc -ninja -j 1 -v +ninja -j 1 -v || exit -1 # Create cache cd .. diff --git a/travis/configuration.sh b/travis/configuration.sh index fa52e793..7b3f5423 100755 --- a/travis/configuration.sh +++ b/travis/configuration.sh @@ -2,7 +2,7 @@ # Stage 1 # Configure QP2 -./configure --install all --config ./config/travis.cfg +./configure --install all --config ./config/travis.cfg || exit -1 # Create cache cd ../ diff --git a/travis/testing.sh b/travis/testing.sh index b2122f5c..f67bd106 100755 --- a/travis/testing.sh +++ b/travis/testing.sh @@ -8,7 +8,7 @@ tar -zxf $HOME/cache/compil.tgz # Configure QP2 cd qp2 source ./quantum_package.rc -qp_test -a && rm $HOME/cache/compil.tgz +exec qp_test -a && rm $HOME/cache/compil.tgz From df6d2a73a40c6dc3130df510a89a59c786061201 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 18:41:41 +0100 Subject: [PATCH 18/22] Merge --- src/determinants/configurations.irp.f | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 192376b9..2c141609 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -585,3 +585,21 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d enddo END_PROVIDER + BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] + + implicit none + BEGIN_DOC + ! psi_configuration_to_psi_det_data(k) -> i : i is the index of the determinant in psi_det. + ! + ! psi_configuration_to_psi_det(1:2,k) gives the first and last index of the + ! determinants of configuration k in array psi_configuration_to_psi_det_data. + END_DOC + + integer :: i + + do i=1,N_configuration + + enddo +END_PROVIDER + From 2785a7f4cd8ab2f2b2bccebcbaad0dd8c396762c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 22:32:52 +0100 Subject: [PATCH 19/22] cfg->det --- src/bitmask/bitmasks_routines.irp.f | 54 +++++++++++ src/determinants/configurations.irp.f | 106 ++++++++++++++++---- src/determinants/create_excitations.irp.f | 112 ++++++++++++++++++++++ 3 files changed, 254 insertions(+), 18 deletions(-) diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 5c4bf347..c34d54dc 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -125,6 +125,41 @@ subroutine bitstring_to_str( output, string, Nint ) output(ibuf:ibuf) = '|' end +subroutine configuration_to_str( output, string, Nint ) + use bitmasks + implicit none + BEGIN_DOC +! Transform the bit string of a configuration to a string for printing + END_DOC + character*(*), intent(out) :: output + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + + integer :: i, j, ibuf + integer(bit_kind) :: itemp + + ibuf = 1 + output = '' + output(ibuf:ibuf) = '|' + ibuf = ibuf+1 + do i=1,Nint + itemp = 1_bit_kind + do j=1,bit_kind_size + if (iand(itemp,string(i,2)) == itemp) then + output(ibuf:ibuf) = '2' + else if (iand(itemp,string(i,1)) == itemp) then + output(ibuf:ibuf) = '1' + else + output(ibuf:ibuf) = '0' + endif + ibuf = ibuf+1 + itemp = shiftl(itemp,1) + enddo + enddo + output(ibuf:ibuf) = '|' +end + + subroutine bitstring_to_hexa( output, string, Nint ) use bitmasks @@ -166,6 +201,25 @@ subroutine debug_det(string,Nint) end +subroutine debug_cfg(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant in '+-' notation and + ! hexadecimal representation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + character*(2048) :: output(2) + call bitstring_to_hexa( output(1), string(1,1), Nint ) + call bitstring_to_hexa( output(2), string(1,2), Nint ) + print *, trim(output(1)) , '|', trim(output(2)) + + call configuration_to_str( output(1), string, Nint ) + print *, trim(output(1)) + +end + subroutine print_det(string,Nint) use bitmasks implicit none diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 42e16880..c17ed04e 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -55,7 +55,7 @@ subroutine configuration_to_dets(o,d,sze,n_alpha,Nint) implicit none BEGIN_DOC ! Generate all possible determinants for a given configuration - ! + ! ! Input : ! o : configuration : (doubly occupied, singly occupied) ! sze : Number of produced determinants, computed by `configuration_to_dets_size` @@ -63,7 +63,7 @@ subroutine configuration_to_dets(o,d,sze,n_alpha,Nint) ! Nint : N_int ! ! Output: - ! d : determinants + ! d : determinants ! END_DOC integer ,intent(in) :: Nint @@ -255,16 +255,13 @@ end endif dup = .True. do k=1,N_int - if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & - .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then - dup = .False. - exit - endif + dup = dup .and. (tmp_array(k,1,i) == tmp_array(k,1,j)) & + .and. (tmp_array(k,2,i) == tmp_array(k,2,j)) enddo if (dup) then duplicate(j) = .True. endif - j+=1 + j = j+1 if (j>N_det) then exit endif @@ -287,7 +284,7 @@ end enddo !- Check -! print *, 'Checking for duplicates in occ pattern' +! print *, 'Checking for duplicates in configuration' ! do i=1,N_configuration ! do j=i+1,N_configuration ! duplicate(1) = .True. @@ -314,6 +311,28 @@ end END_PROVIDER +BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ] + implicit none + BEGIN_DOC + ! Returns the index in psi_configuration of the first cfg with + ! the requested seniority + END_DOC + integer :: i, k, s, sold + cfg_seniority_index(:) = -1 + sold = -1 + do i=1,N_configuration + s = 0 + do k=1,N_int + if (psi_configuration(k,1,i) == 0_bit_kind) cycle + s = s + popcnt(psi_configuration(k,1,i)) + enddo + if (s /= sold) then + sold = s + cfg_seniority_index(s) = i + endif + enddo +END_PROVIDER + BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] implicit none BEGIN_DOC @@ -326,7 +345,8 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] integer*8, allocatable :: bit_tmp(:) integer*8, external :: configuration_search_key - allocate(bit_tmp(N_configuration)) + allocate(bit_tmp(0:N_configuration)) + bit_tmp(0) = 0 do i=1,N_configuration bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int) enddo @@ -341,16 +361,30 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] key = configuration_search_key(occ,N_int) + ! Binary search l = 0 r = N_configuration+1 j = shiftr(r-l,1) do while (j>=1) j = j+l - key2 = configuration_search_key(psi_configuration(1,1,j),N_int) - if (key2 == key) then - det_to_configuration(i) = j - exit - else if (key2 > key) then + if (bit_tmp(j) == key) then + do while (bit_tmp(j) == bit_tmp(j-1)) + j = j-1 + enddo + do while (bit_tmp(j) == key) + found = .True. + do k=1,N_int + found = found .and. (psi_configuration(k,1,j) == occ(k,1)) & + .and. (psi_configuration(k,2,j) == occ(k,2)) + enddo + if (found) then + det_to_configuration(i) = j + exit + endif + j = j+1 + enddo + if (found) exit + else if (bit_tmp(j) > key) then r = j else l = j @@ -418,7 +452,7 @@ END_PROVIDER &BEGIN_PROVIDER [ integer, psi_configuration_sorted_order_reverse, (N_configuration) ] implicit none BEGIN_DOC - ! Configurations sorted by weight + ! Configurations sorted by weight END_DOC integer :: i,j,k integer, allocatable :: iorder(:) @@ -430,8 +464,8 @@ END_PROVIDER call dsort(weight_configuration_average_sorted,iorder,N_configuration) do i=1,N_configuration do j=1,N_int - psi_configuration_sorted(j,1,i) = psi_configuration(j,1,iorder(i)) - psi_configuration_sorted(j,2,i) = psi_configuration(j,2,iorder(i)) + psi_configuration_sorted(j,1,i) = psi_configuration(j,1,iorder(i)) + psi_configuration_sorted(j,2,i) = psi_configuration(j,2,iorder(i)) enddo psi_configuration_sorted_order(iorder(i)) = i psi_configuration_sorted_order_reverse(i) = iorder(i) @@ -551,3 +585,39 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d i += sze enddo END_PROVIDER + + BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] + + implicit none + BEGIN_DOC + ! psi_configuration_to_psi_det_data(k) -> i : i is the index of the + ! determinant in psi_det_sorted_bit. ! + ! + ! psi_configuration_to_psi_det(1:2,k) gives the first and last index of the + ! determinants of configuration k in array psi_configuration_to_psi_det_data. + END_DOC + + integer :: i, k, iorder + integer, allocatable :: confs(:) + allocate (confs(N_det)) + + do i=1,N_det + psi_configuration_to_psi_det_data(i) = i + confs(i) = det_to_configuration(i) + enddo + + call isort(confs, psi_configuration_to_psi_det_data, N_det) + k=1 + psi_configuration_to_psi_det(1,1) = 1 + do i=2,N_det + if (confs(i) /= confs(i-1)) then + psi_configuration_to_psi_det(2,k) = i-1 + k = k+1 + psi_configuration_to_psi_det(1,k) = i + endif + enddo + psi_configuration_to_psi_det(2,k) = N_det + +END_PROVIDER + diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index cec87901..3912240c 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -107,3 +107,115 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin) endif end +subroutine do_single_excitation_cfg(key_in,key_out,i_hole,i_particle,ok) + use bitmasks + implicit none + BEGIN_DOC + ! Applies the signle excitation operator to a configuration + ! If the excitation is possible, ok is True + END_DOC + integer, intent(in) :: i_hole,i_particle + integer(bit_kind), intent(in) :: key_in(N_int,2) + logical , intent(out) :: ok + integer :: k,j,i + integer(bit_kind) :: mask + integer(bit_kind) :: key_out(N_int,2) + + ASSERT (i_hole > 0) + ASSERT (i_particle <= mo_num) + + ok = .True. + key_out(:,:) = key_in(:,:) + + ! hole + k = shiftr(i_hole-1,bit_kind_shift)+1 + j = i_hole-shiftl(k-1,bit_kind_shift)-1 + mask = ibset(0_bit_kind,j) + + ! Check if the position j is singly occupied + ! 1 -> 0 (SOMO) + ! 0 0 (DOMO) + if (iand(key_out(k,1),mask) /= 0_bit_kind) then + key_out(k,1) = ibclr(key_out(k,1),j) + + ! Check if the position j is doubly occupied + ! 0 -> 1 (SOMO) + ! 1 0 (DOMO) + else if (iand(key_out(k,2),mask) /= 0_bit_kind) then + key_out(k,1) = ibset(key_out(k,1),j) + key_out(k,2) = ibclr(key_out(k,2),j) + + ! The position j is unoccupied: Not OK + ! 0 -> 0 (SOMO) + ! 0 0 (DOMO) + else + ok =.False. + return + endif + + + ! particle + k = shiftr(i_particle-1,bit_kind_shift)+1 + j = i_particle-shiftl(k-1,bit_kind_shift)-1 + mask = ibset(0_bit_kind,j) + + ! Check if the position j is singly occupied + ! 1 -> 0 (SOMO) + ! 0 1 (DOMO) + if (iand(key_out(k,1),mask) /= 0_bit_kind) then + key_out(k,1) = ibclr(key_out(k,1),j) + key_out(k,2) = ibset(key_out(k,2),j) + + ! Check if the position j is doubly occupied : Not OK + ! 0 -> 1 (SOMO) + ! 1 0 (DOMO) + else if (iand(key_out(k,2),mask) /= 0_bit_kind) then + ok = .False. + return + + ! Position at j is unoccupied + ! 0 -> 0 (SOMO) + ! 0 0 (DOMO) + else + key_out(k,1) = ibset(key_out(k,1),j) + endif + +end + +subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint) + implicit none + use bitmasks + BEGIN_DOC + ! Generate all single excitation wrt a configuration + ! + ! n_singles : on input, max number of singles : + ! elec_alpha_num * (mo_num - elec_beta_num) + ! on output, number of generated singles + END_DOC + integer, intent(in) :: Nint + integer, intent(inout) :: n_singles + integer(bit_kind), intent(in) :: cfg(Nint,2) + integer(bit_kind), intent(out) :: singles(Nint,2,*) + + integer :: i,k, n_singles_ma, i_hole, i_particle + integer(bit_kind) :: single(Nint,2) + logical :: i_ok + + n_singles = 0 + !TODO + !Make list of Somo and Domo for holes + !Make list of Unocc and Somo for particles + do i_hole = 1, mo_num + do i_particle = 1, mo_num + call do_single_excitation_cfg(cfg,single,i_hole,i_particle,i_ok) + if (i_ok) then + n_singles = n_singles + 1 + do k=1,Nint + singles(k,1,n_singles) = single(k,1) + singles(k,2,n_singles) = single(k,2) + enddo + endif + enddo + enddo +end + From 46ce8a3c175074795e5b7a9fbf076a8a995ecd4b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 22:54:59 +0100 Subject: [PATCH 20/22] Fixed doc --- src/determinants/configurations.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index eb0636dd..adea877f 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -592,7 +592,7 @@ END_PROVIDER implicit none BEGIN_DOC ! psi_configuration_to_psi_det_data(k) -> i : i is the index of the - ! determinant in psi_det_sorted_bit + ! determinant in psi_det ! ! psi_configuration_to_psi_det(1:2,k) gives the first and last index of the ! determinants of configuration k in array psi_configuration_to_psi_det_data. From 0f9c1e6cb3523a52cb8f19ea6843f82efa113e65 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Thu, 28 Jan 2021 02:14:59 +0100 Subject: [PATCH 21/22] Working on building the prototype matrix AIJpq #143. --- src/determinants/configurations.org | 85 +++++++++++++++++++ .../configurations_sigma_vector.irp.f | 10 +++ 2 files changed, 95 insertions(+) create mode 100644 src/determinants/configurations.org create mode 100644 src/determinants/configurations_sigma_vector.irp.f diff --git a/src/determinants/configurations.org b/src/determinants/configurations.org new file mode 100644 index 00000000..f6b3e711 --- /dev/null +++ b/src/determinants/configurations.org @@ -0,0 +1,85 @@ +# -*- mode:org -*- +#+TITLE: CFG-CI +#+AUTHOR: Vijay Gopal Chilkuri +#+FILE: configurations.org +#+EMAIL: vijay.gopal.c@gmail.com +#+OPTIONS: toc:t +#+LATEX_CLASS: article +#+LATEX_HEADER: \usepackage{tabularx} +#+LATEX_HEADER: \usepackage{braket} +#+LATEX_HEADER: \usepackage{minted} + +* Configuration based CI + +Here we write the main functions that perform the functions necessary for +the Configuration based CI. + +There are three main functions required for doing the CI + +- Convert wavefunction from determinant basis to configuration state function (CSF) basis + +- Apply the Hamiltonian to the wavefunction in CSF basis + +- Convert the converged wavefunction back to determinant basis + +** TODO[0/3] Convert basis from DET to CSF + +The conversion of basis is done by going via bonded functions (BFs). +Importantly, all the CSFs of a chosen configuration (CFG) are kept. + +The advantage is that the sigma-vector can be performed efficiently +via BLAS level 3 operations. + + +- [ ] Write a function to calculate the maximum dimensions required + + Prototype array contains the \( \) for all possible + CFGs \( I, J\) and all \(4\) types of excitations for all possible model + orbitals \(p,q\). Note that the orbital ids \(p,q\) here do not refer to + the actual MO ids, they simply refer to the orbitals involved in that spefcific + SOMO, for e.g. an excitation of the type [2 2 2 1 1 1 1 0] -> [ 2 2 1 1 1 1 1] + implies an excitation from orbital \(3\) to orbital \(8\) which are the real MO ids. + However, the prototype only concerns the SOMOs like so [2 1 1 1 1 0] -> [ 1 1 1 1 1 1] + therefore \(p,q\) are model space ids \(1,6\). + + #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f + + BEGIN_PROIDER[ integer, NSOMOMax] + implicit none + BEGIN_DOC + ! Documentation for NSOMOMax + ! The maximum number of SOMOs for the current calculation. + ! required for the calculation of prototype arrays. + END_DOC + NSOMOMax = 8 + NCSFMAx = 14 ! TODO: NCSFs for MS=0 + END_PROVIDER + #+end_src + + The prototype matrix AIJpqMatrixList has the following dimensions + ,,\(\left(NSOMOMax, NSOMOMax, 4, NSOMOMax, NSOMOMax,NCSFMAx,NCSFMax\right)\) where the first two + indices represent the somos in \(I,J\) followed by the type of excitation and + finally the two model space orbitals \(p,q\). + +- [ ] Read the transformation matrix based on the number of SOMOs + + #+begin_src f90 :main no + BEGIN_PROIDER[ double precision, AIJpqMatrixList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,NCSFMax,NCSFMax)] + implicit none + BEGIN_DOC + ! Documentation for AIJpqMatrixList + ! The prototype matrix containing the + ! matrices for each I,J somo pair and orb ids. + END_DOC + do i = 1,NSOMOMax + do j = i-1,i+1 + if(j > NSOMOMax || j ==0) then + continue + end if + end do + end do + END_PROVIDER + + #+end_src + +- [ ] Perform the conversion by matrix-vector BLAS level 2 call diff --git a/src/determinants/configurations_sigma_vector.irp.f b/src/determinants/configurations_sigma_vector.irp.f new file mode 100644 index 00000000..a2913af3 --- /dev/null +++ b/src/determinants/configurations_sigma_vector.irp.f @@ -0,0 +1,10 @@ + BEGIN_PROIDER[ integer, NSOMOMax] + implicit none + BEGIN_DOC + ! Documentation for NSOMOMax + ! The maximum number of SOMOs for the current calculation. + ! required for the calculation of prototype arrays. + END_DOC + NSOMOMax = 8 + NCSFMAx = 14 ! TODO: NCSFs for MS=0 + END_PROVIDER From 1a896cf005b01ac66db682f60ca5f9b9bdd45fca Mon Sep 17 00:00:00 2001 From: v1j4y Date: Thu, 28 Jan 2021 20:20:25 +0100 Subject: [PATCH 22/22] Added provider for getting the dimensions of AIJpq tensor #143. --- src/determinants/configurations.org | 53 +++++++++++++++--- .../configurations_sigma_vector.irp.f | 55 ++++++++++++++++++- 2 files changed, 97 insertions(+), 11 deletions(-) diff --git a/src/determinants/configurations.org b/src/determinants/configurations.org index f6b3e711..fb9329ba 100644 --- a/src/determinants/configurations.org +++ b/src/determinants/configurations.org @@ -44,7 +44,8 @@ via BLAS level 3 operations. #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f - BEGIN_PROIDER[ integer, NSOMOMax] + BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NCSFMax] implicit none BEGIN_DOC ! Documentation for NSOMOMax @@ -52,30 +53,64 @@ via BLAS level 3 operations. ! required for the calculation of prototype arrays. END_DOC NSOMOMax = 8 - NCSFMAx = 14 ! TODO: NCSFs for MS=0 + NCSFMax = 14 ! TODO: NCSFs for MS=0 END_PROVIDER #+end_src The prototype matrix AIJpqMatrixList has the following dimensions - ,,\(\left(NSOMOMax, NSOMOMax, 4, NSOMOMax, NSOMOMax,NCSFMAx,NCSFMax\right)\) where the first two + \(\left(NSOMOMax, NSOMOMax, 4, NSOMOMax, NSOMOMax,NCSFMAx,NCSFMax\right)\) where the first two indices represent the somos in \(I,J\) followed by the type of excitation and finally the two model space orbitals \(p,q\). - [ ] Read the transformation matrix based on the number of SOMOs - #+begin_src f90 :main no - BEGIN_PROIDER[ double precision, AIJpqMatrixList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,NCSFMax,NCSFMax)] + #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f + BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)] + use cfunctions implicit none BEGIN_DOC ! Documentation for AIJpqMatrixList ! The prototype matrix containing the ! matrices for each I,J somo pair and orb ids. END_DOC - do i = 1,NSOMOMax - do j = i-1,i+1 - if(j > NSOMOMax || j ==0) then - continue + integer i,j,k,l + integer*8 Isomo, Jsomo + Isomo = 0 + Jsomo = 0 + integer*8 rows, cols + rows = -1 + cols = -1 + integer*8 MS + MS = 0 + print *,"NSOMOMax = ",NSOMOMax + !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) + do i = 2, NSOMOMax, 2 + Isomo = ISHFT(1,i)-1 + do j = i-2,i+2, 2 + Jsomo = ISHFT(1,j)-1 + if(j .GT. NSOMOMax .OR. j .LE. 0) then + cycle end if + do k = 1,NSOMOMax + do l = k,NSOMOMax + call getApqIJMatrixDims(Isomo, & + Jsomo, & + MS, & + rows, & + cols) + print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols + ! i -> j + AIJpqMatrixDimsList(i,j,1,k,l,1) = rows + AIJpqMatrixDimsList(i,j,1,k,l,2) = cols + AIJpqMatrixDimsList(i,j,1,l,k,1) = rows + AIJpqMatrixDimsList(i,j,1,l,k,2) = cols + ! j -> i + AIJpqMatrixDimsList(j,i,1,k,l,1) = rows + AIJpqMatrixDimsList(j,i,1,k,l,2) = cols + AIJpqMatrixDimsList(j,i,1,l,k,1) = rows + AIJpqMatrixDimsList(j,i,1,l,k,2) = cols + end do + end do end do end do END_PROVIDER diff --git a/src/determinants/configurations_sigma_vector.irp.f b/src/determinants/configurations_sigma_vector.irp.f index a2913af3..44665016 100644 --- a/src/determinants/configurations_sigma_vector.irp.f +++ b/src/determinants/configurations_sigma_vector.irp.f @@ -1,4 +1,5 @@ - BEGIN_PROIDER[ integer, NSOMOMax] + BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NCSFMax] implicit none BEGIN_DOC ! Documentation for NSOMOMax @@ -6,5 +7,55 @@ ! required for the calculation of prototype arrays. END_DOC NSOMOMax = 8 - NCSFMAx = 14 ! TODO: NCSFs for MS=0 + NCSFMax = 14 ! TODO: NCSFs for MS=0 + END_PROVIDER + + BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)] + use cfunctions + implicit none + BEGIN_DOC + ! Documentation for AIJpqMatrixList + ! The prototype matrix containing the + ! matrices for each I,J somo pair and orb ids. + END_DOC + integer i,j,k,l + integer*8 Isomo, Jsomo + Isomo = 0 + Jsomo = 0 + integer*8 rows, cols + rows = -1 + cols = -1 + integer*8 MS + MS = 0 + print *,"NSOMOMax = ",NSOMOMax + !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) + do i = 2, NSOMOMax, 2 + Isomo = ISHFT(1,i)-1 + do j = i-2,i+2, 2 + Jsomo = ISHFT(1,j)-1 + if(j .GT. NSOMOMax .OR. j .LE. 0) then + cycle + end if + do k = 1,NSOMOMax + do l = k,NSOMOMax + call getApqIJMatrixDims(Isomo, & + Jsomo, & + MS, & + rows, & + cols) + print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols + ! i -> j + AIJpqMatrixDimsList(i,j,1,k,l,1) = rows + AIJpqMatrixDimsList(i,j,1,k,l,2) = cols + AIJpqMatrixDimsList(i,j,1,l,k,1) = rows + AIJpqMatrixDimsList(i,j,1,l,k,2) = cols + ! j -> i + AIJpqMatrixDimsList(j,i,1,k,l,1) = rows + AIJpqMatrixDimsList(j,i,1,k,l,2) = cols + AIJpqMatrixDimsList(j,i,1,l,k,1) = rows + AIJpqMatrixDimsList(j,i,1,l,k,2) = cols + end do + end do + end do + end do END_PROVIDER