From 48cb3b3ddc9b6c44ed40da9792b2ded501175713 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Fri, 1 Apr 2016 12:00:03 +0200 Subject: [PATCH] different binaries for mrcepa0/mrsc2/mrsc2sub + corrected reversed index in mrcc_utils --- plugins/MRCC_Utils/davidson.irp.f | 4 +- plugins/MRCC_Utils/mrcc_general.irp.f | 4 +- plugins/MRCC_Utils/mrcc_utils.irp.f | 1 - plugins/mrcepa0/dressing.irp.f | 426 +++++++++++++++----------- plugins/mrcepa0/mrcepa0.irp.f | 2 + plugins/mrcepa0/mrcepa0_general.irp.f | 39 +-- 6 files changed, 275 insertions(+), 201 deletions(-) diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index d278ba13..c9dec40a 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -486,8 +486,8 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) i = idx_ref(ii) do jj = 1, n_det_non_ref j = idx_non_ref(jj) - vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j) - vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i) + vt (i) = vt (i) + delta_ij(istate,jj,ii)*u_0(j) + vt (j) = vt (j) + delta_ij(istate,jj,ii)*u_0(i) enddo enddo !$OMP END DO diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index 50343fdb..647caa63 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -35,8 +35,8 @@ subroutine mrcc_iterations ! lambda = min(1.d0, lambda * 1.1d0) ! endif ! print *, 'energy lambda ', lambda - E_past(j) = E_new - j +=1 +! E_past(j) = E_new +! j +=1 call save_wavefunction if (iteration > 200) then exit diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 80fdd4c7..ec4b5bf8 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -191,4 +191,3 @@ subroutine diagonalize_CI_dressed(lambda) SOFT_TOUCH psi_coef end - diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 30b161d4..4014e789 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -2,19 +2,47 @@ use bitmasks - 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, delta_ij, (N_states,N_det_non_ref,N_det_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ] use bitmasks implicit none + integer :: i, j, i_state -! delta_ij(:,:,:) = delta_ij_old(:,:,:) -! delta_ii(:,:) = delta_ii_old(:,:) - delta_ij(:,:,:) = delta_mrcepa0_ij(:,:,:)! - delta_sub_ij(:,:,:) - delta_ii(:,:)= delta_mrcepa0_ii(:,:)! - delta_sub_ii(:,:) + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + + do i_state = 1, N_states + + if(mrmode == 3) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) + end do + end do + else if(mrmode == 2) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_old(i,i_state) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_old(i,j,i_state) + end do + end do + else if(mrmode == 1) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) + end do + end do + else + stop "invalid mrmode" + end if + end do END_PROVIDER + BEGIN_PROVIDER [ integer, 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_det_non_ref) ] @@ -27,7 +55,6 @@ END_PROVIDER integer i, II, j, k logical, external :: detEq - print *, "provide cepa0" active_sorb(:) = 0_8 nonactive_sorb(:) = not(0_8) @@ -49,10 +76,6 @@ END_PROVIDER call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int) -! do i=1, N_det_non_ref -! print "(B30,B30)", det_noactive(1,1,i), det_noactive(1,2,i) -! end do -! stop do i=1,N_det_ref det_ref_active(i) = iand(psi_ref(1,1,i), active_sorb(1)) det_ref_active(i) = det_ref_active(i) + iand(psi_ref(1,2,i), active_sorb(2)) * 2_8**32_8 @@ -83,17 +106,25 @@ END_PROVIDER use bitmasks implicit none integer :: i,j,k - double precision :: Hjk, Hki, Hij, mat(2,2) - integer i_state + double precision :: Hjk, Hki, Hij + integer i_state, degree provide lambda_mrcc i_state = 1 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) + if(degree /= 2 .and. degree /= 0) cycle delta_cas(i,j,i_state) = 0d0 do k=1,N_det_non_ref +! call get_excitation_degree(psi_ref(1,1,j), psi_non_ref(1,1,k), degree, N_int) +! if(degree /= 2) cycle +! call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,k), degree, N_int) +! if(degree /= 2) cycle + call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) + delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k) end do delta_cas(j,i,i_state) = delta_cas(i,j,i_state) @@ -101,6 +132,9 @@ END_PROVIDER end do END_PROVIDER +!-199.0906497310625 +!-199.0913388716010 + logical function detEq(a,b,Nint) use bitmasks implicit none @@ -120,8 +154,8 @@ 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_old, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_old, (N_det_ref,N_states) ] use bitmasks implicit none @@ -135,11 +169,14 @@ end function integer, allocatable :: idx_sorted_bit(:) integer, external :: get_index_in_psi_det_sorted_bit logical, external :: is_in_wavefunction + double precision, allocatable :: hab(:,:) integer :: II, blok provide det_cepa0_active delta_cas lambda_mrcc + + if(N_int /= 1) then print *, "mrcepa0 experimental N_int==1" stop @@ -148,7 +185,8 @@ end function i_state = 1 delta_mrcepa0_ii(:,:) = 0d0 delta_mrcepa0_ij(:,:,:) = 0d0 - + !allocate(hab(N_det_non_ref, N_det_non_ref)) +!hab(:,:) = 0d0 ! do i=1,N_det_ref ! delta_ii(i,i_state) = delta_cas(i,i,i_state) ! end do @@ -168,7 +206,7 @@ end function made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int) - if (degree > 2 .or. popcnt(made_hole) * popcnt(made_particle) /= degree*2) cycle + if (degree > 2 .or. popcnt(made_hole) + popcnt(made_particle) == 7650) cycle do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP" @@ -185,19 +223,24 @@ end function !!!!! call get_excitation_degree(psi_ref(1,1,J),psi_non_ref(1,1,det_cepa0_idx(k)),degree,N_int) if(degree > 2) stop "BBBB" - !!!!!!!!! -! if(i/=k .and. popcnt(made_hole) /= popcnt(made_particle)) then -! print *, "=================", made_hole, made_particle -! call debug_det(psi_ref(1,1,II),N_int) -! call debug_det(psi_non_ref(1,1,det_cepa0_idx(i)),N_int) -! call debug_det(psi_ref(1,1,J),N_int) -! call debug_det(psi_non_ref(1,1,det_cepa0_idx(k)),N_int) -! print *, "=================" -! end if + call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,II),N_int,Hki) + !call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,Hki) - contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state) + !contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)! * psi_ref_coef(II, I_state)*psi_ref_coef(J, I_state)/(psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2) + contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)! * psi_ref_coef(II, I_state)*psi_ref_coef(J, I_state)/(psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2) + + + ! (psi_ref_coef(II, I_state) * psi_ref_coef(J, I_state)) / (psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2) +! if(hab(det_cepa0_idx(k), det_cepa0_idx(i)) /= 0) then +! print *, "HAB ", contrib, hab(det_cepa0_idx(k), det_cepa0_idx(i)) +! !contrib = 0d0 +! !stop +! else +! hab(det_cepa0_idx(k), det_cepa0_idx(i)) = contrib +! hab(det_cepa0_idx(i), det_cepa0_idx(k)) = contrib +! end if delta_mrcepa0_ij(II, det_cepa0_idx(i), i_state) += contrib ! if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then @@ -222,6 +265,152 @@ END_PROVIDER + BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) + double precision :: contrib, HIIi, HJk + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + + integer :: II, blok + + provide det_cepa0_active delta_cas lambda_mrcc + + + + if(N_int /= 1) then + print *, "mrcepa0 experimental N_int==1" + stop + end if + + i_state = 1 + delta_mrcepa0_ii(:,:) = 0d0 + delta_mrcepa0_ij(:,:,:) = 0d0 + + + provide mo_bielec_integrals_in_map + 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 + !- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_mrcepa0_ij, delta_mrcepa0_ii) + 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 + + made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) + made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) + + + + do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + !if(i==k) cycle + if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP" + if(iand(made_hole, det_cepa0_active(k)) /= 0) cycle + if(iand(made_particle, det_cepa0_active(k)) /= made_particle) cycle + myActive = xor(det_cepa0_active(k), made_hole) + myActive = xor(myActive, made_particle) + if(i==k .and. myActive /= det_ref_active(II)) stop "AAAA" + do J=1,N_det_ref + if(det_ref_active(J) /= myActive) cycle + call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) + contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k)) + delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib + + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state) + end if + end do + end do + end do + end do + end do + + deallocate(idx_sorted_bit) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_exp, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_exp, (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_ + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) + double precision :: contrib, HIIi, HJk + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + + integer :: II, blok + + provide det_cepa0_active delta_cas lambda_mrcc + + + print *, "mrcepa0 experimental" + if(N_int /= 1) then + print *, "mrcepa0 experimental N_int==1" + stop + end if + + i_state = 1 + delta_mrcepa0_ii_exp(:,:) = 0d0 + delta_mrcepa0_ij_exp(:,:,:) = 0d0 + + + provide mo_bielec_integrals_in_map + 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 + !- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_mrcepa0_ii_exp, delta_mrcepa0_ij_exp) + do blok=1,cepa0_shortcut(0) + do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do II=1,N_det_ref + made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) + made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) + do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do J=1,N_det_ref + if(made_hole /= iand(det_ref_active(J), xor(det_cepa0_active(k), det_ref_active(J)))) cycle + if(made_particle /= iand(det_cepa0_active(k), xor(det_cepa0_active(k), det_ref_active(J)))) cycle + call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) + contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k)) + delta_mrcepa0_ij_exp(J, det_cepa0_idx(i), i_state) += contrib + + if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + delta_mrcepa0_ii_exp(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state) + end if + end do + end do + end do + end do + end do + + deallocate(idx_sorted_bit) +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 @@ -260,7 +449,7 @@ END_PROVIDER !$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_sub_ij, delta_sub_ii) do i=1,N_det_non_ref - if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪" + if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪ (sub)" 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 @@ -279,6 +468,8 @@ END_PROVIDER 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 + 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 @@ -295,42 +486,16 @@ END_PROVIDER if(ok) cycle -! call decode_exc(exc_Ii,degree,h1_,p1_,h2_,p2_,s1_,s2_) -! call decode_exc(exc_Ik,degree2,h1,p1,h2,p2,s1,s2) -! -! -! det_tmp(:,:) = 0_bit_kind -! call set_det_bit(det_tmp, p1, s1) -! call set_det_bit(det_tmp, h1, s1) -! call set_det_bit(det_tmp, p1_, s1_) -! call set_det_bit(det_tmp, h1_, s1_) -! if(degree == 2) then -! call set_det_bit(det_tmp, p2_, s2_) -! call set_det_bit(det_tmp, h2_, s2_) -! end if -! if(degree2 == 2) then -! call set_det_bit(det_tmp, p2, s2) -! call set_det_bit(det_tmp, h2, s2) -! end if -! deg = 0 -! do ni = 1, N_int -! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) -! end do -! if(deg == 2*degree2 + 2*degree) cycle - - - - - -! if(degree == -1) 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) 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) +! contrib = delta_IJk * HIl * lambda_mrcc(i_state, l) + contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) !$OMP CRITICAL delta_sub_ij(II, i, i_state) += contrib if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then @@ -376,6 +541,7 @@ implicit none delta_ii_old(:,:) = 0 delta_ij_old(:,:,:) = 0 + i_state = 1 provide mo_bielec_integrals_in_map allocate(idx_sorted_bit(N_det)) @@ -387,7 +553,8 @@ implicit none !$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij_old, delta_ii_old) do i = 1 , N_det_non_ref - if(mod(i,1000) == 0) print *, i, N_det_non_ref + if(lambda_mrcc(i_state, i) == 0d0) cycle + if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪ (old)" do i_I = 1 , N_det_ref call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int) if(degree2 == -1) cycle @@ -396,12 +563,14 @@ implicit none call decode_exc(exc_iI,degree2,h1,p1,h2,p2,s1,s2) call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi) - diI = hIi * lambda_mrcc(i_state,i) + diI = hIi * lambda_mrcc(i_state, i) do J = 1 , N_det_ref !!! call get_excitation(psi_ref(1,1,i_I),psi_ref(1,1,J),exc_IJ,degree,phase_IJ,N_int) call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi) delta_JI = hJi * diI do k = 1 , N_det_non_ref + if(lambda_mrcc(i_state, k) == 0d0) cycle + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int) if(degree == -1) cycle @@ -411,11 +580,6 @@ implicit none det_tmp(:,:) = 0_bit_kind det_tmp2(:,:) = 0_bit_kind - !!!!!!!!!!!!!!! - - - - det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1))) det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1))) @@ -427,60 +591,6 @@ implicit none if(.not. ok) cycle !if(ok) cycle - !!!!!!!!!!!!!! - - - -! call set_det_bit(det_tmp, p1, s1) -! -! call set_det_bit(det_tmp, p1_, s1_) -! -! if(degree == 2) then -! call set_det_bit(det_tmp, p2_, s2_) -! -! end if -! if(degree2 == 2) then -! call set_det_bit(det_tmp, p2, s2) -! end if -! -! x(:) = 0 -! do ni=1,N_int -! x(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1))) -! x(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1))) -! end do -! -! -! !det_tmp(:,:) = 0_bit_kind -! -! call set_det_bit(det_tmp, h1, s1) -! call set_det_bit(det_tmp, h1_, s1_) -! if(degree == 2) then -! call set_det_bit(det_tmp, h2_, s2_) -! end if -! if(degree2 == 2) then -! call set_det_bit(det_tmp, h2, s2) -! end if -! -! y(1) = -x(1) -! y(2) = -x(2) -! do ni=1,N_int -! y(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1))) -! y(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1))) -! end do -! -! ! print *, x, y -! -! if(x(1) * y(1) /= 0) cycle -! if(x(2) * y(2) /= 0) cycle -! -! -! -! deg = 0 -! do ni = 1, N_int -! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) -! end do -! if(deg /= 2*degree2 + 2*degree) cycle - call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) @@ -500,17 +610,22 @@ implicit none l = get_index_in_psi_det_sorted_bit(det_tmp, N_int) if(l == 0) cycle l = idx_sorted_bit(get_index_in_psi_det_sorted_bit(det_tmp, N_int)) - if(l ==-1) cycle + call i_h_j(psi_non_ref(1,1,k), psi_ref(1,1,i_I),N_int,HkI) - dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik - + !dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik * Xref(I_i) + dkI(i_state) = HkI * lambda_mrcc(i_state, k) * phase_Jl * phase_Ik + + !!!! + call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int) + if(degree /= 2 .and. degree /= 0) cycle + !!!!!! + !$OMP CRITICAL contrib = dkI(i_state) * delta_JI - !erro += abs(dkI(i_state) - psi_non_ref_coef(k,i_state) / psi_ref_coef(1,i_state)) delta_ij_old(i_I,l,i_state) += contrib if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(k,i_state) @@ -521,61 +636,24 @@ implicit none enddo enddo !$OMP END PARALLEL DO + +! double precision :: error, acc +! integer :: II +! error = 0d0 +! do i=1, N_det_non_ref +! acc = 0d0 +! do II=1, N_det_ref +! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), N_int, HIi) +! acc += HIi * lambda_mrcc(i_state, i) * Xref(II) * psi_ref_coef(II, i_state) +! end do +! error += (psi_non_ref_coef(i, i_state) - acc)**2 +! end do +! print *, "QUALITY ", error + + deallocate(idx_sorted_bit) END_PROVIDER - - -subroutine apply_excitation(det, exc, res, ok, Nint) - use bitmasks - implicit none - - integer, intent(in) :: Nint - integer, intent(in) :: exc(0:2,2,2) - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: h1,p1,h2,p2,s1,s2,degree - integer :: ii, pos - - - ok = .false. - degree = exc(0,1,1) + exc(0,1,2) - if(.not. (degree > 0 .and. degree <= 2)) then - print *, "apply ex" - STOP - endif - - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - res = det - - ii = (h1-1)/bit_kind_size + 1 - pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s1) = ibclr(res(ii, s1), pos) - - ii = (p1-1)/bit_kind_size + 1 - pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s1) = ibset(res(ii, s1), pos) - - if(degree == 2) then - ii = (h2-1)/bit_kind_size + 1 - pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s2) = ibclr(res(ii, s2), pos) - - ii = (p2-1)/bit_kind_size + 1 - pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s2) = ibset(res(ii, s2), pos) - endif - - ok = .true. -end subroutine - - - diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f index b9eb2fc5..7877cdda 100644 --- a/plugins/mrcepa0/mrcepa0.irp.f +++ b/plugins/mrcepa0/mrcepa0.irp.f @@ -1,5 +1,7 @@ program mrcepa0 implicit none + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + mrmode = 1 if (.not.read_wf) then print *, 'read_wf has to be true.' stop 1 diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index 87ae25f3..b307ca85 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -4,6 +4,10 @@ subroutine run_mrcepa0 call mrcepa0_iterations end +BEGIN_PROVIDER [ integer, mrmode ] + +END_PROVIDER + subroutine mrcepa0_iterations implicit none @@ -11,12 +15,13 @@ subroutine mrcepa0_iterations double precision :: E_new, E_old, delta_e integer :: iteration,i_oscillations - double precision :: E_past(4) + double precision :: E_past(4), lambda E_new = 0.d0 delta_E = 1.d0 iteration = 0 j = 1 i_oscillations = 0 + lambda = 1.d0 do while (delta_E > 1.d-7) iteration += 1 print *, '===========================' @@ -25,29 +30,19 @@ subroutine mrcepa0_iterations print *, '' E_old = sum(ci_energy_dressed) call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy") - call diagonalize_ci_dressed + call diagonalize_ci_dressed(lambda) 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 +! if (E_new > E_old) then +! lambda = lambda * 0.7d0 +! else +! lambda = min(1.d0, lambda * 1.1d0) +! endif +! print *, 'energy lambda ', lambda +! E_past(j) = E_new +! j +=1 call save_wavefunction -! if (i_oscillations > 5) then -! exit -! endif - if (iteration > 100) then + if (iteration > 200) then exit endif print*,'------------' @@ -55,7 +50,7 @@ subroutine mrcepa0_iterations do i = 1, N_det_ref print*,'' print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) - print*,'delta_ii(i,1) = ',delta_cas(i,i,1) + print*,'delta_ii(i,1) = ',delta_ii(i,1) enddo print*,'------------' enddo