From c504518542287b43e243681b7de019f79b448666 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 29 Mar 2016 23:18:26 +0200 Subject: [PATCH] MRCC acceleration --- plugins/MRCC_Utils/H_apply.irp.f | 42 +--- plugins/MRCC_Utils/mrcc_dress.irp.f | 254 +++++++++++++---------- plugins/MRCC_Utils/mrcc_utils.irp.f | 52 ++--- plugins/Psiref_Utils/psi_ref_utils.irp.f | 25 +++ src/Determinants/slater_rules.irp.f | 72 +++---- src/Integrals_Bielec/map_integrals.irp.f | 4 +- 6 files changed, 238 insertions(+), 211 deletions(-) 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_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 79e9ef7c..c41a3431 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -41,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 @@ -50,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) ] @@ -92,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 @@ -200,3 +189,4 @@ subroutine diagonalize_CI_dressed SOFT_TOUCH psi_coef 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