From 2785a7f4cd8ab2f2b2bccebcbaad0dd8c396762c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2021 22:32:52 +0100 Subject: [PATCH] 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 +