diff --git a/config/ifort.cfg b/config/ifort.cfg index cc848cba..2b2fe0a2 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -31,7 +31,7 @@ OPENMP : 1 ; Append OpenMP flags # -ftz : Flushes denormal results to zero # [OPT] -FCFLAGS : -xHost -O2 -ip -ftz -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g # Profiling flags ################# diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 1eb2d45a..3dc1e0f0 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -2,20 +2,20 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply_zmq("FCI") +s = H_apply("FCI") s.set_selection_pt2("epstein_nesbet_2x2") -s.unset_openmp() +#s.unset_openmp() print s s = H_apply_zmq("FCI_PT2") s.set_perturbation("epstein_nesbet_2x2") -s.unset_openmp() +#s.unset_openmp() print s -s = H_apply_zmq("FCI_no_skip") +s = H_apply("FCI_no_skip") s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() -s.unset_openmp() +#s.unset_openmp() print s s = H_apply("FCI_mono") diff --git a/plugins/MRCC_CASSD/EZFIO.cfg b/plugins/MRCC_CASSD/EZFIO.cfg index 21cc5b98..e145c9e0 100644 --- a/plugins/MRCC_CASSD/EZFIO.cfg +++ b/plugins/MRCC_CASSD/EZFIO.cfg @@ -2,3 +2,16 @@ type: double precision doc: Calculated energy interface: ezfio + +[thresh_mrcc] +type: Threshold +doc: Threshold on the convergence of the MRCC energy +interface: ezfio,provider,ocaml +default: 1.e-7 + +[n_it_mrcc_max] +type: Strictly_positive_int +doc: Maximum number of MRCC iterations +interface: ezfio,provider,ocaml +default: 20 + diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 1cafc8de..decc5a75 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -3,19 +3,19 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * s = H_apply("mrcc") -s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" +s.data["parameters"] = ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" s.data["declarations"] += """ - integer, intent(in) :: Ndet_ref,Ndet_non_ref - double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(in) :: delta_ii_(Ndet_ref,*) + integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref + double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, Ndet_ref) + double precision, intent(in) :: delta_ii_(Nstates, Ndet_ref) """ -s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" -s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" -s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" +s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" s.data["decls_main"] += """ - integer, intent(in) :: Ndet_ref,Ndet_non_ref - double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(in) :: delta_ii_(Ndet_ref,*) + integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates + double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref) + double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref) """ s.data["finalization"] = "" s.data["copy_buffer"] = "" @@ -24,27 +24,5 @@ s.data["size_max"] = "3072" print s -s = H_apply("mrcepa") -s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" -s.data["declarations"] += """ - integer, intent(in) :: Ndet_ref,Ndet_non_ref - double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(in) :: delta_ii_(Ndet_ref,*) -""" -s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" -s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" -s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" -s.data["decls_main"] += """ - integer, intent(in) :: Ndet_ref,Ndet_non_ref - double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(in) :: delta_ii_(Ndet_ref,*) -""" -s.data["finalization"] = "" -s.data["copy_buffer"] = "" -s.data["generate_psi_guess"] = "" -s.data["size_max"] = "3072" -# print s - - END_SHELL diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 5747b174..c3f3debf 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -14,14 +14,14 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ] END_PROVIDER -subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) +subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_ref, Ndet_non_ref - double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(inout) :: delta_ii_(Ndet_ref,*) + integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref + double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref) + double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l @@ -32,10 +32,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n integer(bit_kind) :: tq(Nint,2,n_selected) integer :: N_tq, c_ref ,degree - double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision :: hIk, hla, hIl, dIk(Nstates), dka(Nstates), dIa(Nstates) double precision, allocatable :: dIa_hla(:,:) double precision :: haj, phase, phase2 - double precision :: f(N_states), ci_inv(N_states) + double precision :: f(Nstates), ci_inv(Nstates) integer :: exc(0:2,2,2) integer :: h1,h2,p1,p2,s1,s2 integer(bit_kind) :: tmp_det(Nint,2) @@ -46,10 +46,11 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n integer(bit_kind),intent(in) :: key_mask(Nint, 2) integer,allocatable :: idx_miniList(:) integer :: N_miniList, ni, leng + double precision, allocatable :: hij_cache(:) leng = max(N_det_generators, N_det_non_ref) - allocate(miniList(Nint, 2, leng), idx_miniList(leng)) + allocate(miniList(Nint, 2, leng), idx_miniList(leng), hij_cache(N_det_non_ref)) !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) @@ -60,124 +61,157 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) - - allocate (dIa_hla(N_states,Ndet_non_ref)) - + + allocate (dIa_hla(Nstates,Ndet_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) - end if + if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint) + end if do i_alpha=1,N_tq -! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) + ! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) 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 - + + 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) - 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 - ! - ! - 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) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) - 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 - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) - - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree_alpha(k_sd) == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - endif - - ! - 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 - call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) - 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) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 - enddo - 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) = 1.d0/psi_ref_coef(i_I,i_state) - enddo - 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,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 l_sd=1,idx_alpha(0) - k_sd = idx_alpha(l_sd) - do i_state=1,N_states - delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd) - if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then - delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) - else - delta_ii_(i_I,i_state) = 0.d0 - endif - enddo - enddo - call omp_unset_lock( psi_ref_lock(i_I) ) + ! 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,Nstates + 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,Nstates + 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,Nstates + dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + 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 + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree_alpha(k_sd) == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + ! + do i_state=1,Nstates + 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,Nstates + if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then + loop = .False. + exit + endif + enddo + 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,Nstates + dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + enddo + endif + exit + endif + enddo + do i_state=1,Nstates + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,Nstates + ci_inv(i_state) = 1.d0/psi_ref_coef(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,Nstates + dIa_hla(i_state,k_sd) = dIa(i_state) * hla + enddo + enddo + call omp_set_lock( psi_ref_lock(i_I) ) + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + do i_state=1,Nstates + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + 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) + else + delta_ii_(i_state,i_I) = 0.d0 + endif + enddo + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo - deallocate (dIa_hla) + deallocate (dIa_hla,hij_cache) deallocate(miniList, idx_miniList) end diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index c567c76a..5d9acfc1 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -31,23 +31,7 @@ subroutine mrcc_iterations E_past(j) = E_new j +=1 - if(j>4)then - j=1 - endif - if(iteration > 4) then - if(delta_E > 1.d-10)then - if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then - print*,'OSCILLATIONS !!!' - oscillations = .True. - i_oscillations +=1 - lambda_mrcc_tmp = lambda_mrcc - endif - endif - endif call save_wavefunction -! if (i_oscillations > 5) then -! exit -! endif if (iteration > 200) then exit endif diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 1e2f974d..c41a3431 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,102 +1,35 @@ - BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ] - END_PROVIDER - - - BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] -&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ] +BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] implicit none BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) END_DOC - integer :: i,k,j - double precision :: ihpsi(N_states), hii,delta_e_eff,ihpsi_current(N_states),hij - integer :: i_ok,i_pert,i_pert_count - i_ok = 0 - - double precision :: phase_restart(N_states),tmp - do k = 1, N_states - phase_restart(k) = dsign(1.d0,psi_ref_coef_restart(1,k)/psi_ref_coef(1,k)) - enddo - i_pert_count = 0 + integer :: i,k + double precision :: ihpsi(N_states),ihpsi_current(N_states) + integer :: i_pert_count - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref_restart, psi_ref_coef_restart, N_int, N_det_ref,& - size(psi_ref_coef_restart,1), n_states, ihpsi) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) -! TODO --- Test perturbatif ------ - do k=1,N_states - lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - ! TODO : i_h_psi peut sortir de la boucle? - 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) - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 - endif - tmp = psi_non_ref_coef(i,k)/ihpsi_current(k) - i_pert = 0 - ! Perturbation only if 1st order < 0.5 x second order - if((ihpsi(k) * lambda_pert(k,i)) < 0.5d0 * psi_non_ref_coef_restart(i,k) )then - i_pert = 1 - else - do j = 1, N_det_ref - call i_H_j(psi_non_ref(1,1,i),psi_ref(1,1,j),N_int,hij) - ! Perturbation diverges when hij*tmp > 0.5 - if(dabs(hij * tmp).ge.0.5d0)then - i_pert_count +=1 - i_pert = 1 - exit - endif - enddo - endif - if( i_pert == 1)then - pert_determinants(k,i) = i_pert - endif - if(pert_determinants(k,i) == 1)then - i_ok +=1 - lambda_mrcc(k,i) = lambda_pert(k,i) - else - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - endif - enddo -! TODO --- Fin test perturbatif ------ + i_pert_count = 0 + lambda_mrcc = 0.d0 + + 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) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif + if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then + i_pert_count +=1 + else + lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) + endif + enddo enddo -!if(oscillations)then -! print*,'AVERAGING the lambda_mrcc with those of the previous iterations' -! do i = 1, N_det_non_ref -! do k = 1, N_states -! double precision :: tmp -! tmp = lambda_mrcc(k,i) -! lambda_mrcc(k,i) += lambda_mrcc_tmp(k,i) -! lambda_mrcc(k,i) = lambda_mrcc(k,i) * 0.5d0 -! if(dabs(tmp - lambda_mrcc(k,i)).ge.1.d-9)then -! print*,'' -! print*,'i = ',i -! print*,'psi_non_ref_coef(i,k) = ',psi_non_ref_coef(i,k) -! print*,'lambda_mrcc(k,i) = ',lambda_mrcc(k,i) -! print*,' tmp = ',tmp -! endif -! enddo -! enddo -!endif print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of Perturbatively treated determinants = ',i_ok - print*,'i_pert_count = ',i_pert_count + print*,'Number of ignored determinants = ',i_pert_count print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) - END_PROVIDER -BEGIN_PROVIDER [ double precision, lambda_mrcc_tmp, (N_states,psi_det_size) ] - implicit none - lambda_mrcc_tmp = 0.d0 -END_PROVIDER - -BEGIN_PROVIDER [ logical, oscillations ] - implicit none - oscillations = .False. -END_PROVIDER - - !BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] @@ -108,8 +41,22 @@ END_PROVIDER !call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) !END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] +BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] + implicit none + BEGIN_DOC + ! < ref | H | Non-ref > matrix + END_DOC + integer :: i_I, k_sd + do i_I=1,N_det_ref + do k_sd=1,N_det_non_ref + call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,k_sd),N_int,hij_mrcc(k_sd,i_I)) + enddo + enddo + +END_PROVIDER + + 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) ] implicit none BEGIN_DOC ! Dressing matrix in N_det basis @@ -117,32 +64,7 @@ END_PROVIDER integer :: i,j,m delta_ij = 0.d0 delta_ii = 0.d0 - call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref) - double precision :: max_delta - double precision :: accu - integer :: imax,jmax - max_delta = 0.d0 - accu = 0.d0 - do i = 1, N_det_ref - do j = 1, N_det_non_ref - accu += psi_non_ref_coef(j,1) * psi_ref_coef(i,1) * delta_ij(i,j,1) - if(dabs(delta_ij(i,j,1)).gt.max_delta)then - max_delta = dabs(delta_ij(i,j,1)) - imax = i - jmax = j - endif - enddo - enddo - print*,'' - print*,'' - print*,' = ',accu - print*,'MAX VAL OF DRESING = ',delta_ij(imax,jmax,1) - print*,'imax,jmax = ',imax,jmax - print*,'psi_ref_coef(imax,1) = ',psi_ref_coef(imax,1) - print*,'psi_non_ref_coef(jmax,1) = ',psi_non_ref_coef(jmax,1) - do i = 1, N_det_ref - print*,'delta_ii(i,1) = ',delta_ii(i,1) - enddo + 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) ] @@ -159,11 +81,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] enddo do ii = 1, N_det_ref i =idx_ref(ii) - h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) + h_matrix_dressed(i,i,istate) += delta_ii(istate,ii) do jj = 1, N_det_non_ref j =idx_non_ref(jj) - h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate) - h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate) + h_matrix_dressed(i,j,istate) += delta_ij(istate,jj,ii) + h_matrix_dressed(j,i,istate) += delta_ij(istate,jj,ii) enddo enddo enddo @@ -267,3 +189,4 @@ subroutine diagonalize_CI_dressed SOFT_TOUCH psi_coef end + diff --git a/plugins/MRCC_Utils/mrcepa_dress.irp.f b/plugins/MRCC_Utils/mrcepa_dress.irp.f deleted file mode 100644 index 9789b818..00000000 --- a/plugins/MRCC_Utils/mrcepa_dress.irp.f +++ /dev/null @@ -1,260 +0,0 @@ -use omp_lib -use bitmasks - -subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) - use bitmasks - implicit none - - integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_ref, Ndet_non_ref - double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) - double precision, intent(inout) :: delta_ii_(Ndet_ref,*) - - integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) - integer :: i,j,k,l - integer :: degree_alpha(psi_det_size) - integer :: idx_alpha(0:psi_det_size) - logical :: good, fullMatch - - integer(bit_kind) :: tq(Nint,2,n_selected) - integer :: N_tq, c_ref ,degree - - double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) - double precision, allocatable :: dIa_hia(:,:) - 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(bit_kind) :: tmp_det_0(Nint,2) - integer :: iint, ipos - integer :: i_state, i_sd, 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 - integer(bit_kind) :: isum - - double precision :: hia - integer, allocatable :: index_sorted(:) - - - leng = max(N_det_generators, N_det_non_ref) - allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det)) - - !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 - - - call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) - - allocate (dIa_hia(N_states,Ndet_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) - end if - - - do i_alpha=1,N_tq - ! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) - call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) - - integer, external :: get_index_in_psi_det_sorted_bit - index_sorted = huge(-1) - do j=1,idx_alpha(0) - idx_alpha(j) = idx_miniList(idx_alpha(j)) - index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j) - end do - - ! |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 - - !TODO: MR - do i_sd=1,idx_alpha(0) - call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint) - if (degree > 2) then - cycle - endif - call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - tmp_det_0 = 0_bit_kind - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) - - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) - if (degree == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) - endif - - call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia) - - ! |alpha> - do k_sd=1,idx_alpha(0) - 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 - - call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - tmp_det = 0_bit_kind - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - endif - - isum = 0_bit_kind - do iint = 1,N_int - isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) & - + iand(tmp_det(iint,2), tmp_det_0(iint,2)) - enddo - - if (isum /= 0_bit_kind) then - cycle - endif - - ! - ! - 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) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) - 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 - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) - - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - endif - - ! - do i_state=1,N_states - dka(i_state) = 0.d0 - enddo - - -! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) ) -! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint) -! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl) -! do i_state=1,N_states -! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2 -! enddo - - do l_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 - call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) - 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) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 - enddo - 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) = 1.d0/psi_ref_coef(i_I,i_state) - enddo - - k_sd = idx_alpha(i_sd) - do i_state=1,N_states - dIa_hia(i_state,k_sd) = dIa(i_state) * hia - enddo - - call omp_set_lock( psi_ref_lock(i_I) ) - do i_state=1,N_states - delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd) - - if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then - delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) - else - delta_ii_(i_I,i_state) = 0.d0 - endif - enddo - call omp_unset_lock( psi_ref_lock(i_I) ) - enddo - enddo - - enddo - deallocate (dIa_hia,index_sorted) - deallocate(miniList, idx_miniList) -end - - - - - - - - diff --git a/plugins/MRCC_Utils/mrcepa_general.irp.f b/plugins/MRCC_Utils/mrcepa_general.irp.f deleted file mode 100644 index 3479548b..00000000 --- a/plugins/MRCC_Utils/mrcepa_general.irp.f +++ /dev/null @@ -1,97 +0,0 @@ -subroutine run_mrcepa - implicit none - call set_generators_bitmasks_as_holes_and_particles - call mrcepa_iterations -end - -subroutine mrcepa_iterations - implicit none - - integer :: i,j - - double precision :: E_new, E_old, delta_e - integer :: iteration,i_oscillations - double precision :: E_past(4) - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - j = 1 - i_oscillations = 0 - do while (delta_E > 1.d-7) - iteration += 1 - print *, '===========================' - print *, 'MRCEPA Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed) - call write_double(6,ci_energy_dressed(1),"MRCEPA energy") - call diagonalize_ci_dressed - E_new = sum(ci_energy_dressed) - delta_E = dabs(E_new - E_old) - - E_past(j) = E_new - j +=1 - if(j>4)then - j=1 - endif - if(iteration > 4) then - if(delta_E > 1.d-10)then - if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then - print*,'OSCILLATIONS !!!' - oscillations = .True. - i_oscillations +=1 - lambda_mrcc_tmp = lambda_mrcc - endif - endif - endif - call save_wavefunction -! if (i_oscillations > 5) then -! exit -! endif - 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 MRCEPA 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 - integer :: i,k - do k = 1, N_generators_bitmask - do i = 1, N_int - ! Pure single part - generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha - generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha - generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta - generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta - - ! Double excitation - generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha - generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha - generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta - generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta - - generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha - generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha - generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta - generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta - - enddo - enddo - touch generators_bitmask - - - -end diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f index fb45b13d..5115ba23 100644 --- a/plugins/Psiref_Utils/psi_ref_utils.irp.f +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -14,6 +14,31 @@ use bitmasks END_PROVIDER +BEGIN_PROVIDER [ double precision, psi_ref_coef_transp, (n_states,psi_det_size) ] + implicit none + BEGIN_DOC +! Transposed psi_ref_coef + END_DOC + integer :: i,j + do j=1,N_det_ref + do i=1, n_states + psi_ref_coef_transp(i,j) = psi_ref_coef(j,i) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_non_ref_coef_transp, (n_states,psi_det_size) ] + implicit none + BEGIN_DOC +! Transposed psi_non_ref_coef + END_DOC + integer :: i,j + do j=1,N_det_non_ref + do i=1, n_states + psi_non_ref_coef_transp(i,j) = psi_non_ref_coef(j,i) + enddo + enddo +END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ] diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index ec786941..b63ae69f 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -443,7 +443,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral_schwartz + double precision :: get_mo_bielec_integral integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -468,31 +468,31 @@ subroutine i_H_j(key_i,key_j,Nint,hij) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral_schwartz( & + hij = phase*get_mo_bielec_integral( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -510,15 +510,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -537,15 +537,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -579,7 +579,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) integer,intent(out) :: exc(0:2,2,2) integer,intent(out) :: degree - double precision :: get_mo_bielec_integral_schwartz + double precision :: get_mo_bielec_integral integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -604,31 +604,31 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral_schwartz( & + hij = phase*get_mo_bielec_integral( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -646,15 +646,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -673,15 +673,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -715,7 +715,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral_schwartz + double precision :: get_mo_bielec_integral integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -742,31 +742,31 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral_schwartz( & + hij = phase*get_mo_bielec_integral( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral_schwartz( & + hij = phase*(get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral_schwartz( & + get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -784,15 +784,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -811,15 +811,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 84b08715..4041242e 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -324,9 +324,9 @@ double precision function mo_bielec_integral(i,j,k,l) ! Returns one integral in the MO basis END_DOC integer, intent(in) :: i,j,k,l - double precision :: get_mo_bielec_integral_schwartz + double precision :: get_mo_bielec_integral PROVIDE mo_bielec_integrals_in_map - mo_bielec_integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map) + mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) return end diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index de0cd1c8..1be43412 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -155,7 +155,7 @@ function run_all_1h_1p() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -0.762303253805911E+02 1.E-3 + eq $energy -76.2289109271715 1.E-3 }