diff --git a/ocaml/create_git_sha1.sh b/ocaml/create_git_sha1.sh index 7b47e96f..f1fb7fa6 100755 --- a/ocaml/create_git_sha1.sh +++ b/ocaml/create_git_sha1.sh @@ -2,7 +2,7 @@ SHA1=$(git log -1 | head -1 | cut -d ' ' -f 2) DATE=$(git log -1 | grep Date | cut -d ':' -f 2-) -MESSAGE=$(git log -1 | tail -1) +MESSAGE=$(git log -1 | tail -1 | sed 's/"/\\"/g') cat << EOF > Git.ml open Core.Std let sha1 = "$SHA1" |> String.strip diff --git a/plugins/Alavi/.gitignore b/plugins/Alavi/.gitignore deleted file mode 100644 index e4e1a2ab..00000000 --- a/plugins/Alavi/.gitignore +++ /dev/null @@ -1,23 +0,0 @@ -# Automatically created by $QP_ROOT/scripts/module/module_handler.py -.ninja_deps -.ninja_log -AO_Basis -Bitmask -Determinants -Electrons -Ezfio_files -IRPF90_man -IRPF90_temp -Integrals_Bielec -Integrals_Monoelec -MO_Basis -Makefile -Makefile.depend -Nuclei -Pseudo -Utils -alavi_graph -ezfio_interface.irp.f -irpf90.make -irpf90_entities -tags \ No newline at end of file diff --git a/plugins/Alavi/NEEDED_CHILDREN_MODULES b/plugins/Alavi/NEEDED_CHILDREN_MODULES deleted file mode 100644 index aae89501..00000000 --- a/plugins/Alavi/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Determinants diff --git a/plugins/Alavi/README.rst b/plugins/Alavi/README.rst deleted file mode 100644 index f2194755..00000000 --- a/plugins/Alavi/README.rst +++ /dev/null @@ -1,23 +0,0 @@ -===== -alavi -===== - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -`alavi_graph `_ - Undocumented - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -.. image:: tree_dependency.png - -* `Determinants `_ - diff --git a/plugins/Alavi/alavi_graph.irp.f b/plugins/Alavi/alavi_graph.irp.f deleted file mode 100644 index 4e953add..00000000 --- a/plugins/Alavi/alavi_graph.irp.f +++ /dev/null @@ -1,28 +0,0 @@ -program alavi_graph - implicit none - integer :: exc(0:2,2,2),h1,p1,h2,p2,s1,s2 - double precision :: phase - - read_wf = .True. - touch read_wf - - integer :: k,degree - double precision :: hii - - do k=1,N_det - call get_excitation_degree(psi_det(1,1,1),psi_det(1,1,k),degree,N_int) - call i_H_j(psi_det(1,1,k),psi_det(1,1,k),N_int,hii) - print*, k,abs(psi_coef(k,1)), hii,degree - -! if (degree == 2) then -! call get_excitation(psi_det(1,1,1),psi_det(1,1,k),exc,degree,phase,N_int) -! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) -! print*, h1,h2,hii, abs(psi_coef(k,1)) -! endif -! - - - enddo -end - -!plot "test.dat" u (abs($2)):(abs($3)):4 w p palette \ No newline at end of file diff --git a/plugins/Alavi/tree_dependency.png b/plugins/Alavi/tree_dependency.png deleted file mode 100644 index b4f0df8b..00000000 Binary files a/plugins/Alavi/tree_dependency.png and /dev/null differ diff --git a/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f b/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f deleted file mode 100644 index 8adab518..00000000 --- a/plugins/CAS_SD_ZMQ/ezfio_interface.irp.f +++ /dev/null @@ -1,4 +0,0 @@ -! DO NOT MODIFY BY HAND -! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py -! from file /home/scemama/quantum_package/src/CAS_SD_ZMQ/EZFIO.cfg - diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f index a53064b4..0d816f3e 100644 --- a/plugins/Full_CI/full_ci.irp.f +++ b/plugins/Full_CI/full_ci.irp.f @@ -3,6 +3,11 @@ program full_ci integer :: i,k + print *, '====================================================================' + print *, 'This program is slow. Consider using the Full_CI_ZMQ module instead.' + print *, '====================================================================' + call sleep(2) + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) integer :: N_st, degree N_st = N_states diff --git a/plugins/MRPT_Utils/ezfio_interface.irp.f b/plugins/MRPT_Utils/ezfio_interface.irp.f deleted file mode 100644 index 6bd8931d..00000000 --- a/plugins/MRPT_Utils/ezfio_interface.irp.f +++ /dev/null @@ -1,23 +0,0 @@ -! DO NOT MODIFY BY HAND -! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py -! from file /home/scemama/quantum_package/src/MRPT_Utils/EZFIO.cfg - - -BEGIN_PROVIDER [ logical, do_third_order_1h1p ] - implicit none - BEGIN_DOC -! If true, compute the third order contribution for the 1h1p - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrpt_utils_do_third_order_1h1p(has) - if (has) then - call ezfio_get_mrpt_utils_do_third_order_1h1p(do_third_order_1h1p) - else - print *, 'mrpt_utils/do_third_order_1h1p not found in EZFIO file' - stop 1 - endif - -END_PROVIDER diff --git a/plugins/Psiref_threshold/.gitignore b/plugins/Psiref_threshold/.gitignore deleted file mode 100644 index d98a4abc..00000000 --- a/plugins/Psiref_threshold/.gitignore +++ /dev/null @@ -1,29 +0,0 @@ -# Automatically created by $QP_ROOT/scripts/module/module_handler.py -.ninja_deps -.ninja_log -AO_Basis -Bitmask -Determinants -Electrons -Ezfio_files -Generators_full -Hartree_Fock -IRPF90_man -IRPF90_temp -Integrals_Bielec -Integrals_Monoelec -MOGuess -MO_Basis -Makefile -Makefile.depend -Nuclei -Perturbation -Properties -Pseudo -Selectors_full -Utils -ezfio_interface.irp.f -irpf90.make -irpf90_entities -mrcc_general -tags \ No newline at end of file diff --git a/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES b/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES deleted file mode 100644 index 7e790003..00000000 --- a/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Psiref_Utils diff --git a/plugins/Psiref_threshold/README.rst b/plugins/Psiref_threshold/README.rst deleted file mode 100644 index 3d4726e1..00000000 --- a/plugins/Psiref_threshold/README.rst +++ /dev/null @@ -1,24 +0,0 @@ -======================= -Psiref_threshold Module -======================= - - -Reference wave function is defined as all determinants with coefficients -such that | c_i/c_max | > threshold. - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. by the `update_README.py` script. - -.. image:: tree_dependency.png - -* `Psiref_Utils `_ - diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f deleted file mode 100644 index ee69ef5c..00000000 --- a/plugins/Psiref_threshold/psi_ref.irp.f +++ /dev/null @@ -1,41 +0,0 @@ -use bitmasks - - BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] -&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] -&BEGIN_PROVIDER [ integer, N_det_ref ] - implicit none - BEGIN_DOC - ! Reference wave function, defined as determinants with amplitudes > 0.05 - ! idx_ref gives the indice of the ref determinant in psi_det. - END_DOC - integer :: i, k, l - logical :: good - double precision, parameter :: threshold=0.05d0 - double precision :: t(N_states) - N_det_ref = 0 - do l = 1, N_states - t(l) = threshold * abs_psi_coef_max(l) - enddo - do i=1,N_det - good = .False. - do l=1, N_states - psi_ref_coef(i,l) = 0.d0 - good = good.or.(dabs(psi_coef(i,l)) > t(l)) - enddo - if (good) then - N_det_ref = N_det_ref+1 - do k=1,N_int - psi_ref(k,1,N_det_ref) = psi_det(k,1,i) - psi_ref(k,2,N_det_ref) = psi_det(k,2,i) - enddo - idx_ref(N_det_ref) = i - do k=1,N_states - psi_ref_coef(N_det_ref,k) = psi_coef(i,k) - enddo - endif - enddo - call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') - -END_PROVIDER - diff --git a/plugins/Psiref_threshold/tree_dependency.png b/plugins/Psiref_threshold/tree_dependency.png deleted file mode 100644 index 9c2088e1..00000000 Binary files a/plugins/Psiref_threshold/tree_dependency.png and /dev/null differ diff --git a/plugins/mrcc_selected/dressing.irp.f b/plugins/mrcc_selected/dressing.irp.f deleted file mode 100644 index c772e2aa..00000000 --- a/plugins/mrcc_selected/dressing.irp.f +++ /dev/null @@ -1,1076 +0,0 @@ -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) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_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 - delta_ij_s2_mrcc = 0d0 - delta_ii_s2_mrcc = 0d0 - PROVIDE dij - 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, delta_ii_s2_mrcc, delta_ij_s2_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, 1000) == 0) print *, "mrcc ", gen, "/", N_det_generators - do h=1, hh_shortcut(0) - call apply_hole_local(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_local(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) then - call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc, delta_ij_s2_mrcc, delta_ii_s2_mrcc, gen,n,buf,N_int,omask) - endif - - 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_,delta_ij_s2_, delta_ii_s2_,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) - double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref) - double precision, intent(inout) :: delta_ii_s2_(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, sla, dIk(N_states), dka(N_states), dIa(N_states) - double precision, allocatable :: dIa_hla(:,:), dIa_sla(:,:) - 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(:), sij_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), sij_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), dIa_sla(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)) - call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_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) - sla = sij_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 - dIa_sla(i_state,k_sd) = dIa(i_state) * sla - 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.1.d-3)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_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + dIa_sla(i_state,k_sd) - delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(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) + 0.5d0*dIa_hla(i_state,k_sd) - delta_ij_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + 0.5d0*dIa_sla(i_state,k_sd) - enddo - endif - enddo - call omp_unset_lock( psi_ref_lock(i_I) ) - enddo - enddo - deallocate (dIa_hla,dIa_sla,hij_cache,sij_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) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2, (N_states, N_det_ref) ] - use bitmasks - implicit none - integer :: i, j, i_state - - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc - - if(mrmode == 3) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_mrcc(i_state,i) - enddo - do j = 1, N_det_non_ref - do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc(i_state,j,i) - enddo - end do - end do - - ! =-=-= BEGIN STATE AVERAGE -! do i = 1, N_det_ref -! delta_ii(:,i)= delta_ii_mrcc(1,i) -! delta_ii_s2(:,i)= delta_ii_s2_mrcc(1,i) -! do i_state = 2, N_states -! delta_ii(:,i) += delta_ii_mrcc(i_state,i) -! delta_ii_s2(:,i) += delta_ii_s2_mrcc(i_state,i) -! enddo -! do j = 1, N_det_non_ref -! delta_ij(:,j,i) = delta_ij_mrcc(1,j,i) -! delta_ij_s2(:,j,i) = delta_ij_s2_mrcc(1,j,i) -! do i_state = 2, N_states -! delta_ij(:,j,i) += delta_ij_mrcc(i_state,j,i) -! delta_ij_s2(:,j,i) += delta_ij_s2_mrcc(i_state,j,i) -! enddo -! end do -! end do -! delta_ij = delta_ij * (1.d0/dble(N_states)) -! delta_ii = delta_ii * (1.d0/dble(N_states)) - ! =-=-= END STATE AVERAGE - ! - ! 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 - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_old(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_old(i_state,i) - enddo - do j = 1, N_det_non_ref - do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_old(i_state,j,i) - enddo - end do - end do - else if(mrmode == 1) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_ii_s2(i_state,i)= delta_mrcepa0_ii_s2(i,i_state) - enddo - do j = 1, N_det_non_ref - do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_ij_s2(i_state,j,i) = delta_mrcepa0_ij_s2(i,j,i_state) - enddo - end do - end do - else - stop "invalid mrmode" - end if -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) ] -&BEGIN_PROVIDER [ double precision, delta_cas_s2, (N_det_ref, N_det_ref, N_states) ] - use bitmasks - implicit none - integer :: i,j,k - double precision :: Sjk,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,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,delta_cas_s2,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 - delta_cas_s2(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 get_s2(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Sjk) - - delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) - delta_cas_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k) - end do - delta_cas(j,i,i_state) = delta_cas(i,j,i_state) - delta_cas_s2(j,i,i_state) = delta_cas_s2(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) ] -&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_s2, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_s2, (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, contrib2, contrib_s2, contrib2_s2, 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) - - delta_mrcepa0_ii(:,:) = 0d0 - delta_mrcepa0_ij(:,:,:) = 0d0 - delta_mrcepa0_ii_s2(:,:) = 0d0 - delta_mrcepa0_ij_s2(:,:,:) = 0d0 - - do i_state = 1, N_states - !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) & - !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) & - !$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, delta_cas_s2) & - !$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) * dij(J, det_cepa0_idx(k), i_state) - contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) - - if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then - contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) - contrib2_s2 = contrib_s2 / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) - !$OMP ATOMIC - delta_mrcepa0_ii(J,i_state) -= contrib2 - delta_mrcepa0_ii_s2(J,i_state) -= contrib2_s2 - else - contrib = contrib * 0.5d0 - contrib_s2 = contrib_s2 * 0.5d0 - end if - !$OMP ATOMIC - delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib - delta_mrcepa0_ij_s2(J, det_cepa0_idx(i), i_state) += contrib_s2 - - 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 - -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, contrib2, 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, contrib2, 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) - if(dabs(psi_ref_coef(II,i_state)).ge.1.d-3) then - contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) - !$OMP ATOMIC - delta_sub_ii(II,i_state) -= contrib2 - else - contrib = contrib * 0.5d0 - endif - !$OMP ATOMIC - delta_sub_ij(II, i, i_state) += contrib - 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_cache, (N_det_ref,N_det_non_ref) ] -&BEGIN_PROVIDER [ double precision, s2_cache, (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_cache(i,j)) - call get_s2(psi_ref(1,1,i), psi_non_ref(1,1,j), N_int, s2_cache(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/mrcc_selected/dressing_slave.irp.f b/plugins/mrcc_selected/dressing_slave.irp.f deleted file mode 100644 index c2e5dd55..00000000 --- a/plugins/mrcc_selected/dressing_slave.irp.f +++ /dev/null @@ -1,601 +0,0 @@ -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(:,:,:), delta_s2(:,:,:) - - - - 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, contrib_s2, wall, iwall - double precision, allocatable :: dleat(:,:,:), dleat_s2(:,:,:) - 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 (dleat_s2(N_states, N_det_non_ref, 2), delta_s2(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 - n = 0 - delta(:,0,:) = 0d0 - delta(:,:nlink(J),1) = 0d0 - delta(:,:nlink(i_I),2) = 0d0 - delta_s2(:,0,:) = 0d0 - delta_s2(:,:nlink(J),1) = 0d0 - delta_s2(:,: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_cache(J,i) == 0.d0) cycle - if(h_cache(i_I,i) == 0.d0) cycle - - komon(0) += 1 - kn = komon(0) - komon(kn) = i - - do i_state = 1,N_states - dkI = h_cache(J,i) * dij(i_I, i, i_state) - dleat(i_state, kn, 1) = dkI - dleat(i_state, kn, 2) = dkI - - dkI = s2_cache(J,i) * dij(i_I, i, i_state) - dleat_s2(i_state, kn, 1) = dkI - dleat_s2(i_state, kn, 2) = dkI - end do - - end do - - komoned = .true. - end if - - integer :: hpmin(2) - hpmin(1) = 2 - HP(1,k) - hpmin(2) = 2 - HP(2,k) - - do m = 1, komon(0) - - i = komon(m) - if(HP(1,i) <= hpmin(1) .and. HP(2,i) <= hpmin(2) ) then - cycle - end if - - call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) - if(.not. ok) cycle - - do i_state = 1, N_states - contrib = dij(i_I, k, i_state) * dleat(i_state, m, 2) - contrib_s2 = dij(i_I, k, i_state) * dleat_s2(i_state, m, 2) - delta(i_state,ll,1) += contrib - delta_s2(i_state,ll,1) += contrib_s2 - 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) - delta_s2(i_state,0,1) -= contrib_s2 * ci_inv(i_state) * psi_non_ref_coef(l,i_state) - endif - - if(I_i == J) cycle - contrib = dij(J, l, i_state) * dleat(i_state, m, 1) - contrib_s2 = dij(J, l, i_state) * dleat_s2(i_state, m, 1) - delta(i_state,kk,2) += contrib - delta_s2(i_state,kk,2) += contrib_s2 - 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) - delta_s2(i_state,0,2) -= contrib_s2 * 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, delta_s2, 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, delta_s2, 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) - double precision,intent(inout) :: delta_s2(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, delta_s2(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) ! delta_s2(1,0,1) = delta_I delta_s2(1,0,2) = delta_J - if (rc /= (n(kk)+1)*8*N_states) then - print *, irp_here, 'f77_zmq_send( zmq_socket_push, delta_s2, (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, delta_s2, 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) - double precision, intent(inout) :: delta_s2(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, delta_s2(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_s2, (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, idx(1,kk), 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_,delta_ii_s2_,delta_ij_s2_) - 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) - double precision,intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref) - double precision,intent(inout) :: delta_ii_s2_(N_states,N_det_ref) - -! integer :: j,l - integer :: rc - - double precision, allocatable :: delta(:,:,:), delta_s2(:,:,:) - - 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 - delta_ii_s2_(:,:) = 0d0 - delta_ij_s2_(:,:,:) = 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), delta_s2(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, delta_s2, 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) - delta_ij_s2_(i_state,idx(l,1),i_I) += delta_s2(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) - delta_ij_s2_(i_state,idx(l,2),J) += delta_s2(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) - delta_ii_s2_(i_state,i_I) += delta_s2(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) - delta_ii_s2_(i_state,J) += delta_s2(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, delta_s2 ) - - 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) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_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,delta_ii_s2_old,delta_ij_s2_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,delta_ii_s2_old,delta_ij_s2_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/mrcc_selected/ezfio_interface.irp.f b/plugins/mrcc_selected/ezfio_interface.irp.f deleted file mode 100644 index 062af449..00000000 --- a/plugins/mrcc_selected/ezfio_interface.irp.f +++ /dev/null @@ -1,61 +0,0 @@ -! DO NOT MODIFY BY HAND -! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py -! from file /home/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg - - -BEGIN_PROVIDER [ double precision, thresh_dressed_ci ] - implicit none - BEGIN_DOC -! Threshold on the convergence of the dressed CI energy - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrcc_selected_thresh_dressed_ci(has) - if (has) then - call ezfio_get_mrcc_selected_thresh_dressed_ci(thresh_dressed_ci) - else - print *, 'mrcc_selected/thresh_dressed_ci not found in EZFIO file' - stop 1 - endif - -END_PROVIDER - -BEGIN_PROVIDER [ integer, n_it_max_dressed_ci ] - implicit none - BEGIN_DOC -! Maximum number of dressed CI iterations - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrcc_selected_n_it_max_dressed_ci(has) - if (has) then - call ezfio_get_mrcc_selected_n_it_max_dressed_ci(n_it_max_dressed_ci) - else - print *, 'mrcc_selected/n_it_max_dressed_ci not found in EZFIO file' - stop 1 - endif - -END_PROVIDER - -BEGIN_PROVIDER [ integer, lambda_type ] - implicit none - BEGIN_DOC -! lambda type - END_DOC - - logical :: has - PROVIDE ezfio_filename - - call ezfio_has_mrcc_selected_lambda_type(has) - if (has) then - call ezfio_get_mrcc_selected_lambda_type(lambda_type) - else - print *, 'mrcc_selected/lambda_type not found in EZFIO file' - stop 1 - endif - -END_PROVIDER diff --git a/plugins/mrcc_selected/mrcc_selected.irp.f b/plugins/mrcc_selected/mrcc_selected.irp.f deleted file mode 100644 index 91592e62..00000000 --- a/plugins/mrcc_selected/mrcc_selected.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -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/mrcc_selected/mrcepa0_general.irp.f b/plugins/mrcc_selected/mrcepa0_general.irp.f deleted file mode 100644 index e3a2d1f5..00000000 --- a/plugins/mrcc_selected/mrcepa0_general.irp.f +++ /dev/null @@ -1,245 +0,0 @@ - - -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) - - integer :: n_it_mrcc_max - double precision :: thresh_mrcc - double precision, allocatable :: lambda(:) - allocate (lambda(N_states)) - - - thresh_mrcc = thresh_dressed_ci - n_it_mrcc_max = n_it_max_dressed_ci - - 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_mrcepa0_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_mrcepa0_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 print_cas_coefs - implicit none - - integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, (psi_cas_coef(i,j), j=1,N_states) - call debug_det(psi_cas(1,1,i),N_int) - enddo - call write_double(6,ci_energy(1),"Initial CI energy") - -end - - - - -subroutine run_pt2_old(N_st,energy) - implicit none - integer :: i,j,k - integer, intent(in) :: N_st - double precision, intent(in) :: energy(N_st) - double precision :: pt2_redundant(N_st), pt2(N_st) - double precision :: norm_pert(N_st),H_pert_diag(N_st) - - pt2_redundant = 0.d0 - pt2 = 0d0 - !if(lambda_mrcc_pt2(0) == 0) return - - print*,'Last iteration only to compute the PT2' - - print * ,'Computing the redundant PT2 contribution' - - if (mrmode == 1) then - - N_det_generators = lambda_mrcc_kept(0) - N_det_selectors = lambda_mrcc_kept(0) - - do i=1,N_det_generators - j = lambda_mrcc_kept(i) - 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 - - else - - N_det_generators = N_det_non_ref - N_det_selectors = N_det_non_ref - - do i=1,N_det_generators - j = i - 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 - - endif - - 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 - - call H_apply_mrcepa_PT2(pt2_redundant, norm_pert, H_pert_diag, N_st) - - print * ,'Computing the remaining contribution' - - threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) - threshold_generators = max(threshold_generators,threshold_generators_pt2) - - 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 - - call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st) - - - print *, "Redundant PT2 :",pt2_redundant - print *, "Full PT2 :",pt2 - print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1) - pt2 = pt2 - pt2_redundant - - 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_mrcepa0_energy_pt2(energy(1)+pt2(1)) - -end - -subroutine run_pt2(N_st,energy) - implicit none - integer :: i,j,k - integer, intent(in) :: N_st - double precision, intent(in) :: energy(N_st) - double precision :: pt2(N_st) - double precision :: norm_pert(N_st),H_pert_diag(N_st) - - pt2 = 0d0 - !if(lambda_mrcc_pt2(0) == 0) return - - print*,'Last iteration only to compute the PT2' - - N_det_generators = N_det_cas - N_det_selectors = N_det_non_ref - - do i=1,N_det_generators - do k=1,N_int - psi_det_generators(k,1,i) = psi_ref(k,1,i) - psi_det_generators(k,2,i) = psi_ref(k,2,i) - enddo - do k=1,N_st - psi_coef_generators(i,k) = psi_ref_coef(i,k) - enddo - enddo - do i=1,N_det - do k=1,N_int - psi_selectors(k,1,i) = psi_det_sorted(k,1,i) - psi_selectors(k,2,i) = psi_det_sorted(k,2,i) - enddo - do k=1,N_st - psi_selectors_coef(i,k) = psi_coef_sorted(i,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 - - call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st) - -! call ezfio_set_full_ci_energy_pt2(energy+pt2) - - 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_mrcepa0_energy_pt2(energy(1)+pt2(1)) - -end -