diff --git a/config/gfortran.cfg b/config/gfortran.cfg index 694ef0df..c0aa875f 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -35,7 +35,7 @@ OPENMP : 1 ; Append OpenMP flags # -ffast-math and the Fortran-specific # -fno-protect-parens and -fstack-arrays. [OPT] -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Profiling flags ################# diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 72084241..03663eea 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -51,7 +51,7 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [DEBUG] -FCFLAGS : -g -pedantic -msse4.2 +FCFLAGS : -g -msse4.2 # OpenMP flags ################# diff --git a/config/ifort.cfg b/config/ifort.cfg index cb6dc1ef..585e4744 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -32,14 +32,14 @@ OPENMP : 1 ; Append OpenMP flags # [OPT] FC : -traceback -FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g -traceback +FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g # Profiling flags ################# # [PROFILE] FC : -p -g -traceback -FCFLAGS : -xSSE4.2 -O2 -ip -ftz +FCFLAGS : -xSSE4.2 -O2 -ip -ftz # Debugging flags ################# @@ -52,7 +52,7 @@ FCFLAGS : -xSSE4.2 -O2 -ip -ftz # [DEBUG] FC : -g -traceback -FCFLAGS : -xSSE2 -C +FCFLAGS : -xSSE2 -C -fpe0 IRPF90_FLAGS : --openmp # OpenMP flags diff --git a/data/pseudo/tn_df b/data/pseudo/tn_df index 988312b0..79ebf8f5 100644 --- a/data/pseudo/tn_df +++ b/data/pseudo/tn_df @@ -780,7 +780,7 @@ Ar GEN 10 2 -1386.79918148 2 4.23753203 1350.57102634 2 6.12344921 -Ag GEN 36 2 +Ag GEN 36 2 6 11.00000000 1 7.02317516 178.71479273 2 1.36779344 diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 449d262c..0874e584 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -31,11 +31,11 @@ s.set_perturbation("epstein_nesbet_2x2") s.unset_openmp() print s -#s = H_apply_zmq("mrcc_PT2") -#s.energy = "ci_electronic_energy_dressed" -#s.set_perturbation("epstein_nesbet_2x2") -#s.unset_openmp() -#print s +s = H_apply_zmq("mrcepa_PT2") +s.energy = "psi_ref_energy_diagonalized" +s.set_perturbation("epstein_nesbet_2x2") +s.unset_openmp() +print s END_SHELL diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 354ea389..66f4975a 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -47,7 +47,7 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i !$OMP END DO !$OMP DO SCHEDULE(guided) do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(i,istate) + H_jj(idx_ref(i)) += delta_ii(istate,i) enddo !$OMP END DO !$OMP END PARALLEL @@ -269,7 +269,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin to_print(2,k) = residual_norm(k) enddo - write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st) + write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st) call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) if (converged) then exit @@ -487,8 +487,8 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) i = idx_ref(ii) do jj = 1, n_det_non_ref j = idx_non_ref(jj) - vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j) - vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i) + vt (i) = vt (i) + delta_ij(istate,jj,ii)*u_0(j) + vt (j) = vt (j) + delta_ij(istate,jj,ii)*u_0(i) enddo enddo !$OMP END DO diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 1c2e8b74..412c52e2 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -51,9 +51,9 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:) integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) integer :: mobiles(2), smallerlist + logical, external :: is_generable - - + print *, i_generator leng = max(N_det_generators, N_det_non_ref) allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref)) @@ -69,7 +69,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge allocate( microlist(Nint,2,N_minilist*4), & idx_microlist(N_minilist*4)) - if(key_mask(1,1) /= 0) then + if(key_mask(1,1) /= 0_8) then call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) else @@ -87,6 +87,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge ! |alpha> if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint) if(N_minilist == 0) return @@ -117,8 +118,18 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge - do i_alpha=1,N_tq +! ok = .false. +! do i=N_det_generators, 1, -1 +! if(is_generable(psi_det_generators(1,1,i), tq(1,1,i_alpha), Nint)) then +! ok = .true. +! exit +! end if +! end do +! if(.not. ok) then +! cycle +! end if + if(key_mask(1,1) /= 0) then call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) @@ -138,37 +149,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge do j=1,idx_alpha(0) idx_alpha(j) = idx_microlist_zero(idx_alpha(j)) end do - - -! i = 1 -! j = 2 -! do j = 2, idx_alpha_tmp(0) -! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit -! end do -! -! m = j -! -! idx_alpha(0) = idx_alpha_tmp(0) -! -! do l = 1, idx_alpha(0) -! if(j > idx_alpha_tmp(0)) then -! k = i -! i += 1 -! else if(i >= m) then -! k = j -! j += 1 -! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then -! k = i -! i += 1 -! else -! k = j -! j += 1 -! end if -! ! k=l -! idx_alpha(l) = idx_alpha_tmp(k) -! degree_alpha(l) = degree_alpha_tmp(k) -! end do -! else call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) do j=1,idx_alpha(0) @@ -177,12 +157,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge end if -! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) -! do j=1,idx_alpha(0) -! idx_alpha(j) = idx_miniList(idx_alpha(j)) -! end do - !print *, idx_alpha(:idx_alpha(0)) - do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd)) @@ -285,33 +259,31 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge enddo enddo call omp_set_lock( psi_ref_lock(i_I) ) + + do i_state=1,Nstates if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) - delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) enddo else - delta_ii_(i_state,i_I) = 0.d0 + !delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) enddo endif enddo call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo - !deallocate (dIa_hla,hij_cache) - !deallocate(miniList, idx_miniList) + deallocate (dIa_hla,hij_cache) + deallocate(miniList, idx_miniList) end - - - - subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) use bitmasks @@ -360,7 +332,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq endif enddo if (good) then - if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then N_tq += 1 do k=1,N_int tq(k,1,N_tq) = det_buffer(k,1,i) @@ -437,7 +409,7 @@ subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,N endif enddo if (good) then - if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then N_tq += 1 do k=1,N_int tq(k,1,N_tq) = det_buffer(k,1,i) diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index 50343fdb..d356e4b9 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -1,60 +1,3 @@ -subroutine run_mrcc - implicit none - call set_generators_bitmasks_as_holes_and_particles - call mrcc_iterations -end - -subroutine mrcc_iterations - implicit none - - integer :: i,j - - double precision :: E_new, E_old, delta_e - integer :: iteration,i_oscillations - double precision :: E_past(4), lambda - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - j = 1 - i_oscillations = 0 - lambda = 1.d0 - do while (delta_E > 1.d-7) - iteration += 1 - print *, '===========================' - print *, 'MRCC Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed) - call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed(lambda) - E_new = sum(ci_energy_dressed) - delta_E = dabs(E_new - E_old) -! if (E_new > E_old) then -! lambda = lambda * 0.7d0 -! else -! lambda = min(1.d0, lambda * 1.1d0) -! endif -! print *, 'energy lambda ', lambda - E_past(j) = E_new - j +=1 - call save_wavefunction - if (iteration > 200) then - exit - endif - print*,'------------' - print*,'VECTOR' - do i = 1, N_det_ref - print*,'' - print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) - print*,'delta_ii(i,1) = ',delta_ii(i,1) - enddo - print*,'------------' - enddo - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - call save_wavefunction - -end subroutine set_generators_bitmasks_as_holes_and_particles implicit none @@ -81,7 +24,4 @@ subroutine set_generators_bitmasks_as_holes_and_particles enddo enddo touch generators_bitmask - - - end diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 4ac48602..6c2eb133 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,5 +1,13 @@ - BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] +use bitmasks + + BEGIN_PROVIDER [ integer, mrmode ] + mrmode = 0 +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] +&BEGIN_PROVIDER [ integer, lambda_mrcc_pt3, (0:psi_det_size) ] implicit none BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) @@ -8,48 +16,51 @@ double precision :: ihpsi_current(N_states) integer :: i_pert_count double precision :: hii, lambda_pert - integer :: N_lambda_mrcc_pt2 - + integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3 + i_pert_count = 0 lambda_mrcc = 0.d0 N_lambda_mrcc_pt2 = 0 + N_lambda_mrcc_pt3 = 0 lambda_mrcc_pt2(0) = 0 + lambda_mrcc_pt3(0) = 0 - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& - size(psi_ref_coef,1), N_states,ihpsi_current) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 + do i=1,N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& + size(psi_ref_coef,1), N_states,ihpsi_current) + call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif + lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) + lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 + if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then + N_lambda_mrcc_pt2 += 1 + lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i endif - lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) ) - lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then - i_pert_count += 1 - lambda_mrcc(k,i) = 0.d0 - if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then - N_lambda_mrcc_pt2 += 1 - lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i - endif + else + if (lambda_mrcc_pt3(N_lambda_mrcc_pt3) /= i) then + N_lambda_mrcc_pt3 += 1 + lambda_mrcc_pt3(N_lambda_mrcc_pt3) = i endif - enddo + endif enddo - lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 - + enddo + lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 + lambda_mrcc_pt3(0) = N_lambda_mrcc_pt3 print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of ignored determinants = ',i_pert_count print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) print*,'lambda max = ',maxval(dabs(lambda_mrcc)) + print*,'Number of ignored determinants = ',i_pert_count END_PROVIDER - - - - BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] implicit none BEGIN_DOC @@ -74,7 +85,9 @@ END_PROVIDER delta_ij = 0.d0 delta_ii = 0.d0 call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref) + END_PROVIDER + BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] implicit none @@ -201,3 +214,763 @@ subroutine diagonalize_CI_dressed(lambda) end + +logical function is_generable(det1, det2, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) + integer :: degree, f, exc(0:2, 2, 2), t + integer*2 :: h1, h2, p1, p2, s1, s2 + integer, external :: searchExc + logical, external :: excEq + double precision :: phase + + is_generable = .false. + call get_excitation(det1, det2, exc, degree, phase, Nint) + if(degree == -1) return + if(degree == 0) then + is_generable = .true. + return + end if + if(degree > 2) stop "?22??" + + call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) + + if(degree == 1) then + h2 = h1 + p2 = p1 + s2 = s1 + h1 = 0 + p1 = 0 + s1 = 0 + end if + + if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then + f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) + else + f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) + end if + if(f == -1) return + + if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then + f = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) + else + f = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) + end if + + if(f /= -1) is_generable = .true. +end function + + + +integer function searchDet(dets, det, n, Nint) + implicit none + use bitmasks + + integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2) + integer, intent(in) :: nint, n + integer :: l, h, c + integer, external :: detCmp + logical, external :: detEq + + l = 1 + h = n + do while(.true.) + searchDet = (l+h)/2 + c = detCmp(dets(1,1,searchDet), det(1,1), Nint) + if(c == 0) return + if(c == 1) then + h = searchDet-1 + else + l = searchDet+1 + end if + if(l > h) then + searchDet = -1 + return + end if + + end do +end function + + +integer function unsortedSearchDet(dets, det, n, Nint) + implicit none + use bitmasks + + integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2) + integer, intent(in) :: nint, n + integer :: l, h, c + integer, external :: detCmp + logical, external :: detEq + + do l=1, n + if(detEq(det, dets(1,1,l), N_int)) then + unsortedSearchDet = l + return + end if + end do + unsortedSearchDet = -1 +end function + + +integer function searchExc(excs, exc, n) + implicit none + use bitmasks + + integer, intent(in) :: n + integer*2,intent(in) :: excs(4,n), exc(4) + integer :: l, h, c + integer, external :: excCmp + logical, external :: excEq + + l = 1 + h = n + do + searchExc = (l+h)/2 + c = excCmp(excs(1,searchExc), exc(1)) + if(c == 0) return + if(c == 1) then + h = searchExc-1 + else + l = searchExc+1 + end if + if(l > h) then + searchExc = -1 + return + end if + end do +end function + + +subroutine sort_det(key, idx, N_key, Nint) + implicit none + + + integer, intent(in) :: Nint, N_key + integer(8),intent(inout) :: key(Nint,2,N_key) + integer,intent(inout) :: idx(N_key) + integer(8) :: tmp(Nint, 2) + integer :: tmpidx,i,ni + + do i=1,N_key + idx(i) = i + end do + + do i=N_key/2,1,-1 + call tamiser(key, idx, i, N_key, Nint, N_key) + end do + + do i=N_key,2,-1 + do ni=1,Nint + tmp(ni,1) = key(ni,1,i) + tmp(ni,2) = key(ni,2,i) + key(ni,1,i) = key(ni,1,1) + key(ni,2,i) = key(ni,2,1) + key(ni,1,1) = tmp(ni,1) + key(ni,2,1) = tmp(ni,2) + enddo + + tmpidx = idx(i) + idx(i) = idx(1) + idx(1) = tmpidx + call tamiser(key, idx, 1, i-1, Nint, N_key) + end do +end subroutine + + +subroutine sort_exc(key, N_key) + implicit none + + + integer, intent(in) :: N_key + integer*2,intent(inout) :: key(4,N_key) + integer*2 :: tmp(4) + integer :: i,ni + + + do i=N_key/2,1,-1 + call tamise_exc(key, i, N_key, N_key) + end do + + do i=N_key,2,-1 + do ni=1,4 + tmp(ni) = key(ni,i) + key(ni,i) = key(ni,1) + key(ni,1) = tmp(ni) + enddo + + call tamise_exc(key, 1, i-1, N_key) + end do +end subroutine + + +logical function exc_inf(exc1, exc2) + implicit none + integer*2,intent(in) :: exc1(4), exc2(4) + integer :: i + exc_inf = .false. + do i=1,4 + if(exc1(i) < exc2(i)) then + exc_inf = .true. + return + else if(exc1(i) > exc2(i)) then + return + end if + end do +end function + + +subroutine tamise_exc(key, no, n, N_key) + use bitmasks + implicit none + + BEGIN_DOC +! Uncodumented : TODO + END_DOC + integer,intent(in) :: no, n, N_key + integer*2,intent(inout) :: key(4, N_key) + integer :: k,j + integer*2 :: tmp(4) + logical :: exc_inf + integer :: ni + + k = no + j = 2*k + do while(j <= n) + if(j < n) then + if (exc_inf(key(1,j), key(1,j+1))) then + j = j+1 + endif + endif + if(exc_inf(key(1,k), key(1,j))) then + do ni=1,4 + tmp(ni) = key(ni,k) + key(ni,k) = key(ni,j) + key(ni,j) = tmp(ni) + enddo + k = j + j = k+k + else + return + endif + enddo +end subroutine + + +subroutine dec_exc(exc, h1, h2, p1, p2) + implicit none + integer :: exc(0:2,2,2), s1, s2, degree + integer*2, intent(out) :: h1, h2, p1, p2 + + degree = exc(0,1,1) + exc(0,1,2) + + h1 = 0 + h2 = 0 + p1 = 0 + p2 = 0 + + if(degree == 0) return + + call decode_exc_int2(exc, degree, h1, p1, h2, p2, s1, s2) + + h1 += mo_tot_num * (s1-1) + p1 += mo_tot_num * (s1-1) + + if(degree == 2) then + h2 += mo_tot_num * (s2-1) + p2 += mo_tot_num * (s2-1) + if(h1 > h2) then + s1 = h1 + h1 = h2 + h2 = s1 + end if + if(p1 > p2) then + s1 = p1 + p1 = p2 + p2 = s1 + end if + else + h2 = h1 + p2 = p1 + p1 = 0 + h1 = 0 + end if +end subroutine + + + BEGIN_PROVIDER [ integer, N_hh_exists ] +&BEGIN_PROVIDER [ integer, N_pp_exists ] +&BEGIN_PROVIDER [ integer, N_ex_exists ] + implicit none + integer :: exc(0:2, 2, 2), degree, n, on, s, l, i + integer*2 :: h1, h2, p1, p2 + double precision :: phase + logical,allocatable :: hh(:,:) , pp(:,:) + + allocate(hh(0:mo_tot_num*2, 0:mo_tot_num*2)) + allocate(pp(0:mo_tot_num*2, 0:mo_tot_num*2)) + hh = .false. + pp = .false. + N_hh_exists = 0 + N_pp_exists = 0 + N_ex_exists = 0 + + n = 0 + do i=1, N_det_ref + do l=1, N_det_non_ref + call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int) + if(degree == -1) cycle + call dec_exc(exc, h1, h2, p1, p2) + N_ex_exists += 1 + if(.not. hh(h1,h2)) N_hh_exists = N_hh_exists + 1 + if(.not. pp(p1,p2)) N_pp_exists = N_pp_exists + 1 + hh(h1,h2) = .true. + pp(p1,p2) = .true. + end do + end do + N_pp_exists = min(N_ex_exists, N_pp_exists * N_hh_exists) +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref_sorted, (N_int, 2, N_det_non_ref) ] +&BEGIN_PROVIDER [ integer, psi_non_ref_sorted_idx, (N_det_non_ref) ] + implicit none + psi_non_ref_sorted = psi_non_ref + call sort_det(psi_non_ref_sorted, psi_non_ref_sorted_idx, N_det_non_ref, N_int) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1, N_states) ] + implicit none + logical :: ok + integer :: i, j, k, s, II, pp, hh, ind, wk, nex, a_col, at_row + integer, external :: searchDet, unsortedSearchDet + integer(bit_kind) :: myDet(N_int, 2), myMask(N_int, 2) + integer :: N, INFO, AtA_size, r1, r2 + double precision , allocatable:: B(:), AtB(:), AtA_val(:), A_val(:,:), x(:), x_new(:), A_val_mwen(:) + double precision :: t, norm, cx + integer, allocatable :: A_ind(:,:), lref(:), AtA_ind(:), A_ind_mwen(:), col_shortcut(:), N_col(:) + + + + nex = hh_shortcut(hh_shortcut(0)+1)-1 + print *, "TI", nex, N_det_non_ref + allocate(A_ind(N_det_ref+1, nex), A_val(N_det_ref+1, nex)) + allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL ? !!!!!!!! + allocate(x(nex), AtB(nex)) + allocate(A_val_mwen(nex), A_ind_mwen(nex)) + allocate(N_col(nex), col_shortcut(nex), B(N_det_non_ref)) + allocate (x_new(nex)) + + do s = 1, N_states + + A_val = 0d0 + A_ind = 0 + AtA_ind = 0 + AtA_val = 0d0 + x = 0d0 + AtB = 0d0 + A_val_mwen = 0d0 + A_ind_mwen = 0 + N_col = 0 + col_shortcut = 0 + B = 0d0 + x_new = 0d0 + + !$OMP PARALLEL DO schedule(static,10) default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind) & + !$OMP shared(s, hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, psi_non_ref_sorted_idx, psi_ref, N_det_ref) & + !$OMP private(lref, pp, II, ok, myMask, myDet, ind, wk) + do hh = 1, hh_shortcut(0) + do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 + allocate(lref(N_det_non_ref)) + lref = 0 + do II = 1, N_det_ref + call apply_hole(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int) + if(.not. ok) cycle + call apply_particle(myMask, pp_exists(1, pp), myDet, ok, N_int) + if(.not. ok) cycle + ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int) + if(ind /= -1) then + lref(psi_non_ref_sorted_idx(ind)) = II + end if + end do + wk = 0 + do i=1, N_det_non_ref + if(lref(i) /= 0) then + wk += 1 + A_val(wk, pp) = psi_ref_coef(lref(i), s) + A_ind(wk, pp) = i + end if + end do + deallocate(lref) + end do + end do + !$OMP END PARALLEL DO + + AtB = 0d0 + AtA_size = 0 + wk = 0 + col_shortcut = 0 + N_col = 0 + !$OMP PARALLEL DO schedule(dynamic, 100) default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref) & + !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen) & + !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s) + do at_row = 1, nex + wk = 0 + if(mod(at_row, 10000) == 0) print *, "AtA", at_row, "/", nex + do i=1,N_det_ref + if(A_ind(i, at_row) == 0) exit + AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_row), s) * A_val(i, at_row) + end do + do a_col = 1, nex + t = 0d0 + r1 = 1 + r2 = 1 + do while(A_ind(r1, at_row) * A_ind(r2, a_col) /= 0) + if(A_ind(r1, at_row) < A_ind(r2, a_col)) then + r1 += 1 + else if(A_ind(r1, at_row) > A_ind(r2, a_col)) then + r2 += 1 + else + t = t - A_val(r1, at_row) * A_val(r2, a_col) + r1 += 1 + r2 += 1 + end if + end do + + if(a_col == at_row) then + t = (t + 1d0) + end if + if(t /= 0d0) then + wk += 1 + A_ind_mwen(wk) = a_col + A_val_mwen(wk) = t + end if + end do + + if(wk /= 0) then + !$OMP CRITICAL + col_shortcut(at_row) = AtA_size+1 + N_col(at_row) = wk + AtA_ind(AtA_size+1:AtA_size+wk) = A_ind_mwen(:wk) + AtA_val(AtA_size+1:AtA_size+wk) = A_val_mwen(:wk) + AtA_size += wk + !$OMP END CRITICAL + end if + end do + + x = AtB + if(AtA_size > size(AtA_val)) stop "SIZA" + print *, "ATA SIZE", ata_size + integer :: iproc, omp_get_thread_num + iproc = omp_get_thread_num() + do i=1,nex + x_new(i) = 0.D0 + enddo + + do k=0,100000 + !$OMP PARALLEL DO default(shared) + do i=1,nex + x_new(i) = AtB(i) + enddo + + !$OMP PARALLEL DO default(shared) private(cx, i) + do a_col = 1, nex + cx = 0d0 + do i=col_shortcut(a_col), col_shortcut(a_col) + N_col(a_col) - 1 + cx += x(AtA_ind(i)) * AtA_val(i) + end do + x_new(a_col) += cx + end do + !$OMP END PARALLEL DO + double precision :: norm_cas + norm_cas = 0d0 + do i = 1, N_det_ref + norm_cas += psi_ref_coef(i,s)**2 + end do + + norm = 0d0 + t = 0d0 + + do j=1, size(X) + t = t + X_new(j) * X_new(j) + end do + + + t = (1d0 / norm_cas - 1d0) / t + x_new = x_new * sqrt(t) + + do j=1, size(X) + norm += (X_new(j) - X(j))**2 + x(j) = x_new(j) + end do + + + if(mod(k, 100) == 0) then + print *, "residu ", k, norm, "norm t", sqrt(t) + end if + + if(norm < 1d-16) exit + end do + print *, "CONVERGENCE : ", norm + + dIj_unique(:size(X), s) = X(:) + + + end do + + + print *, "done" +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, dij, (N_det_ref, N_det_non_ref, N_states) ] + integer :: s,i,j + print *, "computing amplitudes..." + do s=1, N_states + do i=1, N_det_non_ref + do j=1, N_det_ref + dij(j, i, s) = get_dij_index(j, i, s, N_int) + end do + end do + end do + print *, "done computing amplitudes" +END_PROVIDER + + + + +double precision function get_dij_index(II, i, s, Nint) + integer, intent(in) :: II, i, s, Nint + double precision, external :: get_dij + double precision :: HIi + + if(lambda_type == 0) then + get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) + else + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) + get_dij_index = HIi * lambda_mrcc(s, i) + end if +end function + + +double precision function get_dij(det1, det2, s, Nint) + use bitmasks + implicit none + integer, intent(in) :: s, Nint + integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) + integer :: degree, f, exc(0:2, 2, 2), t + integer*2 :: h1, h2, p1, p2, s1, s2 + integer, external :: searchExc + logical, external :: excEq + double precision :: phase + + get_dij = 0d0 + call get_excitation(det1, det2, exc, degree, phase, Nint) + if(degree == -1) return + if(degree == 0) then + stop "get_dij" + end if + + call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) + + if(degree == 1) then + h2 = h1 + p2 = p1 + s2 = s1 + h1 = 0 + p1 = 0 + s1 = 0 + end if + + if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then + f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) + else + f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) + end if + if(f == -1) return + + if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then + t = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) + else + t = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) + end if + + if(t /= -1) then + get_dij = dIj_unique(t - 1 + hh_shortcut(f), s) + end if +end function + + + BEGIN_PROVIDER [ integer*2, hh_exists, (4, N_hh_exists) ] +&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_hh_exists + 1) ] +&BEGIN_PROVIDER [ integer*2, pp_exists, (4, N_pp_exists) ] + implicit none + integer*2,allocatable :: num(:,:) + integer :: exc(0:2, 2, 2), degree, n, on, s, l, i + integer*2 :: h1, h2, p1, p2 + double precision :: phase + logical, external :: excEq + + allocate(num(4, N_ex_exists+1)) + + hh_shortcut = 0 + hh_exists = 0 + pp_exists = 0 + num = 0 + + n = 0 + do i=1, N_det_ref + do l=1, N_det_non_ref + call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int) + if(degree == -1) cycle + call dec_exc(exc, h1, h2, p1, p2) + n += 1 + num(:, n) = (/h1, h2, p1, p2/) + end do + end do + + call sort_exc(num, n) + + hh_shortcut(0) = 1 + hh_shortcut(1) = 1 + hh_exists(:,1) = (/1_2, num(1,1), 1_2, num(2,1)/) + pp_exists(:,1) = (/1_2, num(3,1), 1_2, num(4,1)/) + s = 1 + do i=2,n + if(.not. excEq(num(1,i), num(1,s))) then + s += 1 + num(:, s) = num(:, i) + pp_exists(:,s) = (/1_2, num(3,s), 1_2, num(4,s)/) + if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. & + hh_exists(4, hh_shortcut(0)) /= num(2,s)) then + hh_shortcut(0) += 1 + hh_shortcut(hh_shortcut(0)) = s + hh_exists(:,hh_shortcut(0)) = (/1_2, num(1,s), 1_2, num(2,s)/) + end if + end if + end do + hh_shortcut(hh_shortcut(0)+1) = s+1 + + do s=2,4,2 + do i=1,hh_shortcut(0) + if(hh_exists(s, i) == 0) then + hh_exists(s-1, i) = 0 + else if(hh_exists(s, i) > mo_tot_num) then + hh_exists(s, i) -= mo_tot_num + hh_exists(s-1, i) = 2 + end if + end do + + do i=1,hh_shortcut(hh_shortcut(0)+1)-1 + if(pp_exists(s, i) == 0) then + pp_exists(s-1, i) = 0 + else if(pp_exists(s, i) > mo_tot_num) then + pp_exists(s, i) -= mo_tot_num + pp_exists(s-1, i) = 2 + end if + end do + end do +END_PROVIDER + + +logical function excEq(exc1, exc2) + implicit none + integer*2, intent(in) :: exc1(4), exc2(4) + integer :: i + excEq = .false. + do i=1, 4 + if(exc1(i) /= exc2(i)) return + end do + excEq = .true. +end function + + +integer function excCmp(exc1, exc2) + implicit none + integer*2, intent(in) :: exc1(4), exc2(4) + integer :: i + excCmp = 0 + do i=1, 4 + if(exc1(i) > exc2(i)) then + excCmp = 1 + return + else if(exc1(i) < exc2(i)) then + excCmp = -1 + return + end if + end do +end function + + +subroutine apply_hole(det, exc, res, ok, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer*2, intent(in) :: exc(4) + integer*2 :: s1, s2, h1, h2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos + + ok = .false. + s1 = exc(1) + h1 = exc(2) + s2 = exc(3) + h2 = exc(4) + res = det + + if(h1 /= 0) then + ii = (h1-1)/bit_kind_size + 1 + pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s1) = ibclr(res(ii, s1), pos) + end if + + ii = (h2-1)/bit_kind_size + 1 + pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s2) = ibclr(res(ii, s2), pos) + + + ok = .true. +end subroutine + + +subroutine apply_particle(det, exc, res, ok, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer*2, intent(in) :: exc(4) + integer*2 :: s1, s2, p1, p2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos + + ok = .false. + s1 = exc(1) + p1 = exc(2) + s2 = exc(3) + p2 = exc(4) + res = det + + if(p1 /= 0) then + ii = (p1-1)/bit_kind_size + 1 + pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) + end if + + ii = (p2-1)/bit_kind_size + 1 + pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s2) = ibset(res(ii, s2), pos) + + + ok = .true. +end subroutine + diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg new file mode 100644 index 00000000..9979f537 --- /dev/null +++ b/plugins/mrcepa0/EZFIO.cfg @@ -0,0 +1,5 @@ +[lambda_type] +type: Strictly_positive_int +doc: lambda type ( 0 = none, 1 = last version ) +interface: ezfio,provider,ocaml +default: 0 diff --git a/plugins/mrcepa0/NEEDED_CHILDREN_MODULES b/plugins/mrcepa0/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..8b6c5a18 --- /dev/null +++ b/plugins/mrcepa0/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ diff --git a/plugins/mrcepa0/README.rst b/plugins/mrcepa0/README.rst new file mode 100644 index 00000000..997d005e --- /dev/null +++ b/plugins/mrcepa0/README.rst @@ -0,0 +1,12 @@ +======= +mrcepa0 +======= + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f new file mode 100644 index 00000000..3a91f42e --- /dev/null +++ b/plugins/mrcepa0/dressing.irp.f @@ -0,0 +1,1004 @@ +use bitmasks + + + + BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] + use bitmasks + implicit none + integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc + integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2) + integer(bit_kind),allocatable :: buf(:,:,:) + logical :: ok + logical, external :: detEq + + delta_ij_mrcc = 0d0 + delta_ii_mrcc = 0d0 + print *, "Dij", dij(1,1,1) + provide hh_shortcut psi_det_size! lambda_mrcc + !$OMP PARALLEL DO default(none) schedule(dynamic) & + !$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) & + !$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc) & + !$OMP private(h, n, mask, omask, buf, ok, iproc) + do gen= 1, N_det_generators + allocate(buf(N_int, 2, N_det_non_ref)) + iproc = omp_get_thread_num() + 1 + if(mod(gen, 10) == 0) print *, "mrcc ", gen, "/", N_det_generators + do h=1, hh_shortcut(0) + call apply_hole(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int) + if(.not. ok) cycle + omask = 0_bit_kind + if(hh_exists(1, h) /= 0) omask = mask + n = 1 + do p=hh_shortcut(h), hh_shortcut(h+1)-1 + call apply_particle(mask, pp_exists(1, p), buf(1,1,n), ok, N_int) + if(ok) n = n + 1 + if(n > N_det_non_ref) stop "MRCC..." + end do + n = n - 1 + + if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask) + + end do + deallocate(buf) + end do + !$OMP END PARALLEL DO +END_PROVIDER + + +! subroutine blit(b1, b2) +! double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref) +! b1 = b1 + b2 +! end subroutine + + +subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) + double precision, intent(inout) :: delta_ii_(N_states,N_det_ref) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l,m + integer,allocatable :: idx_alpha(:), degree_alpha(:) + logical :: good, fullMatch + + integer(bit_kind),allocatable :: tq(:,:,:) + integer :: N_tq, c_ref ,degree + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hla(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer :: iint, ipos + integer :: i_state, k_sd, l_sd, i_I, i_alpha + + integer(bit_kind),allocatable :: miniList(:,:,:) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer,allocatable :: idx_miniList(:) + integer :: N_miniList, ni, leng + double precision, allocatable :: hij_cache(:) + + integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:) + integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) + integer :: mobiles(2), smallerlist + logical, external :: detEq, is_generable + !double precision, external :: get_dij, get_dij_index + + + leng = max(N_det_generators, N_det_non_ref) + allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref)) + allocate(idx_alpha(0:psi_det_size), degree_alpha(psi_det_size)) + !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) + +! if(fullMatch) then +! return +! end if + + allocate(ptr_microlist(0:mo_tot_num*2+1), & + N_microlist(0:mo_tot_num*2) ) + allocate( microlist(Nint,2,N_minilist*4), & + idx_microlist(N_minilist*4)) + + if(key_mask(1,1) /= 0) then + call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) + call filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + else + call filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + end if + + + + deallocate(microlist, idx_microlist) + + allocate (dIa_hla(N_states,N_det_non_ref)) + + ! |I> + + ! |alpha> + + if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint) + if(N_minilist == 0) return + + + if(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!! + allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist)) + + allocate( microlist(Nint,2,N_minilist*4), & + idx_microlist(N_minilist*4)) + call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) + + + do i=0,mo_tot_num*2 + do k=ptr_microlist(i),ptr_microlist(i+1)-1 + idx_microlist(k) = idx_minilist(idx_microlist(k)) + end do + end do + + do l=1,N_microlist(0) + do k=1,Nint + microlist_zero(k,1,l) = microlist(k,1,l) + microlist_zero(k,2,l) = microlist(k,2,l) + enddo + idx_microlist_zero(l) = idx_microlist(l) + enddo + end if + end if + + + do i_alpha=1,N_tq + if(key_mask(1,1) /= 0) then + call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) + + if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then + smallerlist = mobiles(1) + else + smallerlist = mobiles(2) + end if + + + do l=0,N_microlist(smallerlist)-1 + microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l) + idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l) + end do + + call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_microlist_zero(idx_alpha(j)) + end do + + else + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + end do + end if + + + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd)) + enddo + ! |I> + do i_I=1,N_det_ref + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint) + if (degree > 4) then + cycle + endif + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + + ! |alpha> + do k_sd=1,idx_alpha(0) + ! Loop if lambda == 0 + logical :: loop +! loop = .True. +! do i_state=1,N_states +! if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then +! loop = .False. +! exit +! endif +! enddo +! if (loop) then +! cycle +! endif + + call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) + if (degree > 2) then + cycle + endif + + ! + ! + !hIk = hij_mrcc(idx_alpha(k_sd),i_I) + ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) + + do i_state=1,N_states + dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state) + !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state) + enddo + + + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do k=1,N_int + tmp_det(k,1) = psi_ref(k,1,i_I) + tmp_det(k,2) = psi_ref(k,2,i_I) + enddo + logical :: ok + call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) + if(.not. ok) cycle + + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) + if (degree == 0) then + +! loop = .True. +! do i_state=1,N_states +! if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then +! loop = .False. +! exit +! endif +! enddo + loop = .false. + if (.not.loop) then + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + hIl = hij_mrcc(idx_alpha(l_sd),i_I) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2 + !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 + enddo + endif + + exit + endif + enddo + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,N_states + ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) + enddo + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + hla = hij_cache(k_sd) +! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla) + do i_state=1,N_states + dIa_hla(i_state,k_sd) = dIa(i_state) * hla + enddo + enddo + call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + enddo + else + delta_ii_(i_state,i_I) = 0.d0 + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + enddo + endif + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) + enddo + enddo + deallocate (dIa_hla,hij_cache) + deallocate(miniList, idx_miniList) +end + + + + + BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ] + use bitmasks + implicit none + integer :: i, j, i_state + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + + do i_state = 1, N_states + if(mrmode == 3) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) + end do + end do +! +! do i = 1, N_det_ref +! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) +! do j = 1, N_det_non_ref +! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) +! end do +! end do + else if(mrmode == 2) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_old(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) + end do + end do + else if(mrmode == 1) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) + end do + end do + else + stop "invalid mrmode" + end if + end do +END_PROVIDER + + +BEGIN_PROVIDER [ integer, HP, (2,N_det_non_ref) ] + integer :: i + do i=1,N_det_non_ref + call getHP(psi_non_ref(1,1,i), HP(1,i), HP(2,i), N_int) + end do +END_PROVIDER + + BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ] +&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_int,2,N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_int,2,N_det_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0, (N_int,2,N_det_non_ref) ] +&BEGIN_PROVIDER [ integer, nlink, (N_det_ref) ] +&BEGIN_PROVIDER [ integer, linked, (N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ integer, blokMwen, (N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, searchance, (N_det_ref) ] +&BEGIN_PROVIDER [ integer, child_num, (N_det_non_ref,N_det_ref) ] + + use bitmasks + implicit none + + integer(bit_kind),allocatable :: det_noactive(:,:,:) + integer, allocatable :: shortcut(:), idx(:) + integer(bit_kind) :: nonactive_sorb(N_int,2), det(N_int, 2) + integer i, II, j, k, n, ni, blok, degree + logical, external :: detEq + + allocate(det_noactive(N_int, 2, N_det_non_ref)) + allocate(idx(N_det_non_ref), shortcut(0:N_det_non_ref+1)) + print *, "pre start" + active_sorb(:,:) = 0_8 + nonactive_sorb(:,:) = not(0_8) + + if(N_det_ref > 1) then + do i=1, N_det_ref + do k=1, N_int + active_sorb(k,1) = ior(psi_ref(k,1,i), active_sorb(k,1)) + active_sorb(k,2) = ior(psi_ref(k,2,i), active_sorb(k,2)) + nonactive_sorb(k,1) = iand(psi_ref(k,1,i), nonactive_sorb(k,1)) + nonactive_sorb(k,2) = iand(psi_ref(k,2,i), nonactive_sorb(k,2)) + end do + end do + do k=1, N_int + active_sorb(k,1) = iand(active_sorb(k,1), not(nonactive_sorb(k,1))) + active_sorb(k,2) = iand(active_sorb(k,2), not(nonactive_sorb(k,2))) + end do + end if + + + do i=1, N_det_non_ref + do k=1, N_int + det_noactive(k,1,i) = iand(psi_non_ref(k,1,i), not(active_sorb(k,1))) + det_noactive(k,2,i) = iand(psi_non_ref(k,2,i), not(active_sorb(k,2))) + end do + end do + + call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int) + + do i=1,N_det_non_ref + det_cepa0(:,:,i) = psi_non_ref(:,:,det_cepa0_idx(i)) + end do + + cepa0_shortcut(0) = 1 + cepa0_shortcut(1) = 1 + do i=2,N_det_non_ref + if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then + cepa0_shortcut(0) += 1 + cepa0_shortcut(cepa0_shortcut(0)) = i + end if + end do + cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1 + + if(.true.) then + do i=1,cepa0_shortcut(0) + n = cepa0_shortcut(i+1) - cepa0_shortcut(i) + call sort_dets_ab(det_cepa0(1,1,cepa0_shortcut(i)), idx, shortcut, n, N_int) + do k=1,n + idx(k) = det_cepa0_idx(cepa0_shortcut(i)-1+idx(k)) + end do + det_cepa0_idx(cepa0_shortcut(i):cepa0_shortcut(i)+n-1) = idx(:n) + end do + end if + + + do i=1,N_det_ref + do k=1, N_int + det_ref_active(k,1,i) = iand(psi_ref(k,1,i), active_sorb(k,1)) + det_ref_active(k,2,i) = iand(psi_ref(k,2,i), active_sorb(k,2)) + end do + end do + + do i=1,N_det_non_ref + do k=1, N_int + det_cepa0_active(k,1,i) = iand(psi_non_ref(k,1,det_cepa0_idx(i)), active_sorb(k,1)) + det_cepa0_active(k,2,i) = iand(psi_non_ref(k,2,det_cepa0_idx(i)), active_sorb(k,2)) + end do + end do + + do i=1,N_det_non_ref + if(.not. detEq(psi_non_ref(1,1,det_cepa0_idx(i)), det_cepa0(1,1,i),N_int)) stop "STOOOP" + end do + + searchance = 0d0 + child_num = 0 + do J = 1, N_det_ref + nlink(J) = 0 + do blok=1,cepa0_shortcut(0) + do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int) + if(degree <= 2) then + nlink(J) += 1 + linked(nlink(J),J) = k + child_num(k, J) = nlink(J) + blokMwen(nlink(J),J) = blok + searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok))) + end if + end do + end do + end do + print *, "pre done" +END_PROVIDER + + +! BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] +! use bitmasks +! implicit none +! integer :: i,j,k +! double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall +! integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref) +! +! ! provide lambda_mrcc +! npres = 0 +! delta_cas = 0d0 +! call wall_time(wall) +! print *, "dcas ", wall +! do i_state = 1, N_states +! !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref) +! do k=1,N_det_non_ref +! if(lambda_mrcc(i_state, k) == 0d0) cycle +! npre = 0 +! do i=1,N_det_ref +! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) +! if(Hki /= 0d0) then +! !!$OMP ATOMIC +! npres(i) += 1 +! npre += 1 +! ipre(npre) = i +! pre(npre) = Hki +! end if +! end do +! +! +! do i=1,npre +! do j=1,i +! !!$OMP ATOMIC +! delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k) +! end do +! end do +! end do +! !!$OMP END PARALLEL DO +! npre=0 +! do i=1,N_det_ref +! npre += npres(i) +! end do +! !stop +! do i=1,N_det_ref +! do j=1,i +! delta_cas(j,i,i_state) = delta_cas(i,j,i_state) +! end do +! end do +! end do +! +! call wall_time(wall) +! print *, "dcas", wall +! ! stop +! END_PROVIDER + + + BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] + use bitmasks + implicit none + integer :: i,j,k + double precision :: Hjk, Hki, Hij + !double precision, external :: get_dij + integer i_state, degree + + provide lambda_mrcc dIj + do i_state = 1, N_states + !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref,dij) + do i=1,N_det_ref + do j=1,i + call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int) + delta_cas(i,j,i_state) = 0d0 + do k=1,N_det_non_ref + + call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) + call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) + + delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) + !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int) + end do + delta_cas(j,i,i_state) = delta_cas(i,j,i_state) + end do + end do + !$OMP END PARALLEL DO + end do + END_PROVIDER + + + + +logical function isInCassd(a,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2) + integer(bit_kind) :: inac, virt + integer :: ni, i, deg + + + isInCassd = .false. + + deg = 0 + do i=1,2 + do ni=1,Nint + virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i))) + deg += popcnt(iand(virt, a(ni,i))) + if(deg > 2) return + end do + end do + + deg = 0 + do i=1,2 + do ni=1,Nint + inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i))) + deg += popcnt(xor(iand(inac,a(ni,i)), inac)) + if(deg > 2) return + end do + end do + isInCassd = .true. +end function + + +subroutine getHP(a,h,p,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2) + integer, intent(out) :: h, p + integer(bit_kind) :: inac, virt + integer :: ni, i, deg + + + !isInCassd = .false. + h = 0 + p = 0 + + deg = 0 + lp : do i=1,2 + do ni=1,Nint + virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i))) + deg += popcnt(iand(virt, a(ni,i))) + if(deg > 2) exit lp + end do + end do lp + p = deg + + deg = 0 + lh : do i=1,2 + do ni=1,Nint + inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i))) + deg += popcnt(xor(iand(inac,a(ni,i)), inac)) + if(deg > 2) exit lh + end do + end do lh + h = deg + !isInCassd = .true. +end function + + + BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref) + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) + double precision :: contrib, HIIi, HJk, wall + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2) + integer(bit_kind),allocatable :: sortRef(:,:,:) + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit, searchDet + logical, external :: is_in_wavefunction, detEq + !double precision, external :: get_dij + integer :: II, blok + integer*8, save :: notf = 0 + + call wall_time(wall) + allocate(idx_sorted_bit(N_det), sortRef(N_int,2,N_det_ref)) + + sortRef(:,:,:) = det_ref_active(:,:,:) + call sort_det(sortRef, sortRefIdx, N_det_ref, N_int) + + idx_sorted_bit(:) = -1 + do i=1,N_det_non_ref + idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i + enddo + + ! To provide everything + contrib = dij(1, 1, 1) + + do i_state = 1, N_states + delta_mrcepa0_ii(:,:) = 0d0 + delta_mrcepa0_ij(:,:,:) = 0d0 + + !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii) & + !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib) & + !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & + !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) & + !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) + do blok=1,cepa0_shortcut(0) + do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do II=1,N_det_ref + call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int) + if (degree > 2 ) cycle + + do ni=1,N_int + made_hole(ni,1) = iand(det_ref_active(ni,1,II), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II))) + made_hole(ni,2) = iand(det_ref_active(ni,2,II), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II))) + + made_particle(ni,1) = iand(det_cepa0_active(ni,1,i), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II))) + made_particle(ni,2) = iand(det_cepa0_active(ni,2,i), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II))) + end do + + + kloop: do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 !i + !if(lambda_mrcc(i_state, det_cepa0_idx(k)) == 0d0) cycle + + do ni=1,N_int + if(iand(made_hole(ni,1), det_cepa0_active(ni,1,k)) /= 0) cycle kloop + if(iand(made_particle(ni,1), det_cepa0_active(ni,1,k)) /= made_particle(ni,1)) cycle kloop + if(iand(made_hole(ni,2), det_cepa0_active(ni,2,k)) /= 0) cycle kloop + if(iand(made_particle(ni,2), det_cepa0_active(ni,2,k)) /= made_particle(ni,2)) cycle kloop + end do + do ni=1,N_int + myActive(ni,1) = xor(det_cepa0_active(ni,1,k), made_hole(ni,1)) + myActive(ni,1) = xor(myActive(ni,1), made_particle(ni,1)) + myActive(ni,2) = xor(det_cepa0_active(ni,2,k), made_hole(ni,2)) + myActive(ni,2) = xor(myActive(ni,2), made_particle(ni,2)) + end do + + j = searchDet(sortRef, myActive, N_det_ref, N_int) + if(j == -1) then + cycle + end if + j = sortRefIdx(j) + !$OMP ATOMIC + notf = notf+1 + + call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) + !contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k)) + contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + !$OMP ATOMIC + delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib + + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + !$OMP ATOMIC + delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) + end if + + end do kloop + end do + end do + end do + !$OMP END PARALLEL DO + end do + deallocate(idx_sorted_bit) + call wall_time(wall) + print *, "cepa0", wall, notf + !stop +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ + logical :: ok + double precision :: phase_Ji, phase_Ik, phase_Ii + double precision :: contrib, delta_IJk, HJk, HIk, HIl + integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2) + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + + integer :: II, blok + + provide delta_cas lambda_mrcc + allocate(idx_sorted_bit(N_det)) + idx_sorted_bit(:) = -1 + do i=1,N_det_non_ref + idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i + enddo + + do i_state = 1, N_states + delta_sub_ij(:,:,:) = 0d0 + delta_sub_ii(:,:) = 0d0 + + provide mo_bielec_integrals_in_map + + + !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) & + !$OMP private(i, J, k, degree, degree2, l, deg, ni) & + !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) & + !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & + !$OMP private(det_tmp, det_tmp2, II, blok) & + !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) & + !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb) + do i=1,N_det_non_ref + if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref + do J=1,N_det_ref + call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int) + if(degree == -1) cycle + + + do II=1,N_det_ref + call apply_excitation(psi_ref(1,1,II),exc_Ji,det_tmp,ok,N_int) + + if(.not. ok) cycle + l = get_index_in_psi_det_sorted_bit(det_tmp, N_int) + if(l == 0) cycle + l = idx_sorted_bit(l) + + call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl) + + do k=1,N_det_non_ref + if(lambda_mrcc(i_state, k) == 0d0) cycle + call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int) + + det_tmp(:,:) = 0_bit_kind + det_tmp2(:,:) = 0_bit_kind + + ok = .true. + do ni=1,N_int + det_tmp(ni,1) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,k)), not(active_sorb(ni,1))) + det_tmp(ni,2) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,i)), not(active_sorb(ni,1))) + ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2)))) + + det_tmp(ni,1) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,k)), not(active_sorb(ni,2))) + det_tmp(ni,2) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,i)), not(active_sorb(ni,2))) + ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2)))) + end do + + if(ok) cycle + + + call i_h_j(psi_ref(1,1,J), psi_non_ref(1,1,k), N_int, HJk) + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,k), N_int, HIk) + if(HJk == 0) cycle + !assert HIk == 0 + delta_IJk = HJk * HIk * lambda_mrcc(i_state, k) + call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) + if(ok) cycle + contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) + !$OMP ATOMIC + delta_sub_ij(II, i, i_state) += contrib + if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then + !$OMP ATOMIC + delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) + endif + end do + end do + end do + end do + !$OMP END PARALLEL DO + end do + deallocate(idx_sorted_bit) +END_PROVIDER + + +subroutine set_det_bit(det, p, s) + implicit none + integer(bit_kind),intent(inout) :: det(N_int, 2) + integer, intent(in) :: p, s + integer :: ni, pos + + ni = (p-1)/bit_kind_size + 1 + pos = mod(p-1, bit_kind_size) + det(ni,s) = ibset(det(ni,s), pos) +end subroutine + + +BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ] + implicit none + integer :: i,j + do i=1,N_det_ref + do j=1,N_det_non_ref + call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,j), N_int, h_(i,j)) + end do + end do +END_PROVIDER + + + +subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer,allocatable :: degree(:) + integer,allocatable :: idx(:) + logical :: good + + integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators) + integer,intent(in) :: N_miniList + + allocate(degree(psi_det_size)) + allocate(idx(0:psi_det_size)) + N_tq = 0 + + i_loop : do i=1,N_selected + do k=1, N_minilist + if(is_generable(miniList(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + + ! Select determinants that are triple or quadruple excitations + ! from the ref + good = .True. + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo i_loop +end + + +subroutine filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer,allocatable :: degree(:) + integer,allocatable :: idx(:) + logical :: good + + integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + integer(bit_kind),intent(in) :: microlist(Nint,2,*) + integer,intent(in) :: ptr_microlist(0:*) + integer,intent(in) :: N_microlist(0:*) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + + integer :: mobiles(2), smallerlist + + + allocate(degree(psi_det_size)) + allocate(idx(0:psi_det_size)) + N_tq = 0 + + i_loop : do i=1,N_selected + call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint) + if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then + smallerlist = mobiles(1) + else + smallerlist = mobiles(2) + end if + + if(N_microlist(smallerlist) > 0) then + do k=ptr_microlist(smallerlist), ptr_microlist(smallerlist)+N_microlist(smallerlist)-1 + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + end if + + if(N_microlist(0) > 0) then + do k=1, N_microlist(0) + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + end if + + ! Select determinants that are triple or quadruple excitations + ! from the ref + good = .True. + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo i_loop +end + + + + diff --git a/plugins/mrcepa0/dressing_slave.irp.f b/plugins/mrcepa0/dressing_slave.irp.f new file mode 100644 index 00000000..7c9d7fe0 --- /dev/null +++ b/plugins/mrcepa0/dressing_slave.irp.f @@ -0,0 +1,593 @@ +subroutine mrsc2_dressing_slave_tcp(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + call mrsc2_dressing_slave(0,i) +end + + +subroutine mrsc2_dressing_slave_inproc(i) + implicit none + integer, intent(in) :: i + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + call mrsc2_dressing_slave(1,i) +end + +subroutine mrsc2_dressing_slave(thread,iproc) + use f77_zmq + + implicit none + BEGIN_DOC +! Task for parallel MR-SC2 + END_DOC + integer, intent(in) :: thread, iproc +! integer :: j,l + integer :: rc + + integer :: worker_id, task_id + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + double precision, allocatable :: delta(:,:,:) + + + + integer :: i_state, i, i_I, J, k, k2, k1, kk, ll, degree, degree2, m, l, deg, ni, m2 + integer :: n(2) + integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, kn + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al + double precision :: diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv(N_states), cj_inv(N_states) + double precision :: contrib, wall, iwall + double precision, allocatable :: dleat(:,:,:) + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt + integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp + logical, external :: is_in_wavefunction, isInCassd, detEq + integer,allocatable :: komon(:) + logical :: komoned + !double precision, external :: get_dij + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_push = new_zmq_push_socket(thread) + + call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) + + allocate (dleat(N_states, N_det_non_ref, 2), delta(N_states,0:N_det_non_ref, 2)) + allocate(komon(0:N_det_non_ref)) + + do + call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) + if (task_id == 0) exit + read (task,*) i_I, J, k1, k2 + do i_state=1, N_states + ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state) + cj_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state) + end do + !delta = 0.d0 + n = 0 + delta(:,0,:) = 0d0 + delta(:,:nlink(J),1) = 0d0 + delta(:,:nlink(i_I),2) = 0d0 + komon(0) = 0 + komoned = .false. + + + + + do kk = k1, k2 + k = det_cepa0_idx(linked(kk, i_I)) + blok = blokMwen(kk, i_I) + + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int) + + if(J /= i_I) then + call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int) + if(.not. ok) cycle + + l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int) + if(l == -1) cycle + ll = cepa0_shortcut(blok)-1+l + l = det_cepa0_idx(ll) + ll = child_num(ll, J) + else + l = k + ll = kk + end if + + + if(.not. komoned) then + m = 0 + m2 = 0 + + do while(m < nlink(i_I) .and. m2 < nlink(J)) + m += 1 + m2 += 1 + if(linked(m, i_I) < linked(m2, J)) then + m2 -= 1 + cycle + else if(linked(m, i_I) > linked(m2, J)) then + m -= 1 + cycle + end if + i = det_cepa0_idx(linked(m, i_I)) + + if(h_(J,i) == 0.d0) cycle + if(h_(i_I,i) == 0.d0) cycle + + !ok = .false. + !do i_state=1, N_states + ! if(lambda_mrcc(i_state, i) /= 0d0) then + ! ok = .true. + ! exit + ! end if + !end do + !if(.not. ok) cycle +! + + komon(0) += 1 + kn = komon(0) + komon(kn) = i + + +! call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int) +! if(I_i /= J) call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ii,N_int) +! if(I_i == J) phase_Ii = phase_Ji + + do i_state = 1,N_states + dkI = h_(J,i) * dij(i_I, i, i_state)!get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int) + !dkI = h_(J,i) * h_(i_I,i) * lambda_mrcc(i_state, i) + dleat(i_state, kn, 1) = dkI + dleat(i_state, kn, 2) = dkI + end do + + end do + + komoned = .true. + end if + + + do m = 1, komon(0) + + i = komon(m) + + call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) + if(.not. ok) cycle + if(HP(1,i) + HP(1,k) <= 2 .and. HP(2,i) + HP(2,k) <= 2) then +! if(is_in_wavefunction(det_tmp, N_int)) cycle + cycle + end if + + !if(isInCassd(det_tmp, N_int)) cycle + + do i_state = 1, N_states + !if(lambda_mrcc(i_state, i) == 0d0) cycle + + + !contrib = h_(i_I,k) * lambda_mrcc(i_state, k) * dleat(i_state, m, 2)! * phase_al + contrib = dij(i_I, k, i_state) * dleat(i_state, m, 2) + delta(i_state,ll,1) += contrib + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then + delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state) + endif + + if(I_i == J) cycle + !contrib = h_(J,l) * lambda_mrcc(i_state, l) * dleat(i_state, m, 1)! * phase_al + contrib = dij(J, l, i_state) * dleat(i_state, m, 1) + delta(i_state,kk,2) += contrib + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state) + end if + enddo !i_state + end do ! while + end do ! kk + + + call push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) + +! end if + + enddo + + deallocate(delta) + + call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + +end + + +subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) + use f77_zmq + implicit none + BEGIN_DOC +! Push integrals in the push socket + END_DOC + + integer, intent(in) :: i_I, J + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + double precision,intent(inout) :: delta(N_states, 0:N_det_non_ref, 2) + integer, intent(in) :: task_id + integer :: rc , i_state, i, kk, li + integer,allocatable :: idx(:,:) + integer ::n(2) + logical :: ok + + allocate(idx(N_det_non_ref,2)) + rc = f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + + do kk=1,2 + n(kk)=0 + if(kk == 1) li = nlink(j) + if(kk == 2) li = nlink(i_I) + do i=1, li + ok = .false. + do i_state=1,N_states + if(delta(i_state, i, kk) /= 0d0) then + ok = .true. + exit + end if + end do + + if(ok) then + n(kk) += 1 +! idx(n,kk) = i + if(kk == 1) then + idx(n(1),1) = det_cepa0_idx(linked(i, J)) + else + idx(n(2),2) = det_cepa0_idx(linked(i, i_I)) + end if + + do i_state=1, N_states + delta(i_state, n(kk), kk) = delta(i_state, i, kk) + end do + end if + end do + + rc = f77_zmq_send( zmq_socket_push, n(kk), 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, n, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + if(n(kk) /= 0) then + rc = f77_zmq_send( zmq_socket_push, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) ! delta(1,0,1) = delta_I delta(1,0,2) = delta_J + if (rc /= (n(kk)+1)*8*N_states) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_send( zmq_socket_push, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) + if (rc /= n(kk)*4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, delta, 8*n(kk), ZMQ_SNDMORE)' + stop 'error' + endif + end if + end do + + + rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, task_id, 4, 0)' + stop 'error' + endif + +! ! Activate is zmq_socket_push is a REQ +! integer :: idummy +! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) +! if (rc /= 4) then +! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)' +! stop 'error' +! endif +end + + + +subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) + use f77_zmq + implicit none + BEGIN_DOC +! Push integrals in the push socket + END_DOC + + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + integer, intent(out) :: i_I, J, n(2) + double precision, intent(inout) :: delta(N_states, 0:N_det_non_ref, 2) + integer, intent(out) :: task_id + integer :: rc , i, kk + integer,intent(inout) :: idx(N_det_non_ref,2) + logical :: ok + + rc = f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + do kk = 1, 2 + rc = f77_zmq_recv( zmq_socket_pull, n(kk), 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, n, 4, ZMQ_SNDMORE)' + stop 'error' + endif + + if(n(kk) /= 0) then + rc = f77_zmq_recv( zmq_socket_pull, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) + if (rc /= (n(kk)+1)*8*N_states) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)' + stop 'error' + endif + + rc = f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) + if (rc /= n(kk)*4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, delta, n(kk)*4, ZMQ_SNDMORE)' + stop 'error' + endif + end if + end do + + rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)' + stop 'error' + endif + + +! ! Activate is zmq_socket_pull is a REP +! integer :: idummy +! rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0) +! if (rc /= 4) then +! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)' +! stop 'error' +! endif +end + + + +subroutine mrsc2_dressing_collector(delta_ii_,delta_ij_) + use f77_zmq + implicit none + BEGIN_DOC +! Collects results from the AO integral calculation + END_DOC + + double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) + double precision,intent(inout) :: delta_ii_(N_states,N_det_ref) + +! integer :: j,l + integer :: rc + + double precision, allocatable :: delta(:,:,:) + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer*8 :: control, accu + integer :: task_id, more + + integer :: I_i, J, l, i_state, n(2), kk + integer,allocatable :: idx(:,:) + + delta_ii_(:,:) = 0d0 + delta_ij_(:,:,:) = 0d0 + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + + allocate ( delta(N_states,0:N_det_non_ref,2) ) + + allocate(idx(N_det_non_ref,2)) + more = 1 + do while (more == 1) + + call pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) + + + do l=1, n(1) + do i_state=1,N_states + delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1) + end do + end do + + do l=1, n(2) + do i_state=1,N_states + delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2) + end do + end do + + +! +! do l=1,nlink(J) +! do i_state=1,N_states +! delta_ij_(i_state,det_cepa0_idx(linked(l,J)),i_I) += delta(i_state,l,1) +! delta_ij_(i_state,det_cepa0_idx(linked(l,i_I)),j) += delta(i_state,l,2) +! end do +! end do +! + if(n(1) /= 0) then + do i_state=1,N_states + delta_ii_(i_state,i_I) += delta(i_state,0,1) + end do + end if + + if(n(2) /= 0) then + do i_state=1,N_states + delta_ii_(i_state,J) += delta(i_state,0,2) + end do + end if + + + if (task_id /= 0) then + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) + endif + + + enddo + deallocate( delta ) + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + +end + + + + + BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ] + implicit none + + integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2 + integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, nex, nzer, ntot +! integer, allocatable :: linked(:,:), blokMwen(:, :), nlink(:) + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states) + double precision :: contrib, wall, iwall ! , searchance(N_det_ref) + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt + integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp + logical, external :: is_in_wavefunction, isInCassd, detEq + character*(512) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer :: KKsize = 1000000 + + + call new_parallel_job(zmq_to_qp_run_socket,'mrsc2') + + + call wall_time(iwall) +! allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref)) + + +! searchance = 0d0 +! do J = 1, N_det_ref +! nlink(J) = 0 +! do blok=1,cepa0_shortcut(0) +! do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 +! call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int) +! if(degree <= 2) then +! nlink(J) += 1 +! linked(nlink(J),J) = k +! blokMwen(nlink(J),J) = blok +! searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok))) +! end if +! end do +! end do +! end do + + + +! stop +nzer = 0 +ntot = 0 + do nex = 3, 0, -1 + print *, "los ",nex + do I_s = N_det_ref, 1, -1 +! if(mod(I_s,1) == 0) then +! call wall_time(wall) +! wall = wall-iwall +! print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall +! end if + + + do J_s = 1, I_s + + call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int) + if(degree /= nex) cycle + if(nex == 3) nzer = nzer + 1 + ntot += 1 +! if(degree > 3) then +! deg += 1 +! cycle +! else if(degree == -10) then +! KKsize = 100000 +! else +! KKsize = 1000000 +! end if + + + + if(searchance(I_s) < searchance(J_s)) then + i_I = I_s + J = J_s + else + i_I = J_s + J = I_s + end if + + KKsize = nlink(1) + if(nex == 0) KKsize = int(float(nlink(1)) / float(nlink(i_I)) * (float(nlink(1)) / 64d0)) + + !if(KKsize == 0) stop "ZZEO" + + do kk = 1 , nlink(i_I), KKsize + write(task,*) I_i, J, kk, int(min(kk+KKsize-1, nlink(i_I))) + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + end do + + ! do kk = 1 , nlink(i_I) + ! k = linked(kk,i_I) + ! blok = blokMwen(kk,i_I) + ! write(task,*) I_i, J, k, blok + ! call add_task_to_taskserver(zmq_to_qp_run_socket,task) + ! + ! enddo !kk + enddo !J + + enddo !I + end do ! nex + print *, "tasked" +! integer(ZMQ_PTR) ∷ collector_thread +! external ∷ ao_bielec_integrals_in_map_collector +! rc = pthread_create(collector_thread, mrsc2_dressing_collector) + print *, nzer, ntot, float(nzer) / float(ntot) + provide nproc + !$OMP PARALLEL DEFAULT(none) SHARED(delta_ii_old,delta_ij_old) PRIVATE(i) NUM_THREADS(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call mrsc2_dressing_collector(delta_ii_old,delta_ij_old) + else + call mrsc2_dressing_slave_inproc(i) + endif + !$OMP END PARALLEL + +! rc = pthread_join(collector_thread) + call end_parallel_job(zmq_to_qp_run_socket, 'mrsc2') + + +END_PROVIDER + + + diff --git a/plugins/mrcepa0/mrcc.irp.f b/plugins/mrcepa0/mrcc.irp.f new file mode 100644 index 00000000..91592e62 --- /dev/null +++ b/plugins/mrcepa0/mrcc.irp.f @@ -0,0 +1,19 @@ +program mrsc2sub + implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + mrmode = 3 + + read_wf = .True. + SOFT_TOUCH read_wf + call print_cas_coefs + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) +end + diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f new file mode 100644 index 00000000..34d3dec5 --- /dev/null +++ b/plugins/mrcepa0/mrcepa0.irp.f @@ -0,0 +1,19 @@ +program mrcepa0 + implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + mrmode = 1 + + read_wf = .True. + SOFT_TOUCH read_wf + call print_cas_coefs + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) +end + diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f new file mode 100644 index 00000000..53a0822d --- /dev/null +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -0,0 +1,169 @@ + + +subroutine run(N_st,energy) + implicit none + + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration + double precision :: E_past(4), lambda + + integer :: n_it_mrcc_max + double precision :: thresh_mrcc + + + + thresh_mrcc = 1d-7 + n_it_mrcc_max = 10 + + if(n_it_mrcc_max == 1) then + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef ci_energy_dressed + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + energy(:) = ci_energy_dressed(:) + else + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + lambda = 1.d0 + do while (delta_E > thresh_mrcc) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA0 Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy") + call diagonalize_ci_dressed(lambda) + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + call save_wavefunction + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + if (iteration > n_it_mrcc_max) then + exit + endif + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") + energy(:) = ci_energy_dressed(:) + endif +end + + +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + double precision :: pt3(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + pt3 = 0d0 + !if(lambda_mrcc_pt2(0) == 0) return + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + + N_det_generators = lambda_mrcc_pt3(0) + N_det_ref + N_det_selectors = lambda_mrcc_pt3(0) + N_det_ref + + psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + + do i=N_det_ref+1,N_det_generators + j = lambda_mrcc_pt3(i-N_det_ref) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + psi_selectors(k,1,i) = psi_non_ref(k,1,j) + psi_selectors(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + psi_selectors_coef(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + + SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized! psi_coef_energy_diagonalized + call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st) + + N_det_generators = N_det_non_ref + N_det_ref + N_det_selectors = N_det_non_ref + N_det_ref + + psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref) + psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:) + + do i=N_det_ref+1,N_det_generators + j = i-N_det_ref + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + psi_selectors(k,1,i) = psi_non_ref(k,1,j) + psi_selectors(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + psi_selectors_coef(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + + SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized! psi_coef_energy_diagonalized + call H_apply_mrcepa_PT2(pt3, norm_pert, H_pert_diag, N_st) + + +!!!!!!!!!!!!!!!! + + + + print *, "2-3 :",pt2, pt3 + print *, lambda_mrcc_pt3(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1) + pt2 = pt2 - pt3 + + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + + + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) + +end + + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + call write_double(6,ci_energy(1),"Initial CI energy") + +end + + diff --git a/plugins/mrcepa0/mrsc2.irp.f b/plugins/mrcepa0/mrsc2.irp.f new file mode 100644 index 00000000..d0f44a33 --- /dev/null +++ b/plugins/mrcepa0/mrsc2.irp.f @@ -0,0 +1,19 @@ +program mrsc2 + implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + mrmode = 2 + read_wf = .True. + SOFT_TOUCH read_wf + call print_cas_coefs + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) +end + + diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 5d9198c4..1f03908b 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -11,7 +11,7 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) PROVIDE N_int PROVIDE N_det - + $declarations @@ -180,7 +180,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl $initialization - + $omp_parallel !$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 7a54bdbc..2f53c799 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -165,7 +165,7 @@ logical function is_connected_to(key,keys,Nint,Ndet) integer :: i, l integer :: degree_x2 - + logical, external :: is_generable_cassd ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -183,12 +183,35 @@ logical function is_connected_to(key,keys,Nint,Ndet) if (degree_x2 > 4) then cycle else +! if(.not. is_generable_cassd(keys(1,1,i), key(1,1), Nint)) cycle !!!Nint==1 !!!!! is_connected_to = .true. return endif enddo end + +logical function is_generable_cassd(det1, det2, Nint) !!! TEST Cl HARD !!!!! + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) + integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t + double precision :: phase + + is_generable_cassd = .false. + call get_excitation(det1, det2, exc, degree, phase, Nint) + if(degree == -1) return + if(degree == 0) then + is_generable_cassd = .true. + return + end if + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(degree == 1 .and. h1 <= 11) is_generable_cassd = .true. + if(degree == 2 .and. h1 <= 11 .and. h2 <= 11) is_generable_cassd = .true. +end function + + logical function is_connected_to_by_mono(key,keys,Nint,Ndet) use bitmasks implicit none diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index 7c1f43c2..a4166e10 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -241,8 +241,8 @@ subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) END_DOC integer, intent(in) :: Nint, N_key integer(bit_kind),intent(inout) :: key(Nint,2,N_key) - integer,intent(out) :: idx(N_key) - integer,intent(out) :: shortcut(0:N_key+1) + integer,intent(inout) :: idx(N_key) + integer,intent(inout) :: shortcut(0:N_key+1) integer(bit_kind) :: tmp(Nint, 2) integer :: tmpidx,i,ni diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index a7727cda..52d2cc53 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -664,3 +664,44 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde end +logical function detEq(a,b,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detEq = .false. + do i=1,2 + do ni=1,Nint + if(a(ni,i) /= b(ni,i)) return + end do + end do + detEq = .true. +end function + + +integer function detCmp(a,b,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detCmp = 0 + do i=1,2 + do ni=Nint,1,-1 + + if(a(ni,i) < b(ni,i)) then + detCmp = -1 + return + else if(a(ni,i) > b(ni,i)) then + detCmp = 1 + return + end if + + end do + end do +end function + + diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 98d4f71d..133d9e52 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -139,6 +139,72 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) end +subroutine decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) + use bitmasks + implicit none + BEGIN_DOC + ! Decodes the exc arrays returned by get_excitation. + ! h1,h2 : Holes + ! p1,p2 : Particles + ! s1,s2 : Spins (1:alpha, 2:beta) + ! degree : Degree of excitation + END_DOC + integer, intent(in) :: exc(0:2,2,2),degree + integer*2, intent(out) :: h1,h2,p1,p2,s1,s2 + ASSERT (degree > 0) + ASSERT (degree < 3) + + select case(degree) + case(2) + if (exc(0,1,1) == 2) then + h1 = exc(1,1,1) + h2 = exc(2,1,1) + p1 = exc(1,2,1) + p2 = exc(2,2,1) + s1 = 1 + s2 = 1 + else if (exc(0,1,2) == 2) then + h1 = exc(1,1,2) + h2 = exc(2,1,2) + p1 = exc(1,2,2) + p2 = exc(2,2,2) + s1 = 2 + s2 = 2 + else + h1 = exc(1,1,1) + h2 = exc(1,1,2) + p1 = exc(1,2,1) + p2 = exc(1,2,2) + s1 = 1 + s2 = 2 + endif + case(1) + if (exc(0,1,1) == 1) then + h1 = exc(1,1,1) + h2 = 0 + p1 = exc(1,2,1) + p2 = 0 + s1 = 1 + s2 = 0 + else + h1 = exc(1,1,2) + h2 = 0 + p1 = exc(1,2,2) + p2 = 0 + s1 = 2 + s2 = 0 + endif + case(0) + h1 = 0 + p1 = 0 + h2 = 0 + p2 = 0 + s1 = 0 + s2 = 0 + end select +end + + subroutine get_double_excitation(det1,det2,exc,phase,Nint) use bitmasks implicit none @@ -915,7 +981,6 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis fullMatch = .false. N_miniList = 0 N_subList = 0 - l = popcnt(key_mask(1,1)) + popcnt(key_mask(1,2)) do ni = 2,Nint l = l + popcnt(key_mask(ni,1)) + popcnt(key_mask(ni,2)) @@ -948,8 +1013,13 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis miniList(ni,2,N_minilist) = fullList(ni,2,i) enddo else if(k == 0) then - fullMatch = .true. - return + N_minilist += 1 + do ni=1,Nint + miniList(ni,1,N_minilist) = fullList(ni,1,i) + miniList(ni,2,N_minilist) = fullList(ni,2,i) + enddo +! fullMatch = .true. +! return end if end do end if @@ -1761,4 +1831,3 @@ subroutine apply_excitation(det, exc, res, ok, Nint) ok = .true. end subroutine - diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index 8d5726f5..2eec0dea 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -10,7 +10,7 @@ integer*8 function spin_det_search_key(det,Nint) use bitmasks implicit none BEGIN_DOC -! Return an integer*8 corresponding to a determinant index for searching +! Return an integer(8) corresponding to a determinant index for searching END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: det(Nint) @@ -64,9 +64,9 @@ BEGIN_TEMPLATE integer :: i,j,k integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8 :: last_key - integer*8, external :: spin_det_search_key + integer(8), allocatable :: bit_tmp(:) + integer(8) :: last_key + integer(8), external :: spin_det_search_key logical,allocatable :: duplicate(:) allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) ) @@ -149,8 +149,8 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) integer(bit_kind), intent(in) :: key(Nint) integer :: i, ibegin, iend, istep, l - integer*8 :: det_ref, det_search - integer*8, external :: spin_det_search_key + integer(8) :: det_ref, det_search + integer(8), external :: spin_det_search_key logical :: in_wavefunction in_wavefunction = .False. @@ -231,8 +231,8 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) integer(bit_kind), intent(in) :: key(Nint) integer :: i, ibegin, iend, istep, l - integer*8 :: det_ref, det_search - integer*8, external :: spin_det_search_key + integer(8) :: det_ref, det_search + integer(8), external :: spin_det_search_key logical :: in_wavefunction in_wavefunction = .False. @@ -305,10 +305,10 @@ end subroutine write_spindeterminants use bitmasks implicit none - integer*8, allocatable :: tmpdet(:,:) + integer(8), allocatable :: tmpdet(:,:) integer :: N_int2 integer :: i,j,k - integer*8 :: det_8(100) + integer(8) :: det_8(100) integer(bit_kind) :: det_bk((100*8)/bit_kind) equivalence (det_8, det_bk) diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 11de1270..00f61101 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -148,10 +148,10 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) integer, intent(in) :: LDA, ldc, n, m double precision, intent(in) :: overlap(lda,n) double precision, intent(inout) :: C(ldc,n) - double precision :: U(ldc,n) - double precision :: Vt(lda,n) - double precision :: D(n) - double precision :: S_half(lda,n) + double precision, allocatable :: U(:,:) + double precision, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + double precision, allocatable :: S_half(:,:) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j, k @@ -159,6 +159,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) return endif + allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n)) + call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) !$OMP PARALLEL DEFAULT(NONE) & @@ -203,6 +205,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) + deallocate(U,Vt,S_half,D) end