diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f new file mode 100644 index 00000000..0171adf4 --- /dev/null +++ b/src/determinants/configurations.irp.f @@ -0,0 +1,623 @@ +use bitmasks +subroutine configuration_of_det(d,o,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Transforms a determinant to a configuration + ! + ! occ(:,1) : Single occupations + ! + ! occ(:,2) : Double occupations + ! + END_DOC + integer ,intent(in) :: Nint + integer(bit_kind),intent(in) :: d(Nint,2) + integer(bit_kind),intent(out) :: o(Nint,2) + + integer :: k + + do k=1,Nint + o(k,1) = ieor(d(k,1),d(k,2)) + o(k,2) = iand(d(k,1),d(k,2)) + enddo +end + + +subroutine configuration_to_dets_size(o,sze,n_alpha,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Number of possible determinants for a given configuration + END_DOC + integer ,intent(in) :: Nint, n_alpha + integer(bit_kind),intent(in) :: o(Nint,2) + integer, intent(out) :: sze + integer :: amax,bmax,k + double precision, external :: binom_func + + bmax = 0 + amax = n_alpha + do k=1,Nint + bmax += popcnt( o(k,1) ) + 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 + sze = int(binom_int(bmax, amax),4) +end + + +subroutine configuration_to_dets(o,d,sze,n_alpha,Nint) + use bitmasks + 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` + ! n_alpha : Number of $\alpha$ electrons + ! Nint : N_int + ! + ! Output: + ! d : determinants + ! + END_DOC + integer ,intent(in) :: Nint + integer ,intent(in) :: n_alpha ! Number of alpha electrons + integer ,intent(inout) :: sze ! Dimension of the output dets + integer(bit_kind),intent(in) :: o(Nint,2) ! Configurations + integer(bit_kind),intent(out) :: d(Nint,2,sze) ! Output determinants + + integer :: i, k, n, ispin, ispin2 + + ! Extract list of singly occupied MOs as (int,pos) pairs + ! ------------------------------------------------------ + + integer :: iint(2*n_alpha), ipos(2*n_alpha) + integer(bit_kind) :: v, t, tt, diff, v_prev + integer :: n_alpha_in_single + + n=0 + n_alpha_in_single = n_alpha + do i=1,Nint + v = o(i,1) + do while(v /= 0_bit_kind) + n = n+1 + iint(n) = i + ipos(n) = trailz(v) + v = iand(v,v-1) + enddo + n_alpha_in_single = n_alpha_in_single - popcnt( o(i,2) ) + enddo + + v = shiftl(1,n_alpha_in_single) - 1 + + ! Initialize first determinant + d(:,1,1) = o(:,2) + d(:,2,1) = o(:,2) + + do k=1,n_alpha_in_single + d(iint(k),1,1) = ibset( d(iint(k),1,1), ipos(k) ) + enddo + + do k=n_alpha_in_single+1,n + d(iint(k),2,1) = ibset( d(iint(k),2,1), ipos(k) ) + enddo + + sze = int(binom_int(n,n_alpha_in_single),4) + + if ( (shiftl(n_alpha_in_single,1) == n).and.n>0 ) then + + ! Time reversal symmetry + d(:,1,2) = d(:,2,1) + d(:,2,2) = d(:,1,1) + + do i=3,sze,2 + ! Generate next permutation with Anderson's algorithm + v_prev = v + t = ior(v,v-1) + tt = t+1 + v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) ) + + ! Find what has changed between v_prev and v + diff = ieor(v,v_prev) + + ! Initialize with previous determinant + d(:,1,i) = d(:,1,i-2) + d(:,2,i) = d(:,2,i-2) + + ! Swap bits only where they have changed from v_prev to v + do while (diff /= 0_bit_kind) + k = trailz(diff)+1 + if (btest(v,k-1)) then + d(iint(k),1,i) = ibset( d(iint(k),1,i), ipos(k) ) + d(iint(k),2,i) = ibclr( d(iint(k),2,i), ipos(k) ) + else + d(iint(k),1,i) = ibclr( d(iint(k),1,i), ipos(k) ) + d(iint(k),2,i) = ibset( d(iint(k),2,i), ipos(k) ) + endif + diff = iand(diff,diff-1_bit_kind) + enddo + + ! Time reversal symmetry + d(:,1,i+1) = d(:,2,i) + d(:,2,i+1) = d(:,1,i) + + enddo + + else + + do i=2,sze + ! Generate next permutation with Anderson's algorithm + v_prev = v + t = ior(v,v-1) + tt = t+1 + v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) ) + + ! Find what has changed between v_prev and v + diff = ieor(v,v_prev) + + ! Initialize with previous determinant + d(:,1,i) = d(:,1,i-1) + d(:,2,i) = d(:,2,i-1) + + ! Swap bits only where they have changed from v_prev to v + do while (diff /= 0_bit_kind) + k = trailz(diff)+1 + if (btest(v,k-1)) then + d(iint(k),1,i) = ibset( d(iint(k),1,i), ipos(k) ) + d(iint(k),2,i) = ibclr( d(iint(k),2,i), ipos(k) ) + else + d(iint(k),1,i) = ibclr( d(iint(k),1,i), ipos(k) ) + d(iint(k),2,i) = ibset( d(iint(k),2,i), ipos(k) ) + endif + diff = iand(diff,diff-1_bit_kind) + enddo + + enddo + + endif + +end + + + BEGIN_PROVIDER [ integer(bit_kind), psi_configuration, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_configuration ] + implicit none + BEGIN_DOC + ! Array of the configurations present in the wave function. + ! + ! psi_configuration(:,1,j) = j-th configuration of the wave function : represents all the single occupations + ! + ! psi_configuration(:,2,j) = j-th configuration of the wave function : represents all the double occupations + ! + ! The occ patterns are sorted by :c:func:`configuration_search_key` + END_DOC + integer :: i,j,k + + ! create + do i = 1, N_det + do k = 1, N_int + psi_configuration(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i)) + psi_configuration(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i)) + enddo + enddo + + ! Sort + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: configuration_search_key + integer(bit_kind), allocatable :: tmp_array(:,:,:) + logical,allocatable :: duplicate(:) + logical :: dup + + + allocate ( iorder(N_det), duplicate(N_det), bit_tmp(N_det), tmp_array(N_int,2,N_det) ) + + do i=1,N_det + iorder(i) = i + bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int) + enddo + + call i8sort(bit_tmp,iorder,N_det) + + + !$OMP PARALLEL DEFAULT(shared) PRIVATE(i,j,k,dup) + + !$OMP DO + do i=1,N_det + do k=1,N_int + tmp_array(k,1,i) = psi_configuration(k,1,iorder(i)) + tmp_array(k,2,i) = psi_configuration(k,2,iorder(i)) + enddo + duplicate(i) = .False. + enddo + !$OMP END DO + + ! Find duplicates + !$OMP DO + do i=1,N_det-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j+=1 + if (j>N_det) then + exit + endif + cycle + endif + dup = .True. + do k=1,N_int + 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 = j+1 + if (j>N_det) then + exit + endif + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + ! Copy filtered result + N_configuration=0 + do i=1,N_det + if (duplicate(i)) then + cycle + endif + N_configuration += 1 + do k=1,N_int + psi_configuration(k,1,N_configuration) = tmp_array(k,1,i) + psi_configuration(k,2,N_configuration) = tmp_array(k,2,i) + enddo + enddo + +!- Check +! print *, 'Checking for duplicates in configuration' +! do i=1,N_configuration +! do j=i+1,N_configuration +! duplicate(1) = .True. +! do k=1,N_int +! if (psi_configuration(k,1,i) /= psi_configuration(k,1,j)) then +! duplicate(1) = .False. +! exit +! endif +! if (psi_configuration(k,2,i) /= psi_configuration(k,2,j)) then +! duplicate(1) = .False. +! exit +! endif +! enddo +! if (duplicate(1)) then +! call debug_det(psi_configuration(1,1,i),N_int) +! call debug_det(psi_configuration(1,1,j),N_int) +! stop 'DUPLICATE' +! endif +! enddo +! enddo +! print *, 'No duplicates' +!- + deallocate(iorder,duplicate,bit_tmp,tmp_array) + +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 + ! Returns the index of the configuration for each determinant + END_DOC + integer :: i,j,k,r,l + integer*8 :: key, key2 + integer(bit_kind) :: occ(N_int,2) + logical :: found + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: configuration_search_key + + 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 + + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP PRIVATE(i,k,j,r,l,key,found,occ) + do i=1,N_det + do k = 1, N_int + occ(k,1) = ieor(psi_det(k,1,i),psi_det(k,2,i)) + occ(k,2) = iand(psi_det(k,1,i),psi_det(k,2,i)) + enddo + + 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 + 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 + endif + j = shiftr(r-l,1) + enddo + + enddo + !$OMP END PARALLEL DO + deallocate(bit_tmp) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, psi_configuration_Hii, (N_configuration) ] + implicit none + BEGIN_DOC + ! $\langle I|H|I \rangle$ where $|I\rangle$ is a configuration. + ! This is the minimum $H_{ii}$, where the $|i\rangle$ are the + ! determinants of $|I\rangle$. + END_DOC + integer :: j, i + + psi_configuration_Hii(:) = huge(1.d0) + do i=1,N_det + j = det_to_configuration(i) + psi_configuration_Hii(j) = min(psi_configuration_Hii(j), psi_det_Hii(i)) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, weight_configuration, (N_configuration,N_states) ] + implicit none + BEGIN_DOC + ! Weight of the configurations in the wave function + END_DOC + integer :: i,j,k + weight_configuration = 0.d0 + do i=1,N_det + j = det_to_configuration(i) + do k=1,N_states + weight_configuration(j,k) += psi_coef(i,k) * psi_coef(i,k) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, weight_configuration_average, (N_configuration) ] + implicit none + BEGIN_DOC + ! State-average weight of the configurations in the wave function + END_DOC + integer :: i,j,k + weight_configuration_average(:) = 0.d0 + do i=1,N_det + j = det_to_configuration(i) + do k=1,N_states + weight_configuration_average(j) += psi_coef(i,k) * psi_coef(i,k) * state_average_weight(k) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_configuration_sorted, (N_int,2,N_configuration) ] +&BEGIN_PROVIDER [ double precision, weight_configuration_average_sorted, (N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order, (N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order_reverse, (N_configuration) ] + implicit none + BEGIN_DOC + ! Configurations sorted by weight + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + allocate ( iorder(N_configuration) ) + do i=1,N_configuration + weight_configuration_average_sorted(i) = -weight_configuration_average(i) + iorder(i) = i + enddo + 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)) + enddo + psi_configuration_sorted_order(iorder(i)) = i + psi_configuration_sorted_order_reverse(i) = iorder(i) + weight_configuration_average_sorted(i) = -weight_configuration_average_sorted(i) + enddo + + deallocate(iorder) + +END_PROVIDER + + +subroutine make_s2_eigenfunction + implicit none + integer :: i,j,k + integer :: smax, s + integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:) + integer :: N_det_new, ithread, omp_get_thread_num + integer, parameter :: bufsze = 1000 + logical, external :: is_in_wavefunction + logical :: update + + update=.False. + call write_int(6,N_configuration,'Number of configurations') + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(N_configuration, psi_configuration, elec_alpha_num,N_int,update) & + !$OMP PRIVATE(s,ithread, d, det_buffer, smax, N_det_new,i,j,k) + N_det_new = 0 + call configuration_to_dets_size(psi_configuration(1,1,1),s,elec_alpha_num,N_int) + allocate (d(N_int,2,s+64), det_buffer(N_int,2,bufsze) ) + smax = s + ithread=0 + !$ ithread = omp_get_thread_num() + !$OMP DO SCHEDULE (dynamic,1000) + do i=1,N_configuration + call configuration_to_dets_size(psi_configuration(1,1,i),s,elec_alpha_num,N_int) + s += 1 + if (s > smax) then + deallocate(d) + allocate ( d(N_int,2,s+64) ) + smax = s + endif + call configuration_to_dets(psi_configuration(1,1,i),d,s,elec_alpha_num,N_int) + do j=1,s + if ( is_in_wavefunction(d(1,1,j), N_int) ) then + cycle + endif + update = .true. + N_det_new += 1 + det_buffer(:,:,N_det_new) = d(:,:,j) + if (N_det_new == bufsze) then + call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread) + N_det_new = 0 + endif + enddo + enddo + !$OMP END DO NOWAIT + + if (N_det_new > 0) then + call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,ithread) + endif + !$OMP BARRIER + deallocate(d,det_buffer) + !$OMP END PARALLEL + + if (update) then + call copy_H_apply_buffer_to_wf + TOUCH N_det psi_coef psi_det psi_configuration N_configuration + endif + call write_time(6) + +end + + + +BEGIN_PROVIDER [ integer, dominant_cfg, (N_states) ] + implicit none + BEGIN_DOC + ! Configuration of the determinants with the largest weight, for each state + END_DOC + integer :: k + do k=1,N_states + dominant_cfg(k) = det_to_configuration(dominant_det(k)) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer, N_dominant_dets_of_cfgs ] + implicit none + BEGIN_DOC + ! Number of determinants in all dominant determinants + END_DOC + integer :: k, sze + + N_dominant_dets_of_cfgs = 0 + do k=1,N_states + call configuration_to_dets_size( & + psi_configuration(1,1,dominant_cfg(k)), & + sze, elec_alpha_num, N_int) + N_dominant_dets_of_cfgs += sze + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_dets_of_cfgs) ] + implicit none + BEGIN_DOC + ! Configuration of the determinants with the largest weight, for each state + END_DOC + integer :: i,k,sze + i=1 + do k=1,N_states + sze = N_dominant_dets_of_cfgs + call configuration_to_dets( & + psi_configuration(1,1,dominant_cfg(k)), & + dominant_dets_of_cfgs(1,1,i), & + sze,elec_alpha_num,N_int) + 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. + ! + ! 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 +